Geographic–genomic and geographic–phenotypic differentiation of the Aquilegia viridiflora complex

Abstract

How species diverge into different lineages is a central issue in evolutionary biology. Despite the increasing evidence indicating that such divergences do not need geographic isolation, the correlation between lineage divergence and the adaptive ecological divergence of phenotype corresponding to distribution is still unknown. In addition, gene flow has been widely detected during and through such diverging processes. We used one widely distributed Aquilegia viridiflora complex as a model system to examine genomic differentiation and corresponding phenotypic variations along geographic gradients. Our phenotypic analyses of 20 populations from northwest to northeast China identified two phenotypic groups along the geographic cline. All examined traits are distinct from each other, although a few intermediate individuals occur in their contacting regions. We further sequenced the genomes of representative individuals of each population. However, four distinct genetic lineages were detected based on nuclear genomes. In particular, we recovered numerous genetic hybrids in the contact regions of four lineages. Gene flow is widespread and continuous between four lineages but much higher between contacting lineages than geographically isolated lineages. Gene flow and natural selection might result in inconsistency between heredity and phenotype. Moreover, many genes with fast lineage-specific mutations were identified to be involved in local adaptation. Our results suggest that both geographic isolation and local selection exerted by the environment and pollinators may together create geographic distributions of phenotypic variations as well as the underlying genomic divergences in numerous lineages.

Acknowledgements

We are grateful to Jianquan Liu for suggestions and comments to improve the manuscript; we also thank reviewers for their comments on the manuscript; Mingzhou Sun and Shu Wang for collecting the plant materials; and Peng Peng for providing photographs. The research was supported by the National Natural Science Foundation of China (32070244) and Fundamental Research Funds for the Central Universities.

Author contributions

X.H.X. and W.H.Y. designed the study and evaluated the results; W.H.Y. and Z.W. collected the materials; Z.W., Z.T.J., F.X.X., and L.M.Y participated in data analysis; Z.W. and W.H.Y. prepared the manuscript; all authors read and approved the final manuscript.

Data availability

