Computational protocol: Multilocus Analyses Reveal Postglacial Demographic Shrinkage of Juniperus morrisonicola (Cupressaceae), a Dominant Alpine Species in Taiwan

Similar protocols

Protocol publication

[…] Nucleotide sequences were aligned with the program Clustal X []. Haplotypes of each locus were assigned with the DnaSP v.5.10. Neighbor-joining trees of individual genes were generated using the Kimura 2-parameter model with the program MEGA 5 []. Numbers of mutations between DNA sequence haplotypes based on pairwise comparisons were obtained with mega 5 [], and were used to construct a minimum spanning network with the program minspnet []. [...] Levels of genetic diversity within populations and within species were quantified with measures of nucleotide diversity (π) and haplotype diversity (Hd) using DnaSP v5 []. Patterns of geographical subdivision and levels of genetic differentiation among mountain systems were estimated hierarchically with DnaSP v5. We also calculated Tajima’s D for detecting population growth []. A linear regression was also used to investigate the correlation between average genetic diversity across loci and altitude gradients. The significance of the correlation was determined using a regression F test in the spss program. The samova program was used to identify groups of geographically adjacent populations of J. morrisonicola that were maximally differentiated []. We performed the analyses based on 100 simulated annealing procedures for K values from 2 to 8 and examined maximum indicators of differentiation (FCT values). [...] For microsatellites, genetic diversity within populations/mountain systems was assessed by calculating allele number (A), observed (Ho) and expected heterozygosity (He), and population differentiation (FST) with Arlequin v. 3.5 []. To detect recent genetic bottlenecks, we conducted a Wilcoxon’s sign rank test on the excess of heterozygosity. Bottleneck tests were performed with the bottleneck 1.2.02 using the infinite alleles model (IAM), stepwise mutation model (SMM), and two-phase model (TPM; 95% SMM and 5% IAM) []. We then looked for evidence of genetic structure at a finer resolution by analyzing the data with program structure v. 2.3 []. The program structure v. 2.3 uses a Bayesian method to infer the number of clusters (K) without using prior information of individual sampling locations. The program was run for K = 1 to K = 11 clusters with 10 separate runs for each number of clusters to evaluate the consistency of the results. Each run was pursued for 1 million MCMC interactions, with an initial burn-in of 100,000 and an ancestry model that allowed for admixture. The final posterior probability of K, Pr (X|K), was computed using the runs with highest probability for each K value []. However, as indicated in the structure documentation, Pr (X|K) usually plateaus or increases slightly after the “right” K is reached. Thus, following Evanno et al. [], we calculated ΔK, where the modal value of the distribution is located at the right K. [...] For detecting demographic fluctuations, Bayesian skyline plot (BSP) analyses were used to infer historical demographic dynamics of J. morrisonicola using beast v. 1.7.5 []. For BSP analyses, which use standard MCMC sampling procedures calculated from a sample of molecular sequences for estimating a posterior distribution of effective population size without dependence on a pre-specified parametric model of demographic history, multiple loci can be incorporated to estimate the effective population size through time. In this study, we used beast v. 1.7.5 to determine the demographic history that the species experienced in nucleotide sequence loci and microsatellites, respectively. The substitution model for each gene was set as suggested by JMODELTEST v. 2.1.4 [–]: GTR for trnS-trnG, GTR+I for Maldehy, GTR+I+G for Needly, HKY for trnT-trnL, coxI, and coxIII, and HKY+I for Chs, Myb, and Pgi. The generation time of J. morrisonicola was set as 50 years. Linear and stepwise models were explored using a strict clock with substitution rates of 1.94×10−10 and 1.12×10−10 mutations per site per year for nuclear and organellar loci, respectively [–]. For microsatellite loci, a substitution rate of 2×10−6 was applied to the genetic analyses. Runs consisted of 50 million generations, with trees sampled every 1000 generations. The BSP was visualized in the program Tracer v1.5 [], which summarizes the posterior distribution of population size over time. In this study, we also re-analyzed sequences in species of lower elevations or with a wide altitudinal range in Taiwan, including Castanopsis carlesii (Hemsl.) Hayata (300–700 m a.s.l.) [], Lithocarpus formosana (Hayata) Hayata (ca. 400 m) and L. dodonaeifolius (Hayata) Hayata (600–1,200 m) [], Cyclobalanopsis glauca (Thunb.) Oerst. (ca. 500 m) [], Cycas taitungensis C. F. Shen et al. (400–800 m) [], Trochodendron aralioides Siebold & Zucc. (300–2,500 m) [], Cunninghamia konishii Hayata (1,300–2,000 m) [], Pinus taiwanensis Hayata (750–2,500 m) [], Michelia formosana (Kaneh.) Masam. & Suzuki (100–2,200 m) [], Machilus kusanoi Hayata and M. thunbergii Siebold & Zucc (both ca. 2,300 m) [], Euphrasia transmorrisonensis Hayata (2,200–3,900 m) [], Abies kawakamii (Hayata) Ito (3,000–3,500 m) [], Rhododendron pseudochrysanthum Hayata (3,000–3,900 m) [] ().To further test the demographic hypotheses suggested by the BSP analyses, an approximate Bayesian computation (ABC) was conducted with the program DIYABC v.2.1.0 []. According to the glacial records in Taiwan with the southernmost distribution in Mt. Sanchashan [], individuals of J. morrisonicola were thereby divided into the southern group, with populations of TK and KS, both southern from Shanchashan, and the northern group with all other populations. Four hypothetical scenarios based on the timing of population split and the occurrence/absence of recent expansion were tested: (1) the southern group splitting from the northern group prior to the demographic shrinkage of the northern populations, whereas without subsequent expansions (2) the split occurring after the shrinkage without subsequent expansion, (3) the split occurring after the shrinkage, followed by a recent expansion of the southern populations, and (4) the split occurring prior to the shrinkage, followed by a recent expansion of the southern populations (). Nucleotide sequences and microsatellite data were analyzed separately. For sequence data, nuclear and organellar loci were assigned to different groups; mutation models were both set as Hasegawa-Kishino-Yano, prior distribution of mutation rate was set as uniform with 1×10−9 to 1×10−6 per site per generation, and other parameters were set with the default values. For microsatellites, all loci were assigned to a single group with the default values in all settings. The prior distributions for demographic scenarios and the selected summary statistics were detailed in . Here, we simulated 1 million and 4.5 million datasets for sequences and microsatellites, respectively. The “pre-evaluate scenario-prior combinations” function in DIYABC was used to confirm if all demographic scenarios generated approximate the observed data. The “compute posterior probability of scenarios” function was used to find the best scenario based on the posterior probability estimated by a logistic regression method. […]

Pipeline specifications