Computational protocol: Gene Turnover in the Avian Globin Gene Families and Evolutionary Changes in Hemoglobin Isoform Expression

Similar protocols

Protocol publication

[…] The 52 avian genome assemblies that we analyzed were derived from multiple sources (; ; ; ; ; ). Details regarding all examined genome assemblies and scaffolds are provided in supplementary table S1, Supplementary Material online. We identified contigs containing α- and β-globin genes in the genome assembly of each species by using BLAST () to make comparisons with the α- and β-globin genes of chicken. We then annotated globin genes within the identified gene clusters of each species by using GENSCAN () in combination with BLAST2 () to make comparisons with known exon sequences, following Hoffmann et al. (; ). Due to incomplete sequence coverage in the genome assemblies of several species, there were some cases where it was not possible to ascertain the full extent of conserved synteny.To estimate rates of gene turnover in the globin gene clusters of birds and mammals, we used a stochastic birth–death model of gene family evolution (), as implemented in the program CAFE v3.1 (). As input for the program we used ultrametric trees based on well-resolved phylogenies of birds () and eutherian mammals (). The analysis was based on data for 24 bird species for which we had complete sequence coverage of the α- and β-globin gene clusters. For comparison, we used genome sequence assemblies for 22 mammal species representing each of the major eutherian lineages, as reported in . These sets of birds and eutherian mammals span a similar range of divergence times, and therefore provide a good basis for comparing lineage-specific rates of gene family evolution. [...] We used phylogenetic reconstructions to resolve orthology for all annotated globin genes, including truncated sequences. The α- and β-type globin sequences were aligned separately and, when possible, amino acid sequence alignments were based on conceptual translations of nucleotide sequence. All alignments were performed using the L-INS-i strategy from MAFFT v. 6.8 (), and phylogenies were estimated using both maximum-likelihood and Bayesian approaches. We used Treefinder (v. March 2011; []) for the maximum-likelihood analyses, and MrBayes () for Bayesian estimates of each phylogeny. In the maximum-likelihood analyses, we used the “propose model” routine of Treefinder to select the best-fitting models of nucleotide substitution, with an independent model for each codon position in analyses based on nucleotide sequences. Model selection was based on the Akaike information criterion with correction for small sample size. In the case of MrBayes, we ran six simultaneous chains for 2 × 107 generations, sampling every 2.5 × 103 generations, using default priors. A given run was considered to have reached convergence once the likelihood scores reached an asymptotic value and the average standard deviation of split frequencies remained less than 0.01. We discarded all trees that were sampled prior to convergence, and we evaluated support for the nodes and parameter estimates from a majority rule consensus of the last 2,500 trees.For the β-type globin genes, we analyzed four separate sets of sequences: 1) an inclusive data set that included all genes, apparent pseudogenes, and truncated sequences (containing at least one full-length exon) that were derived from the genome assemblies, 2) a data set restricted to genes with full-length coding sequences and intact reading frames, and 3) the set of full-length coding sequences plus 34 βA-globin cDNA sequences from avian taxa that do not currently have fully sequenced genomes. We used each of these data sets to perform likelihood-based topology tests ().For the β-type globin genes, we inferred orthologous relationships by using an approach that combined the “orthology-by-content” and “orthology-by-context” criteria (, ). We applied the orthology-by-content criterion to assign genes to the βH, βA, or ρ+ε groups (based on phylogenetic relationships), and we then used the orthology-by-context criterion to assign genes in the latter group as ρ- or ε-globin (based on positional homology). [...] To measure variation in functional constraint among the avian α- and β-type globin genes and to test for evidence of positive selection, we estimated ω (=dN/dS, the ratio of the rate of nonsynonymous substitution per nonsynonymous site [dN] to the rate of synonymous substitution per synonymous site [dS] site) using a maximum-likelihood approach () implemented in the CODEML program, v. 4.8 (). In all cases, we used LRTs to compare nested sets of models (). For each alignment, we first compared models that allow ω to vary among codons (M0 vs. M3, M1a vs. M2a). The latter LRT constitutes a test of positive selection, as the alternative model (M2a) allows a subset of sites to have ω > 1, in contrast to the null model (M1a) that only includes site classes with ω < 1. These analyses revealed statistically significant variation in ω among sites in the alignment of βA-globin sequences (), so we then used the BEB approach () to calculate the posterior probability that a given site belongs to the site class with ω > 1, which suggests a possible history of positive selection. [...] For the comparative analysis of HbA/HbD expression levels, we combined our experimental data for 68 species with published data for an additional 54 species (supplementary table S3, Supplementary Material online), yielding data for a phylogenetically diverse set of 122 species in total. We used the final supertree to conduct comparative analyses under a Brownian motion model of evolution in R v. 2.15.1. Maximum-likelihood estimates of ancestral states for percent expression of HbD were generated using APE (). Phylogenetic inertia in expression level (percent) of HbD was quantified using Blomberg’s K () and Pagel’s Lambda (, calculated using GEIGER (). We compiled body mass data for each species from Museum of Southwestern Biology specimens and Dunning (, and we tested for an association between body mass and expression level of HbD using a phylogenetic generalized least-squares model as implemented in APE (). We log-transformed body mass and arcsin-squareroot-transformed HbD percent prior to all analyses. […]

Pipeline specifications

Software tools GENSCAN, CAFE, MAFFT, MrBayes, PAML, APE, GEIGER
Application Phylogenetics