Computational protocol: Gene flow during glacial habitat shifts facilitates character displacement in a Neotropical flycatcher radiation

Similar protocols

Protocol publication

[…] We used DNA extracts of samples that had been sourced for previous studies by Rheindt et al. [, ]. We sequenced DNA for 46 tissue samples across nine Elaenia species from various museums (Fig. ; Additional file : Table S3), including all species belonging to a monophyletic clade of Elaenia flycatchers of primarily montane distribution, termed the ‘montane Elaenia clade’ by Rheindt et al. [, ], as well as a sample of E. chiriquensis (the sister taxon of the ‘montane Elaenia clade’) as an outgroup. Taxonomic treatment follows the recommendations by Rheindt et al. [, ].Fig. 1 In this study we sequenced nine loci (eight nuclear and one mitochondrial). Five of these (q4, q6, q8, q25, q26) are from among a panel of 26 anonymous loci distributed throughout the genome, primers for which were randomly generated using established protocols [] and trialed for this project (Additional file : Table S4). These five loci yielded clear PCR products and were hence chosen over the remaining anonymous loci. One of our nine loci (02401) was chosen because it amplified best among a panel of 11 loci selected from a list of 242 anonymous markers proposed by Backström et al. [] for population-genetic studies in birds. Another three loci, mitochondrial NADH dehydrogenase subunit 2 (ND2), β-fibrinogen intron 5 (Fib5) and tyrosinase related protein 1 (Tyrp1), are more conventional population-genetic and phylogenetic markers, two of which have amplified Elaenia DNA well in the past (ND2, Fib5; [, ]). Recent genomic data reveal that Tyrp1 is sex linked (present in the Z chromosome) in birds so sequence information of this locus was analyzed accordingly []. General PCR conditions largely followed Rheindt et al. [], with primers and annealing temperatures for the nine loci provided in Additional file : Table S4. We manually edited all sequence chromatograms following Rheindt et al. []. Heterozygous sites were scored as degenerate sites based on established conventions of the International Union of Pure and Applied Chemistry (IUPAC), and were retained as one sequence per individual. All sequences generated in this study are available on GenBank (KY448474 - KY448797). We calculated the number of variable sites and parsimony-informative sites for our dataset in MEGA6 []. [...] We performed genetic analysis for recombinant detection (GARD) [] in order to assess the presence of intragenic recombination for the nuclear loci within our dataset. We conducted separate GARD analyses for each locus in the Datamonkey webserver []. We also performed the Hudson-Kreitman-Aguadé (HKA) test [] as implemented in the program HKA https://bio.cst.temple.edu/~hey/software to assess deviations from neutrality for all loci in our dataset. We used the program SITES https://bio.cst.temple.edu/~hey/software [] to generate input for the HKA test. The HKA test was carried out for the dataset consisting of the Elaenia taxa chilensis, albiceps and pallatangae only as this was our focal group for testing gene flow and isolation. [...] We reconstructed the phylogeny of the ‘montane Elaenia clade’ using both sequence concatenation as well as species tree methods []. We concatenated sequences of all nuclear loci using Sequence Matrix v1.7.8 [] and employed RAxML GUI v1.5 [] to reconstruct a maximum likelihood (ML) tree under the GTR + gamma model, performing 100 runs using a bootstrap algorithm with 1000 replicates. The resulting tree was visualized in FigTree v1.4.2 []. We performed species tree reconstructions using the maximum pseudo-likelihood coalescent method MP-EST [] implemented in the STRAW [] webserver, as well as a Bayesian method implemented in *BEAST v1.8.2 []. The anonymous locus q25 could not be sequenced for E. fallax and hence was removed from species tree analyses.For the MP-EST analysis we first obtained individual gene trees with 1000 bootstraps for the eight loci using RAxML following the same parameters as for the concatenated analysis. We further rooted each gene tree with the E. chiriquensis outgroup using the ReRoot module in the STRAW webserver. The species tree was then estimated in MP-EST. In the Bayesian species tree reconstruction we considered each locus as a separate partition. For each locus, we determined the substitution model through jModelTest v 2.1.7 [, ] (Additional file : Table S5). We used a separate relaxed molecular clock model for each locus and estimated relative clock rates, applying a Yule speciation process with random starting gene trees under two independent runs of 100 million generations each and sampling every 100th generation. We checked parameter convergence of both runs in Tracer v 1.6 [] and combined them in LogCombiner 1.8.2 with a 50% burnin while resampling once per 1000 generations. Subsequently we constructed a maximum clade credibility tree using median heights in TreeAnnotator 1.8.2 while allowing for a further 10% burnin using a posterior probability limit of 0.95. The final tree was visualized in FigTree. In order to explore reticulations and possible species tree topologies, we also used the combined trees obtained from LogCombiner 1.8.2 to generate a representation of topological uncertainty of species tree in DensiTree [] with a 10% burnin.We further explored the mitochondrial phylogeny by determining mitochondrial ND2 haplotypes and constructing a phylogenetic network with the program PopArt 1.7 [] using the TCS method []. [...] We investigated gene flow among the three taxa of the ‘albiceps complex’ within the ‘montane Elaenia clade’ (i.e., chilensis, albiceps, pallatangae) with Bayesian clustering and admixture analyses in BAPS 5.4 [–]. We performed separate analyses on the alignment of the autosomal nuclear loci and the mitochondrial gene ND2. We first performed mixture analyses with ten iterations for a specified genetic cluster (K) and obtained estimates of ancestry coefficients of each individual. Based on our preliminary observations from phylogenetic analyses we performed mixture analyses for both K = 2 and K = 3. We further used the results from mixture analysis to perform population assignment and admixture analysis at both K = 2 and K = 3 using default settings. BAPS requires phased data for nuclear loci, hence we phased our data using PHASE v2.1 [, ] with default settings as implemented in DnaSP 5 []. [...] To further investigate gene flow among the three taxa chilensis, pallatangae and albiceps, we performed coalescent simulations of various evolutionary models in fastsimcoal2 v 2.5.2.2.21 [] and selected the best model using an Approximate Bayesian Computation (ABC) framework in the R package ABC []. Fastsimcoal simulates data following the Wright-Fisher model of evolution assuming neutrality of the genetic markers, no recombination within loci and free recombination between loci and random mating. ABC based approaches can efficiently differentiate between genetic patterns arising from gene flow versus shared ancestral polymorphism [–]. Coalescent simulations were only carried out on loci that showed no significant evidence of intragenic recombination and deviation from neutrality. We simulated a no-gene flow model and various gene flow models in a three-population framework and used alternative topological arrangements resulting from species tree analyses as well as leveraging hypothesized earth-historic events to construct evolutionary models reflecting these scenarios.We also removed the sex-linked locus Tyrp1 for this analysis and from the remaining set of eight loci (seven nuclear autosomal and one mitochondrial) we simulated sequences of the same length as obtained from our sequencing and alignments. We allowed for complete recombination between loci but no intra-locus recombination. We performed 1,000,000 simulations for each model and calculated summary statistics from both the observed data as well as simulated datasets using arlsumstats v. 3.5.2 []. We performed separate simulations for the alignment of nuclear loci and the mitochondrial gene ND2. We set priors for population size (10 to 100,000) based on information obtained from Rheindt et al. []. Following the species tree topology and net interspecies mitochondrial distances observed in this study, and accounting for the uncertainty of divergence estimates between these lineages, we applied broad priors for divergence time parameters (divergence between chilensis and pallatangae: 10,000 to 2,000,000; divergence between albiceps and ancestor of chilensis and pallatangae: 1,000,000 to 7,000,000). We chose a broad prior for migration based on our understanding of species biology (probability of migration: 1e-8 to 1e-4 per generation), with all priors following a uniform distribution. We performed all simulations assuming a constant population size. In our first set of models we assumed various continuous, symmetric and low gene flow scenarios alongside a no-gene flow model. We plotted prior and posterior distributions and verified that the simulations had effectively sampled from the prior distributions for our parameters of interest. We always considered chilensis and pallatangae as sister clades with a much more recent coalescence than the timing of coalescence between albiceps and the ancestral population of chilensis-pallatangae. Whenever necessary, we further refined our models to test more complex scenarios of asymmetric gene flow as well as of isolation and migration patterns.We chose summary statistics that can explain differentiation between populations and their shared variability (Additional file : Table S6). Before performing model selection, we quantified model misclassification using both hard (confusion matrix; see []) and soft classification (mean posterior probability) schemes as implemented in the ABC package []. We estimated the posterior probability of each model from a subset of simulations closest to the observed data (tolerance level) using multinomial logistic regression in the ABC package []. [...] We performed ancestral character state reconstruction of underparts coloration using BayesTraits v 2.0 []. We used a discrete two state coding (yellow versus white) for the terminal species nodes on the tree. The traits for each taxon are: pallatangae (yellow), chilensis (white), albiceps (white), frantzii (white), olivina (yellow), fallax (white), cherriei (white), martinica (white), chiriquensis (white). We used the multistate module in BayesTraits and performed a Bayesian MCMC run with 2,000,000 iterations, including a 100,000 iteration burnin, sampling every 500 generations, and estimated the ancestral state of each node in the final Bayesian species tree. We used an exponential distribution with a hyperprior as suggested by the program manual. We checked for parameter convergence using the coda package [] in R. […]

Pipeline specifications