Population Genetic Structure and Demographic History
We estimated the phylogenetic relationships of other Aquilegiaspecies with A. viridiflora to determine whether the A.
viridiflora materials used in our study shared an MRCA using IQ-TREE
multicore version 1.6.12 (Nguyen, Schmidt, Von Haeseler, & Minh, 2015)
and MEGA X (Kumar, Stecher, Li, Knyaz, & Tamura, 2018) with 1000
bootstrap replicates. Both the ML tree and NJ tree indicated that 20
populations of A. viridiflora shared an MRCA. Therefore, 20
populations of A. viridiflora were used for subsequent population
genetic structure analysis. To infer population genetic structure,
principal component analysis (PCA) was performed using EIGENSOFT v.6.1.4
(Price et al., 2006). ADMIXTURE v.1.3.0 (Alexander, Novembre, & Lange,
2009) was used to investigate the maximum likelihood of the ancestry of
each individual or population. We set the K values from 2 to 10 with 10
replicates for each K value and examined the optimum K value according
to the lowest value of the error rate. By combining the results of the
phylogenetic relationship and genetic structure analysis, we divided 20
populations into four groups.
The Generalized Phylogenetic Coalescent Sampler (G-PhoCS) (Gronau,
Hubisz, Gulko, Danko, & Siepel, 2011) was used to infer the demographic
history of A. viridiflora , including group divergence times,
ancestral population size and migration rates based on neutral loci. The
neutral loci were obtained according to the method of Wang et al.
(2016), and the parameters were automatically set by setting the
‘find-finetunes’ attribute to ‘TRUE’. Because G-PhoCS often shows
limitations in recognizing complex migration models, we performed
inferences under different gene flows between each group based on the
allele frequency by applying the ABBA-BABA statistic using the qpDstat
module in AdmixTools 7.0 (Patterson et al., 2012). Next, we set up a
scenario involving gene flow by combining the results for the ABBA-BABA
statistic. Each Markov chain was run for 2,000,000 generations while
sampling parameter values every 20th iteration. To obtain stable and
reliable results, we independently ran the analysis five times. The
burn-in and convergence of each run were determined using the boa
(Smith, 2007) module in R packages. A neutral mutation rate of
10-8 and a generation time of 1 year were used to
calculate the effective population size, divergence times and migration
rates (M. Li et al., 2019).