Computational protocol: The complexity of selection at the major primate β-defensin locus

Similar protocols

Protocol publication

[…] In alignments and figures primate species names were abbreviated to two letters as follows: Cercopithecus preussi (Preuss's monkey) (cp), Cercopithecus aethiops (Vervet monkey) (ca), Cercopithecus erythrogaster (Red-bellied monkey) (ce), Presbytis cristata (Silvered langur) (pc), Presbytis obscurus (Spectacled langur) (po), Presbytis melalophos (Banded langur) (pm), Macaca mulatta (Rhesus Macaque) (mm), Macaca fascicularis (crab-eating macaque) (mf), Papio anubis (olive baboon) (pa), Hylobates lar (Lar gibbon) (hl), Hylobates moloch (Silvery gibbon) (hm), Hylobates concolor (crested gibbon) (hc), Callithrix jacchus (common marmoset) (cj), Saguinus oedipus (cotton-top tamarin) (so), Pan troglodytes (chimpanzee) (pt), Gorilla gorilla (gorilla) (gg), Pongo pygmaeus (orangutan) (pp), Homo sapiens (human) (hs). The Cercopithecidae are represented by cp, ca, ce, pc, po, pm, mm, mf and pa; the Hylobatidae by hl, hm and hc; the Callitrichidae by cj and so; and the Hominidae by pt, gg, pp and hs. Note that sequences from every species were not available for each primate gene (see Figure and Additional Files).All mouse sequences were recently published by Zaballos et al. []. The H. sapiens and P. anubis sequences were as previously published []. Previously published sequence data for primate DEFB1 [], DEFB4 [] and DEFB103 [] were combined with the following novel data. Whole genome shotgun reads from the M. mulatta genome representing DEFB104 (69840222, 73807150), DEFB105 (73492526, 74381588, 72060044) and DEFB108 (71259620, 72564100, 71889652, 72776644, RHQRA66TR) were identified using BLAST [] from the Ensembl Trace server . Genomic sequence assembly contigs from the P. troglodytes genome were obtained from Ensembl in the same way for DEFB105 (AADA01159356) and DEFB107 (AADA01159356). The published Rattus norvegicus (rat) genome assembly [] and the full Canis familiaris (dog) ~7.6X coverage whole genome shotgun data (downloaded September 2004 from the Ensembl Trace Server: ) were searched using TBLASTN [] with default settings.PCR of novel second exon sequences from P. anubis, C. aethiops, P. pygmaeus,C. jacchus DNA was achieved using primers designed to the human exon 2 flanking sequence. PCR programmes were used with a relaxed annealing temperature that revealed a single species by gel electrophoresis. PCR products were cloned and several clones were analysed for each PCR. At least two clones were sequenced in both directions. Primers were as follows with forward primer sequence preceding the reverse primer sequence for each gene. DEFB103: 5'GTGCTGTTTTGTCATTGCAG, 5'GATTTAAAAAAAAAAATCAAGCTC; DEFB104: 5'CAGTGCCATATCCTGTTATCTAG, 5'GCTGCTAGGCCGCAGGAAGG; DEFB105: 5'GCAGCTCTTTCTTGGCAGAG, 5'GCTGGTCTGGTTTGTCAGATC; DEFB106: 5'TGGCTCCTTCCCTGTGTAG, 5'CACTTGACAAACTGAGCAAAG; DEFB107: 5'CTGCTTTCTTTACTTAGCCA, 5'GTGCTTAGTTTTTAATGTTTCTTTC; DEFB108: 5'CAATAACCCCTTCTGCATGTAG, 5'CTCAATTCTTGGTTGATGCCC. Novel sequences were deposited in GenBank under accession numbers: AY831729, AY831730, AY831731, AY831732, AY831733, AY831734, AY831735, AY831736, AY831737, AY831738, AY831739, AY831740, AY831741, AY831742, AY831743, AY831744, AY831745, AY831746. [...] Protein sequences were aligned using T-Coffee []. An alignment of coding sequences corresponding to the protein alignment was constructed using the tranalign program from the EMBOSS package []. A phylogenetic tree was constructed using MEGA (version 2.1 []) by the neighbour joining (NJ) method [] using p-distance estimates, which is thought to be the most reliable method for constructing NJ trees of closely related sequences []. The tree was rooted with chicken gallinacin 1 (P46156) and the reliability of each node was assessed with 1000 bootstrap replications. All of the best supported branches (>= 50% of replications) were also observed in an equivalent NJ tree constructed with the gamma distribution model implemented to account for rate heterogeneity among sites. The shape parameter of the gamma distribution (α) was estimated using baseml from the PAML package (version 3.13 []). The same alignment was used as the basis for trees constructed by maximum parsimony (using MEGA version 2.1 with default settings) and maximum likelihood (using PHYLIP version 3.6 [] with default settings). In both cases the trees produced shared substantively similar topology with the NJ trees. (All trees are available on demand.) Nodes within the primate lineage were dated according to a widely accepted phylogenetic analysis []. The divergence of rat and mouse was taken to be 12–24 MYA [] and the last common ancestor of mammals was assumed to be 92 MYA [].Likelihood ratio tests (LRTs) were performed using codeml with the site-specific models of the PAML package and the tree constructed as above. The six site-specific models recommended by Anisimova et al. [] were tested: M0 (one-ratio), M1 (neutral), M2 (selection), M3 (discrete), M7 (beta), and M8 (beta+ω). These LRTs indicate whether the substitutions inferred from an alignment are best explained by one of two models of ω = dN/dS, where dN and dS are the nonsynonymous and synonymous substitution rates respectively. When parameter estimates under a model allowing positive selection suggest the presence of a number of sites with ω > 1, Bayes theorem was used to calculate the posterior probability that a given site is one of those that are selected []. It is worth noting that PAML LRTs are reported to be conservative for short sequences, though the Bayesian prediction of sites under positive selection is largely unaffected by sequence length [].Concerns have been raised over the reliability of particular sites inferred to be subject to positive selection using PAML [] and so corroborating evidence was sought from independent methods. Evidence for positively and negatively selected sites was sought using ADAPTSITE (version 1.3) according to the procedure recommended by Suzuki []. Specifically, equal equilibrium codon frequencies were assumed, an NJ tree based upon p-distance was used as above but was unrooted, the transition/transversion rate ratio was taken to be 1.02 for both primate and mouse datasets (estimated using MEGA) and the significance level for detecting selection was 5%. All three approaches accommodated within ADAPTSITE were employed: maximum parsimony, a distance-based Bayesian method and maximum likelihood. The maximum likelihood analysis was run with two different initial values of ω (0.00001 and 1) to ensure the results were robust to such differences, and only sites where the estimated number of nucleotide substitutions was greater than 15 were considered []. The selection operating in different regions of the sequences and within different branches of the phylogenetic tree under study were also estimated using SWAPSC with 1000 simulated data sets []. This program uses the differences between the estimated and expected numbers of synonymous and nonsynonymous substitutions to evaluate various hypotheses []. Firstly it seeks evidence for regions that have suffered the saturation of synonymous sites: where the number of synonymous substitutions is significantly smaller than expected. In addition it seeks mutational hotspots: regions where the number of synonymous and nonsynonymous nucleotide substitutions are greater than expected under neutrality. Remaining regions where the number of nonsynonymous nucleotide substitutions is smaller than expected (or where ω is significantly smaller than the mean ω estimated for the alignment) are identified as under negative selection. Positive selection is inferred where the estimated number of nonsynonymous nucleotide substitutions is greater than expected by chance and where ω is significantly greater than 1. Where regions have an estimated number of nonsynonymous substitutions greater than expected but ω < 1 or where ω > 1 but there is evidence for saturation of synonymous sites such regions are said to have accelerated rates of nonsynonymous substitutions. Thus SWAPSC seeks to avoid inferring positive selection where there is insufficient data to support it or where saturation may cause bias. A by-product of the SWAPSC analysis is substitution rate estimates for all branches of the tree under study, within overlapping 3 codon windows across the alignment. Each site identified by PAML or ADAPSITE as subject to selection was checked against the synonymous rate estimates made by SWAPSC to ensure that these sites were not saturated at any branch of the tree. [...] In order to establish the relevance of these finding to the solution structure of murine-defensins, we mapped the adaptive sites onto the nuclear magnetic resonance (NMR) structures of mouse Defb7 [] and human DEFB1 []. For each structures were downloaded as PDB files from the Brookhaven protein databank, . The structures were viewed using VMD . […]

Pipeline specifications

Software tools TBLASTN, T-Coffee, EMBOSS, MEGA, PAML, PHYLIP, VMD
Applications Phylogenetics, Protein structure analysis, Amino acid sequence alignment
Organisms Mus musculus, Homo sapiens