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).