Raw sequence data are available from the National Center for Biotechnology Information’s (NCBI) Sequence Read Archive (SRA) under the submission PRJNA666554. Data are available from the Dryad Digital Repository (https://doi.org/10.5061/dryad.cnp5hqc8g).

Supplementary data

Supplementary data is available at Horticulture Research online.

References

[1] Wiens JJ, Graham CH . Niche conservatism: integrating evolution, ecology, and conservation biology. Annu Rev Ecol Evol Syst. 2005; 36: 519-39.

[2] Schluter D. Ecology and the origin of species. Trends Ecol Evol. 2001; 16: 372-80.

[3] Sobel JM, Chen GF, Watt LR et al. The biology of speciation. Evolution. 2010; 64: 295-315.

[4] Wu CI . The genic view of the process of speciation. J Evol Biol. 2001; 14: 851-65.

[5] Arias CF, Van Belleghem S, McMillan WO . Genomics at the evolving species boundary. Curr Opin Insect Sci. 2016; 13: 7-15.

[6] Tavares H, Whibley A, Field DL et al. Selection and gene flow shape genomic islands that control floral guides. Proc Natl Acad Sci USA. 2018; 115: 11006-11.

[7] Edelman NB, Frandsen PB, Miyagi M et al. Genomic architecture and introgression shape a butterfly radiation. Science. 2019; 366: 594-9.

[8] Han F, Lamichhaney S, Grant BR et al. Gene flow, ancient polymorphism, and ecological adaptation shape the genomic landscape of divergence among Darwin’s finches. Genome Res. 2017; 27: 1004-15.

[9] Ma T, Wang K, Hu Q et al. Ancient polymorphisms and divergence hitchhiking contribute to genomic islands of divergence within a poplar species complex. Proc Natl Acad Sci USA. 2018; 115: E236-43.

[10] Li LF, Cushman SA, He YX et al. Genome sequencing and population genomics modeling provide insights into the local adaptation of weeping forsythia. Hortic Res. 2020; 7: 130.

[11] Zhang T, Chen J, Zhang J et al. Phenotypic and genomic adaptations to the extremely high elevation in plateau zokor (Myospalax baileyi). Mol Ecol. 2021; 30: 5765-79.

[12] Ayoola AO, Zhang BL, Meisel RP et al. Population genomics reveals incipient speciation, introgression, and adaptation in the African mona monkey (Cercopithecus mona). Mol Biol Evol. 2021; 38: 876-90.

[13] Coyne JA . Theodosius Dobzhansky on hybrid sterility and speciation. Genetics. 2016; 202: 5-7.

[14] Feder JL, Opp SB, Wlazlo B et al. Host fidelity is an effective premating barrier between sympatric races of the apple maggot fly. Proc Natl Acad Sci USA. 1994; 91: 7990-4.

[15] Hodges SA, Whittall JB, Fulton M et al. Genetics of floral traits influencing reproductive isolation between Aquilegia formosa and Aquilegia pubescens. Am Nat. 2002; 159: S51-60.

[16] Kuriya S, Hattori M, Nagano Y et al. Altitudinal flower size variation correlates with local pollinator size in a bumblebeepollinated herb, Prunella vulgaris L. (Lamiaceae). J Evol Biol. 2015; 28: 1761-9.

[17] Quattrocchio F, Wing J, van der Woude K et al. Molecular analysis of the anthocyanin2 gene of petunia and its role in the evolution of flower color. Plant Cell. 1999; 11: 1433-44.

[18] Schwinn K, Venail J, Shang Y et al. A small family of MYBregulatory genes controls floral pigmentation intensity and patterning in the genus Antirrhinum. Plant Cell. 2006; 18: 831-51.

[19] Des Marais DL, Rausher MD . Parallel evolution at multiple levels in the origin of hummingbird pollinated flowers in Ipomoea. Evolution. 2010; 64: 2044-54.

[20] Anderson SAS, Weir JT . The role of divergent ecological adaptation during allopatric speciation in vertebrates. Science. 2022; 378: 1214-8.

[21] Krak K, Vít P, Belyayev A et al. Allopolyploid origin of Chenopodium album s. str. (Chenopodiaceae): a molecular and cytogenetic insight . PLoS One. 2016; 11: e0161063.

[22] Hirota SK, Yasumoto AA, Nitta K et al. Evolutionary history of Hemerocallis in Japan inferred from chloroplast and nuclear phylogenies and levels of interspecific gene flow . Mol Phylogenet Evol. 2021; 164: 107264.

[23] SuarezGonzalez A, Hefer CA, Lexer C et al. Scale and direction of adaptive introgression between black cottonwood (Populus trichocarpa) and balsam poplar (P. balsamifera). Mol Ecol. 2018; 27: 1667-80.

[24] Schluter D. Evidence for ecological speciation and its alternative. Science. 2009; 323: 737-41.

[25] Ficetola GF, Mazel F, Thuiller W . Global determinants of zoogeographical boundaries. Nat Ecol Evol. 2017; 1: 1-7.

[26] ArecesBerazain F, Hinsinger DD, Strijk JS . Genomewide supermatrix analyses of maples (Acer, Sapindaceae) reveal recurring intercontinental migration, mass extinction, and rapid lineage divergence . Genomics. 2021; 113: 681-92.

[27] Song X, Milne RI, Fan X et al. Blow to the northeast? Intraspecific differentiation of Populus davidiana suggests a northeastward skew of a phylogeographic break in East Asia . J Biogeogr. 2021; 48: 187-201.

[28] Wu S, Wang Y, Wang Z et al. Species divergence with gene flow and hybrid speciation on the Qinghai-Tibet plateau. New Phytol. 2022; 234: 392-404.

[29] Wu Z, Raven P, Hong D . Flora of China, Vol. 6 . Beijing: Science Press; 2001.

[30] Erst A, Shaulo D, Schmakov A . Aquilegia kamelinii (Ranunculaceae) - a new species from North Asia . Turczaninowia. 2013; 16: 19-24.

[31] Erst AS, Wang W, Yu SX et al. Two new species and four new records of Aquilegia (Ranunculaceae) from China . Phytotaxa. 2017; 316: 121-37.

[32] Kramer EM . Aquilegia: a new model for plant development, ecology, and evolution . Annu Rev Plant Biol. 2009; 60: 261-77.

[33] Fior S, Li M, Oxelman B et al. Spatiotemporal reconstruction of the Aquilegia rapid radiation through nextgeneration sequencing of rapidly evolving cpDNA regions . New Phytol. 2013; 198: 579-92.

[34] Li M, Wang HY, Ding N et al. Rapid divergence followed by adaptation to contrasting ecological niches of two closely related columbine species Aquilegia japonica and A. oxysepala. Genome Biol Evol. 2019; 11: 919-30.

[35] Zhang W, Wang H, Dong J et al. Comparative chloroplast genomes and phylogenetic analysis of Aquilegia. Appl Plant Sci. 2021; 9: e11412.

[36] Bates D, Mächler M, Bolker B et al. Fitting linear mixedeffects models using lme4. arXiv 2014;1406.5823.

[37] Doyle JJ, Doyle JL . A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem Bull 1987; 19: 11-15.

[38] Andrew S. FastQC. A quality control tool for high throughput sequence data. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/?elqTrackId=3d05e68be90a4fe99b4bdd67ed523fb5&elqaid=2041&elqat=2 (2010).

[39] Patel RK, Jain M . NGS QC Toolkit: a toolkit for quality control of next generation sequencing data. PloS one. 2012; 7: e30619.

[40] Filiault DL, Ballerini ES, Mandáková T et al. The Aquilegia genome provides insight into adaptive radiation and reveals an extraordinarily polymorphic chromosome with a unique history . eLife. 2018; 7: e36426.

[41] Li H, Durbin R . Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009; 25: 1754-60.

[42] Lu T, Li MR, Ding N et al. Genetic and epigenetic mechanisms underpinning the adaptive radiation of Aquilegia species . bioRxiv. 2019; 782821.

[43] Li H, Handsaker B, Wysoker A et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009; 25: 2078-9.

[44] McKenna A, Hanna M, Banks E et al. The genome analysis toolkit: a MapReduce framework for analyzing nextgeneration DNA sequencing data. Genome Res. 2010; 20: 1297-303.

[45] Danecek P, Auton A, Abecasis G et al. The variant call format and VCFtools. Bioinformatics. 2011; 27: 2156-8.

[46] Nguyen LT, Schmidt HA, Von Haeseler A et al. IQTREE: a fast and effective stochastic algorithm for estimating maximumlikelihood phylogenies. Mol Biol Evol. 2015; 32: 268-74.

[47] Alexander DH, Novembre J, Lange K . Fast modelbased estimation of ancestry in unrelated individuals. Genome Res. 2009; 19: 1655-64.

[48] Patterson N, Moorjani P, Luo Y et al. Ancient admixture in human history. Genetics. 2012; 192: 1065-93.

[49] Durand EY, Patterson N, Reich D et al. Testing for ancient admixture between closely related populations. Mol Biol Evol. 2011; 28: 2239-52.

[50] Pickrell J, Pritchard J . Inference of population splits and mixtures from genomewide allele frequency data. Nat Precedings. 2012.

[51] Fitak RR . OptM: estimating the optimal number of migration edges on population trees using Treemix. Biol Methods Protoc. 2021; 6: bpab017.

[52] Huson DH, Bryant D . Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006; 23: 254-67.

[53] Price AL, Patterson NJ, Plenge RM et al. Principal components analysis corrects for stratification in genomewide association studies. Nat Genet. 2006; 38: 904-9.

[54] Blischak PD, Chifman J, Wolfe AD et al. HyDe: a python package for genomescale hybridization detection. Syst Biol. 2018; 67: 821-9.

[55] Leigh JW, Bryant D . POPART: fullfeature software for haplotype network construction. Methods Ecol Evol. 2015; 6: 1110-6.

[56] Excoffier L, Marchi N, Marques DA et al. fastsimcoal2: demographic inference under complex evolutionary scenarios. Bioinformatics. 2021; 37: 4882-5.

[57] Terhorst J, Kamm JA, Song YS . Robust and scalable inference of population history from hundreds of unphased whole genomes. Nat Genet. 2017; 49: 303-9.

[58] Zhang C, Dong SS, Xu JY et al. PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files. Bioinformatics. 2019; 35: 1786-8.

[59] Gilbert KJ, Whitlock MC . QST-FST comparisons with unbalanced halfsib designs. Mol Ecol Resour. 2015; 15: 262-7.

[60] Foll M, Gaggiotti O . A genomescan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics. 2008; 180: 977-93.

[61] Lischer HE, Excoffier L . PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics. 2012; 28: 298-9.

[62] Yu G, Wang LG, Han Y et al. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012; 16: 284-7.

[63] Kim E, Donohue K . Local adaptation and plasticity of Erysimum capitatum to altitude: its implications for responses to climate change . J Ecol. 2013; 101: 796-805.

[64] Dalgleish HJ, Koons DN, Hooten MB et al. Climate influences the demography of three dominant sagebrush steppe plants. Ecology. 2011; 92: 75-85.

[65] Cahill AE, AielloLammens ME, Caitlin FisherReid M et al. Causes of warmedge range limits: systematic review, proximate factors and implications for climate change. J Biogeogr. 2014; 41: 429-42.

[66] Oksanen J, Blanchet FG, Kindt RL et al. vegan: Community Ecology Package , version 2 (2013).

[67] Rousset F . Genetic differentiation and estimation of gene flow from Fstatistics under isolation by distance. Genetics. 1997; 145: 1219-28.

[68] Clarke RT, RotheryP, Raybould AF . Confidence limits for regression relationships between distance matrices: estimating gene flow with distance. Journal of Agricultural, Biological, and Environmental Statistics. 2002; 7: 361-372.

[69] Guillot G, Rousset F . Dismantling the mantel tests. Methods Ecol Evol. 2013; 4: 336-44.

[70] Legendre P, Fortin MJ, Borcard D . Should the Mantel test be used in spatial analysis? Methods Ecol Evol. 2015; 6: 1239-47.

[71] Coop G, Witonsky D, Di Rienzo A et al. Using environmental correlations to identify loci underlying local adaptation. Genetics. 2010; 185: 1411-23.

[72] Galimberti A, Assandri G, Maggioni D et al. Italian odonates in the Pandora’s box: a comprehensive DNA barcoding inventory shows taxonomic warnings at the Holarctic scale. Mol Ecol Resour. 2021; 21: 183-200.

[73] Parisy V, Poinssot B, Owsianowski L et al. Identification of PAD2 as a γglutamylcysteine synthetase highlights the importance of glutathione in disease resistance of Arabidopsis. Plant J. 2007; 49: 159-72.

[74] Mehdi RM, Reza YA, Hassan MM et al. Disease prevention, genetic selection, and vaccination based on BoLADRB3.2 polymorphism: a model for immunogenetic studies. Journal of Advanced Biomedical Sciences. 2022; 12: 1-11.

[75] Eulgem T, Tsuchiya T, Wang XJ et al. EDM2 is required for RPP7dependent disease resistance in Arabidopsis and affects RPP7 transcript levels . Plant J. 2007; 49: 829-39.

[76] Wu L, Chang Y, Wang L et al. Genetic dissection of drought resistance based on root traits at the bud stage in common bean. Theor Appl Genet. 2021; 134: 1047-61.

[77] Manohara, Reddy B, Kumari VN, Johnson A et al. Scarecrow like protein 1,(CtSCL1) involved in drought stress tolerance by interacting with SWI3B component of chromatin modelling complex in cluster bean, Cyamopsis tetragonaloba (L.) Taub. Int J Res Anal Rev. 2018; 5: 304-12.

[78] Rahman T, Shao M, Pahari S et al. Dissecting the roles of cuticular wax in plant resistance to shoot dehydration and lowtemperature stress in Arabidopsis. Int J Mol Sci. 2021; 22: 1554.

[79] Bielsa B, ÁvilaAlonso JI, Fernández i Martí Á et al. Gene expression analysis in cold stress conditions reveals BBX20 and CLO as potential biomarkers for cold tolerance in almond. Horticulturae. 2021; 7: 527.

[80] Beltramino M, Ercoli MF, Debernardi JM et al. Robust increase of leaf size by Arabidopsis thaliana GRF3like transcription factors under different growth conditions . Sci Rep. 2018; 8: 1-13.

[81] Irigoyen S, Karlsson PM, Kuruvilla J et al. The sinkspecific plastidic phosphate transporter PHT4; 2 influences starch accumulation and leaf size in Arabidopsis. Plant Physiol. 2011; 157: 1765-77.

[82] Tyagi S, Mazumdar PA, Mayee P et al. Natural variation in Brassica FT homeologs influences multiple agronomic traits including flowering time, silique shape, oil profile, stomatal morphology and plant height in B. juncea. Plant Sci. 2018; 277: 251-66.

[83] Lim MH, Kim J, Kim YS et al. A new Arabidopsis gene, FLK, encodes an RNA binding protein with K homology motifs and regulates flowering time via FLOWERING LOCUS C . Plant Cell. 2004; 16: 731-40.

[84] Zhu L, Guan Y, Liu Y et al. Regulation of flowering time in chrysanthemum by the R2R3 MYB transcription factor CmMYB2 is associated with changes in gibberellin metabolism. Hortic Res. 2020; 7: 96.

[85] Henderson IR, Liu F, Drea S et al. An allelic series reveals essential roles for FY in plant development in addition to floweringtime control. Development. 2005; 132; 3597-607.

[86] Bollier N, Sicard A, Leblond J et al. AtMINI ZINC FINGER2 and SlINHIBITOR OF MERISTEM ACTIVITY, a conserved missing link in the regulation of floral meristem termination in Arabidopsis and tomato . Plant Cell. 2018; 30: 83-100.

[87] Bartrina I, Otto E, Strnad M et al. Cytokinin regulates the activity of reproductive meristems, flower organ size, ovule formation, and thus seed yield in Arabidopsis thaliana. Plant Cell. 2011; 23: 69-80.

[88] Jiao F, Zhao L, Wu X et al. Metabolome and transcriptome analyses of the molecular mechanisms of flower color mutation in tobacco. BMC Genomics. 2020; 21: 1-10.

[89] Yang X, Bu Y, Niu F et al. Comprehensive analysis of LIM gene family in wheat reveals the involvement of TaLIM2 in pollen development. Plant Sci. 2022; 314: 111101.

[90] Marković V, Cvrčková F, Potocký M et al. EXO70A2 is critical for exocyst complex function in pollen development. Plant Physiol. 2020; 184: 1823-39.

[91] Khosa J, Bellinazzo F, Kamenetsky Goldstein R et al. PHOSPHATIDYLETHANOLAMINEBINDING PROTEINS: the conductors of dual reproduction in plants with vegetative storage organs. J Exp Bot. 2021; 72: 2845-56.

[92] Zhang B, Li C, Li Y et al. Mobile TERMINAL FLOWER1 determines seed size in Arabidopsis. Nature Plants. 2020; 6: 1146-57.

[93] Zhu Y, Klasfeld S, Wagner D . Molecular regulation of plant developmental transitions and plant architecture via PEPB family proteins: an update on mechanism of action. J Exp Bot. 2021; 72: 2301-11.

[94] Ao H, Rohling EJ, Stringer C et al. Twostage midBrunhes climate transition and midPleistocene human diversification. Earth Sci Rev. 2020; 210: 103354.

[95] Rieseberg LH, Soltis D . Phylogenetic consequences of cytoplasmic gene flow in plants. Evol Trends Plants. 1991; 5: 65-84.

[96] Stegemann S, Keuthe M, Greiner S et al. Horizontal transfer of chloroplast genomes between plant species. Proc Natl Acad Sci USA. 2012; 109: 2434-8.

[97] Funk DJ, Omland KE . Specieslevel paraphyly and polyphyly: frequency, causes, and consequences, with insights from animal mitochondrial DNA. Annu Rev Ecol Evol Syst. 2003; 34: 397-423.

[98] Rundell RJ, Price TD . Adaptive radiation, nonadaptive radiation, ecological speciation and nonecological speciation. Trends Ecol Evol. 2009; 24: 394-9.

[99] Jia H, Liu G, Li J et al. Genome resequencing reveals demographic history and genetic architecture of seed salinity tolerance in Populus euphratica. J Exp Bot. 2020; 71: 4308-20.

[100] Hou Z, Li A, Zhang J . Genetic architecture, demographic history, and genomic differentiation of Populus davidiana revealed by wholegenome resequencing . Evol Appl. 2020; 13: 2582-96.

[101] Zhang Q, Zhang X, Yang Y et al. Genetic diversity of Juglans mandshurica populations in Northeast China based on SSR markers . Front Plant Sci. 2022; 13: 931578-8.

[102] Tomasello S, Karbstein K, Hodač L et al. Phylogenomics unravels quaternary vicariance and allopatric speciation patterns in temperatemontane plant species: a case study on the Ranunculus auricomus species complex . Mol Ecol. 2020; 29: 2031-49.

[103] Li JL, Zhong LL, Wang J et al. Genomic insights into speciation history and local adaptation of an alpine aspen in the Qinghai-Tibet plateau and adjacent highlands. J Syst Evol. 2021; 59: 1220-31.

[104] Gao K, Hua Y, Xing LX et al. Speciation of the coldadapted scorpionfly Cerapanorpa brevicornis (Mecoptera: Panorpidae) via interglacial refugia . Insect Cons Divers. 2022; 15: 114-27.

[105] Nei M, Maruyama T, Chakraborty R . The bottleneck effect and genetic variability in populations. Evolution. 1975; 29: 1-10.

[106] Dong T, Yin X, Wang H et al. ABAINDUCED expression 1 is involved in ABAinhibited primary root elongation via modulating ROS homeostasis in Arabidopsis. Plant Sci. 2021; 304: 110821.

[107] Seo M, Aoki H, Koiwai H et al. Comparative studies on the Arabidopsis aldehyde oxidase (AAO) gene family revealed a major role of AAO3 in ABA biosynthesis in seeds . Plant Cell Physiol. 2004; 45: 1694-703.

[108] Xu LL, Yu RM, Lin XR et al. Different rates of pollen and seed gene flow cause branchlength and geographic cytonuclear discordance within Asian butternuts. New Phytol. 2021; 232: 388-403.

[109] Xiao JH, Ding X, Li L et al. Miocene diversification of a goldenthread nanmu tree species (Phoebe zhennan, Lauraceae) around the Sichuan Basin shaped by the east Asian monsoon . Ecology and evolution. 2020; 10: 10543-57.

[110] Hu H, Yang Y, Li A et al. Genomic divergence of Stellera chamaejasme through local selection across the QinghaiTibet plateau and northern China . Mol Ecol. 2022; 31: 4782-96.

[111] Liu S, Li Q, Liu Z . Genomewide identification, characterization and phylogenetic analysis of 50 catfish ATPbinding cassette (ABC) transporter genes. PLoS One. 2013; 8: e63895.

[112] Matsuda S, Funabiki A, Furukawa K et al. Genomewide analysis and expression profiling of halfsize ABC protein subgroup G in rice in response to abiotic stress and phytohormone treatments. Mol Gen Genomics. 2012; 287: 819-35.

[113] Zhang XD, Zhao KX, Yang ZM . Identification of genomic ATP binding cassette (ABC) transporter genes and cdresponsive ABCs in Brassica napus. Gene. 2018; 664: 139-51.

[114] Xu L, Yuan K, Yuan M et al. Regulation of rice tillering by RNAdirected DNA methylation at miniature invertedrepeat transposable elements. Mol Plant. 2020; 13: 851-63.

[115] He Z, Webster S, He SY . Growth-defense tradeoffs in plants. Curr Biol. 2022; 32: R634-9.