\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,romanian,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Mechanism of rapid speciation and the role of environmental factors: A
case study of Aquilegia viridiflora}



\author[1]{Wei Zhang}%
\author[1]{Hua-Ying Wang}%
\author[1]{Tengjiao Zhang}%
\author[1]{Xiaoxue Fang}%
\author[1]{Hongxing Xiao}%
\affil[1]{Northeast Normal University}%


\vspace{-1em}



  
  \date{May 20, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Clarifying the role of key genes in the early stage of speciation is
critical to understand the origin and evolution of species. Aquilegia
viridiflora, in the early stage of speciation, has a wide distribution
range, and shows obvious phenotypic variation among its different
populations. Here, we analyzed the adaptive evolution mechanism of
phenotypic differences in the early stages of speciation using the
phenotypes, genomes, and transcriptomes of different populations of A.
viridiflora. Our results indicated that A. viridiflora originated in
northwestern China, and the environmental changes caused by the uplift
of the northeastern Qinghai--Tibet Plateau in the late Miocene may have
caused its differentiation. Within its distribution range, flower size
was significantly negatively correlated with the inflorescence number
and leaf area. Aqcoe5G459400 expression was reduced in environments with
large temperature differences between day and night, possibly resulting
in longer spurs and promoting the formation of reproductive isolation.
Additionally, the cytochrome P450 (CYP 450) superfamily may have driven
species differentiation in the early stage of speciation. Our study
reveals the genetic basis of the adaptive evolution of the phenotype in
the early stage of speciation and provides new evidence of the rapid
evolution of angiosperms.%
\end{abstract}%



\sloppy


\textbf{INTRODUCTION}

Elucidating the factors contributing to population divergence and
speciation is the core question in evolutionary biology and has
attracted interest from biologists for a long time (Darwin, 1859;
Dobzhansky, 1982; Mayr, 1947; Turelli, Barton, \& Coyne, 2001).
Speciation is a gradual process occurring through the cumulative action
of natural selection and genetic drift in which populations diverge into
groups that cease or limit allele exchange (Coyne \& Orr, 2004; Futuyma,
1998; Raman, 1998). For example, some models of speciation with gene
flow suggest that reproductive isolation (RI) begins at a single locus
but will progressively expand and facilitate divergence of other linked
genes throughout speciation (Rieseberg \& Blackman, 2010; Wu, 2001).

The narrowly differentiated peak in the genome can arise at multiple
prezygotic and postzygotic life-history stages (Coyne \& Orr, 2004;
Haerty \& Singh, 2006; Presgraves, 2010). During population divergence,
prezygotic isolation mechanisms are more important than postzygotic
barriers because they act earlier before hybrid zygotes are formed.
Additionally, morphological divergence is expected to be associated with
prezygotic reproductive barriers. When taxa begin to evolve
morphological differences, they can directly or indirectly reduce
mattings, sperm transfer, or fertilization (Coyne, 2016; Feder et al.,
1994). For example, Turbek et al. found that assortative mating based on
both plumage coloration and song can maintain premating isolation
between \emph{Sporophila iberaensis} and \emph{S. hypoxantha} (Turbek et
al., 2021)\emph{.} This view is also applicable to plants. For many
taxa, including \emph{Aquilegia} , flower color differences leading to
pollinator shifts likely contribute to RI between populations (Des
Marais \& Rausher, 2010; Hodges, Whittall, Fulton, \& Yang, 2002;
Kuriya, Hattori, Nagano, \& Itino, 2015; Quattrocchio et al., 1999;
Schwinn et al., 2006). Additional differences between species that
result in further RI can also accumulate after speciation is already
complete. Furthermore, distinguishing between genetic changes that play
determinant roles in speciation and those that arise after speciation is
complete may be challenging. Thus, organisms in the early stages of
speciation provide an opportunity to identify causal differentiation in
speciation of the genomic landscape.

With the publication of new species, \emph{Aquilegia} currently includes
approximately 110 species worldwide, with approximately 70 species based
on the record of Munz, and is widely distributed in the northern
temperate zone (A. S. Erst et al., 2020; A. S. Erst et al., 2017; Luo,
Erst, Yang, Deng, \& Li, 2018; Munz, 1946). Lu et al. analyzed the
genomes and methylomes of 10 \emph{Aquilegia} plants worldwide and
speculated that the specific genetic and epigenetic structure
endows\emph{Aquilegia} with a high degree of adaptability and promotes
adaptive radiation speciation of the genus (Lu et al., 2019). However,
these species are mature in the late stage of differentiation. In
addition to the gene regions that began to differentiate in the early
stage of speciation, various factors have led to obvious differentiation
of the genome among species. Unlike its congeners, the phenotypes of
different populations of \emph{A. viridiflora} showed high diversity,
which occurs in the early stage of speciation. \emph{A. viridiflora} is
an ideal material to study the key genes that lead to RI in the early
stage of speciation.

\emph{Aquilegia viridiflora} shows rich morphological variation in
different populations, so Andrey S. Erst et al. identified \emph{A.
viridiflora} with laminae purple or lilac-blue petals as the new
species\emph{A. kamelinii} (A. Erst, Shaulo, \& Schmakov, 2013) and
\emph{A. viridiflora} with dark purple petals in North China as the new
species\emph{A. hebeica} (A. S. Erst et al., 2017). However, the
publication of these new species lacks support from molecular data.
According to field surveys, the phenotypic variation of \emph{A.
viridiflora} is continuous (see Figure 1). Therefore, \emph{A.
viridiflora} , \emph{A. kamelinii}and \emph{A. hebeica} were considered
the same species (\emph{A. viridiflora} ) in our study. \emph{A.
viridiflora} is primarily distributed in northern China, extending
between 75.14-132.58\selectlanguage{ngerman}° E and 33.24-53.28° N, in which the east-west span
is wide, and the environment changes gradually. These factors suggest
that environment-mediated divergent selection may be the driving force
for the early differentiation of \emph{A. viridiflora} . Additionally,
other species of\emph{Aquilegia} are produced through ecological
speciation, such as\emph{A. formosa} and \emph{A. pubescens} (Hodges et
al., 2002; Whittall, 2005). Therefore, this species provided a test case
to reveal the effects of adapting to different ecological environments
in the early stage of the speciation process.

In the present study, we collected population samples covering the main
distribution range of \emph{A. viridiflora} for genomic analyses and
morphological characteristic measurements to explore the following: (1)
the morphological divergence during the early stages of speciation, (2)
the genomic basis of such phenotypic traits, and (3) the role of the
environment in maintaining population divergence. Answering these
questions will reveal the speciation process of \emph{A. viridiflora}and
the environmental factors that lead to its differentiation.

\textbf{MATERIALS AND METHODS}

\textbf{Sampling}

Seeds of \emph{Aquilegia viridiflora} were collected from 20 locations
covering its current distribution range, and all voucher specimens were
identified by Dr Hongxing Xiao and deposited in the Northeast Normal
University Herbaria (Figure 1; Table S1). Each of the maternal plants
was separated from the others by \textgreater{} 50 m. Furthermore, to
determine whether the \emph{A. viridiflora} materials collected shared
the most recent common ancestor (MRCA), we also collected \emph{A.
amurensis} , \emph{A. ecalcarata} , \emph{A. japonica} , \emph{A.
oxysepala} var. \emph{kansuensis} , \emph{A. oxysepala}
var.\emph{oxysepala} , and \emph{A. yabeana} seeds, as reported in our
previous study (W. Zhang, Wang, Dong, Zhang, \& Xiao, 2021).
Additionally, seeds of \emph{Paraquilegia microphylla} were collected as
an outgroup. All the seeds were grown in the experimental field of
Northeast Normal University. Fresh leaves and buds were used to extract
DNA and RNA, respectively. These materials were quickly frozen in liquid
nitrogen and stored at -80 \selectlanguage{english}.

\textbf{Phenotypic Measurements and Statistics}

We measured 12 morphological traits of different \emph{A.
viridiflora}populations in the experimental field during the full
flowering stage. To ensure the accuracy of the measurements, we randomly
selected 5 plants from each population to measure the number of
inflorescences and plant height; among these 5 plants, 6 flowers and 6
leaves were randomly selected from each plant, and the corolla diameter,
pistil length, stamen length, leaf area, leaf perimeter and chlorophyll
content were measured and recorded. Among these 6 flowers, we randomly
selected 3 petals and calyxes to measure the petal length (not including
spurs), calyx length, spur length, angle between the petals and spurs
(referred to simply as angle hereafter) and calyx length.

To reduce the error caused by different measurement batches, we used a
mixed linear model to evaluate traits in the lme4 (Bates, Machler,
Bolker, \& Walker, 2014) package in R according to the following
regression model equation:

\begin{equation}
Y_{i}=\mu+\beta_{1}+\beta_{2}P+\varepsilon_{i}\nonumber \\
 \end{equation}

where \(Y_{i}\) represents the traits of different populations,
\(\mu\)is the actual measurement, \(\beta_{1}\) and
\(\beta_{2}\) are regression coefficients, A represents the
measurement batches, P is the person conducting the measurement, and
\(\varepsilon_{i}\) is the residual variance. The evaluation results were
used for ANOVA and K-means cluster analysis in R.

\textbf{DNA Sequencing and SNP Calling}

Three to five individuals were randomly selected from each population
whose phenotype was measured to extract DNA. Total genomic DNA was
extracted using a modified cetyltrimethylammonium bromide (CTAB) method
(Doyle \& Doyle, 1987). Illumina Xten from Biomarker Technologies, Inc.
(Beijing, China) was used for genomic library generation and sequencing
with 2 \selectlanguage{ngerman}× 150 bp paired reads. Furthermore, the raw reads of
\emph{Semiaquilegia adoxoide} s (SRR437677) were downloaded from the
NCBI SRA database (http://www.ncbi.nlm.nih.gov/sra) to be used as an
outgroup. To obtain high-quality genomes, all the reads were assessed
using FastQC (Andrew, 2010) and filtered as follows: reads with adapters
and reads with more than a 10\% N content or more than 50\% low-quality
bases (quality value of less than 10) were removed. Low-quality reads
were removed using NGStoolkit (Mulcare, 2004). The clean reads of each
sample were mapped to the genome sequence of \emph{A. coerulea}
downloaded from the previous study of Filiault et al. (Filiault et al.,
2018) using BWA v.0.7.12 with default parameters (H. Li \& Durbin,
2009). SAMtools v.0.1.18 was used to sort reads (H. Li et al., 2009).
The HaplotypeCaller, GenotypeGVCFs and CombineGVCFs modules in GATK
v.4.1.8.0 were used to produce accurate SNP calls (McKenna et al.,
2010). To improve the quality of SNPs, the VariantFiltration module in
GATK v4.1.8.0 was used for filtration with the following parameters:
--filter-name FilterQual --filter-expression ''QUAL \textless{} 30.0''
--filter-name FilterQD --filter-expression ''QD \textless{} 2.0''
--filter-name FilterMQ --filter-expression ''MQ \textless{} 40.0''
--filter-name FilterFS --filter-expression ''FS \textgreater{} 60.0''
-window 5 -cluster 2. Next, VCFtools v0.1.13 (Danecek et al., 2011) was
used to remove variants that 1) showed a minor allele frequency (MAF) of
0.02 or less, 2) were not balletic variants, 3) showed a sequencing
depth of less than 5, and 4) showed a missing rate exceeding 0.5. The
filtered data were used for subsequent analysis. After filtering, each
SNP was annotated using SnpEff v.5.0 (Cingolani et al., 2012).

\textbf{Population Genetic Structure and Demographic History}

We estimated the phylogenetic relationships of other
\emph{Aquilegia}species with \emph{A. viridiflora} to determine whether
the \emph{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 \emph{A. viridiflora} shared an MRCA. Therefore,
20 populations of \emph{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 \emph{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\textsuperscript{-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).

\textbf{Nucleotide Variation Pattern and Selected Analysis}

Nucleotide diversity (\selectlanguage{greek}π\selectlanguage{english}) was calculated for the four groups using a 100
kb nonoverlapping sliding window and VCFtools (Danecek et al., 2011).
Additionally, the fixation index (Fst) between each of the four groups
was calculated by VCFtools (Danecek et al., 2011) using SNPs. PopLDdecay
(C. Zhang, Dong, Xu, He, \& Yang, 2019) software was applied to compute
linkage disequilibrium (LD) decay among different groups and chromosomes
separately. To avoid the influence of LD on subsequent association and
selection analysis, we employed PLINK v.1.9 (Purcell et al., 2007) to
filter SNPs with the following parameters: --indep-pairwise 50 10 0.2.
The filtered data with 151,577 SNPs were phased using Beagle v.3.3.2
(Browning \& Browning, 2007). To explore the effect of local adaptation
on 12 quantitative traits, we used the single-phenotype Qst-Fst test
with the R package `QstFstComp' (Gilbert \& Whitlock, 2015). If Qst
\textgreater{} Fst, the differentiation of traits is the major effect of
divergent selection and shows local adaptation. We used the half-sib dam
model and 10000 resampling steps for each QstFstComp analysis. To
identify candidate loci under natural selection from 151,577 SNPs, we
used BayeScan v.2.1 (Foll \& Gaggiotti, 2008) software with default
parameters, and PGDSpider (Lischer \& Excoffier, 2012) was used to
produce an input file for BayeScan. SNPs with a q value lower than 0.05
were considered potentially selected loci.

\textbf{Environmental Association Analysis}

Temperature and precipitation are considered the driving factors that
affect the population growth rate and limit the distribution of species
(Cahill et al., 2014; Dalgleish, Koons, Hooten, Moffet, \& Adler, 2011;
Kim \& Donohue, 2013). Therefore, nineteen current bioclimatic variables
with a spatial resolution of 30 s were collected from WorldClim
(http://www.worldclim.org/). At the same time, we recorded the GPS
information of the sampling locations and downloaded the GPS information
of \emph{A. viridiflora}from CVH (http://www.cvh.ac.cn) and GBIF
(http://www.gbif.org). ArcMap v.10.4 was used to limit the spatial
extent according to the buffer radius (5 km) around each occurrence
record. We used \textbar{}r\textbar{} \textless{} 0.8 (Pearson
correlation coefficient) as a cutoff to remove highly correlated
variables. The seven retained current bioclimatic factors (Bio1: annual
mean temperature; Bio2: mean diurnal range; Bio3: isothermality; Bio4:
temperature seasonality; Bio8: mean temperature of the wettest quarter;
Bio15: precipitation seasonality; Bio17: precipitation of the driest
quarter) were used for subsequent analysis. To identify the loci related
to the seven retained current bioclimatic factors, BAYENV2 (Coop,
Witonsky, Di Rienzo, \& Pritchard, 2010) was used with 1,000,000
iterations and run three times separately. For each bioclimatic factor,
the SNPs among the top 1\% according to BF and among the top 5\%
according to the absolute Spearman's \selectlanguage{greek}ρ \selectlanguage{english}were considered candidates.

\textbf{RNA Sequencing and Differential Expression Analysis}

The Qst-Fst test indicated that more floral characteristics showed local
adaptation than other plant characteristics; thus, a bud stage was
selected for transcriptome sequencing. Considering the abovementioned
population genetic analysis results and clustering results for floral
characteristics, ten populations were selected for transcriptome
sequencing. Three samples were collected from each population as
biological replicates to ensure the stability and reliability of the
sequencing results (Table S1). Total RNA was extracted from the samples
using TRIzol (Invitrogen, USA) according to the manufacturer's protocol.
A NanoDrop 2000 spectrophotometer (Thermo Scientific, USA) and gel
electrophoresis were used to estimate the RNA quality. An Illumina Hiseq
system at Novogene Technologies, Inc. (Beijing, China) was used for
genomic library generation and sequencing with 2 \selectlanguage{ngerman}× 150 bp paired reads.
Raw reads were assessed using FastQC (Andrew, 2010) and filtered using
Trimmomatic v.0.39 (Bolger, Lohse, \& Usadel, 2014). The trimmed reads
of each sample were mapped to the same reference genome employed above
using STAR v.2.7.5c (Dobin et al., 2013), and the expression levels were
estimated using RSEM v.1.3.0 (B. Li \& Dewey, 2011). Based on the
obtained read counts, the R package DESeq2 was used to calculate
differential gene expression via pairwise comparisons (Love, Anders, \&
Huber, 2014). Transcripts with an absolute log2FoldChange (LFC) value
greater than 1 were considered significant DEGs. All the DEGs were
analyzed using the clusterProfiler package in R for Kyoto Encyclopedia
of Genes and Genomes (KEGG) enrichment analysis (Yu, Wang, Han, \& He,
2012).

\textbf{Weighted Gene Correlation Network Analysis (WGCNA)}

Groups with the greatest genetic distance and DEGs were selected for
WGCNA. The analysis was conducted using the WGCNA R package following
the provided tutorials (Langfelder \& Horvath, 2008). The soft
thresholding power was determined according to the principle of a
nonscale network, and the lowest power when the correlation coefficient
reached the plateau period was used as a parameter in subsequent
analysis. A gene clustering tree was constructed according to the
correlations between the expression levels of genes and coexpression
modules were identified by dynamic tree cutting using a minimum module
size of 30. When the correlation between the modules and traits was
greater than 0.75 (p \textless{} 0.05), the modules were significantly
related to the traits. The R package clusterProfiler was used for KEGG
enrichment analysis of genes in the phenotype-related modules (Yu et
al., 2012).

\textbf{Identification of Key Genes Related to Speciation and
Adaptation}

To explore the differentiation mechanism in the early stage of
speciation, SNPs related to selection were identified. Candidate SNPs
were mapped to the corresponding genes to check whether these genes were
included in DEGs. The genes obtained using the above filtration
procedure were considered driver genes in the early speciation
of\emph{A. viridiflora} . Next, the clusterProfiler package in R was
used to perform KEGG enrichment analysis on candidate genes (Yu et al.,
2012). Additionally, environment-related genes that were selected and
differentially expressed were regarded as candidate genes that
differentiated by adapting to the environment in the early stages of
speciation. The coding sequence (CDS) regions of candidate genes were
analyzed to identify whether the mutations in the CDS regions were
synonymous or nonsynonymous. Additionally, Multiple Em for Motif
Elicitation (MEME) v.5.3.3 (https://meme-suite.org/meme/tools/meme)
(Bailey \& Elkan, 1994) was used to identify transcription factor
binding sites in the 2 kb region upstream and downstream of the gene,
and Motif Comparison Tool (Tomtom) v.5.3.3
(https://meme-suite.org/meme/tools/tomtom) (Gupta, Stamatoyannopoulos,
Bailey, \& Noble, 2007) was used to annotate the motif functions.

\textbf{RESULTS}

\textbf{Phylogeny and Population Structure}

After filtering, we obtained 1,064,089 high-quality SNPs from 80
individuals, 55.604\% of which were located upstream and downstream,
17.454\% were located in intronic regions, and 7.07\% were located in
the exonic regions of genes. According to the above SNPs, ML and NJ
trees were constructed to infer the phylogenetic relationships
among\emph{Aquilegia} species. Both topologies indicated that different
populations of \emph{A. viridiflora} shared a recently common ancestor
with strong support (Figure S1). Therefore, population genetics analysis
could be performed subsequently on the 20 populations of \emph{A.
viridiflora} , including 672,439 high-quality SNPs. The individuals
of\emph{A. viridiflora} were divided into eastern and western lineages,
each of which contained two groups. For convenience, the four groups
were designated as NE (including the WD, AE, SZ and LF populations), EL
(including the LT, YS, MS and TS populations), CN (including the YT, QB,
TL and YM populations), and NW (including the XW, HL, HH, HD, YS, SF, ZG
and JQ populations) (Figure 2). The NE and EL groups belong to eastern
lineages, and the CN and NW groups belong to western lineages.

PCA was performed to determine the grouping of populations, and
principal component 1 (PC1) and principal component 2 (PC2) explained
6.61\% and 4.04\% of the observed variation, respectively (Figure S2A).
The groups obtained according to the PCA plot were the same as those
indicated in the phylogenetic tree. ADMIXTURE analysis showed that the
optimal \emph{K} value was 4 because \emph{K} = 4 resulted in the lowest
cross-validation error (Figure S2B). Additionally, populations SZ, LT,
YM, XW, HL, HH and HD showed multiple ancestral compositions (Figure 2).

\textbf{Group Dynamic History Inference}

To determine the cause of the multiple ancestral compositions shown in
the ADMIXTURE results, gene flow among the four groups was detected
using the ABBA-BABA test. Tests involving NE and EL, NE and CN, and EL
and CN showed a significant deviation of D-stat from zero (absolute
value of Z-score greater than 3), indicating gene flow between these
groups. Interestingly, no gene flow was detected between the NW and
other groups (Figure 3A and Table S2). Combined with the results of the
ABBA-BABA test, the gene flow pattern and divergence time between
different groups and effective population size (Ne) of each group were
all inferred by G-PhoCS. The gene flow between NE and EL and between NE
and CN was greatest, there was less gene flow from NE to the other two
groups than in the opposite direction, and the gene flow between EL and
CN was smaller. The effective population size was 23,946,300 for the
ancestors. The ancestry population was differentiated into eastern and
western lineages at approximately 6.45 Mya, and their effective
population size were 7,654,150 and 13,578,175, respectively. Next, CN
and NW separated at approximately 6.04 Mya, and NE and EL separated at
5.82 Mya. All the divergence times were in the late Miocene. The current
effective population sizes of NE, EL, CN and NW were 6,566,575,
1,650,775, 1,534,025 and 2,618,550, respectively (Figure 3B and Table
S3). All the Ne values of the other groups were lower than those of
their common ancestor.

The nucleotide diversity (\selectlanguage{greek}π\selectlanguage{english}) of the NE, EL, CN and NW populations was
calculated throughout the genome. Among the four groups, NW showed the
highest nucleotide diversity, and EL showed the lowest nucleotide
polymorphism (Figure S3A). F\textsubscript{ST} was calculated for the
four groups to infer population genetic differentiation. At the overall
level of the genome, the F\textsubscript{ST} values between NE and EL,
NE and CN, NE and NW, EL and CN, EL and NW, and CN and NW were 0.160,
0.151, 0.152, 0.159, 0.193 and 0.161, respectively. These results
indicated that the four groups of \emph{A. viridiflora} had a moderate
level of differentiation, among which EL and NW were the most
differentiated. The F\textsubscript{ST} values calculated for the four
groups using 10 kb windows across genomes were consistent with that at
the overall level of the genome (Figure S3B). Linkage disequilibrium
analysis showed that CN and EL presented a greater degree of LD, while
NE and NW showed less linkage disequilibrium (indicated by
r\textsuperscript{2}). When r\textsuperscript{2} = 0.1, the decay
distances of NE, EL, CN and NW were 10 kb, 37 kb, 44 kb and 6.9 kb,
respectively (Figure S3C). The rapid decay of NE and NW may have been
due to the higher genetic diversity relative to that of EL and CN.

\textbf{Local Phenotypic Adaptation}

ANOVA indicated that each phenotype differed significantly among
different groups (Table 1). K-means cluster analysis of phenotypes
showed that Dim1 and Dim2 could explain 61.3\% and 25.6\% of the
observed variation, respectively, and the cumulative contribution to the
observed phenotypic variation was close to 90\%. All individuals of CN
clustered into one clade; the individuals of the population MS and TS in
the EL group clustered into one clade, which was clearly separated from
the other two populations in EL; the phenotypes of the individuals in NW
and NE showed small differences and could not be clearly distinguished;
and the YS and LT populations in EL were grouped with the individuals of
NW and NE (Figure S4).

Among the four groups, CN showed the smallest corolla diameter, the
shortest petal length, spur length, pistil length, stamen length and
calyx length and the lowest chlorophyll content, followed by EL, while
NW and NE showed the largest values of these phenotypes (the difference
between the last two groups was small). By contrast, plant height and
the number of inflorescences showed the reverse order of the above
traits. EL showed the largest leaf area and perimeter, followed by CN,
while the leaf area of NW and NE presented little difference, but the
leaf perimeter of NW was larger than that of NE. The angle of CN was the
largest, followed by those of NE and NW (the difference between them was
not much), and the angle of EL was the smallest (Figure 4).

The floral characteristics except for the stamen length and angle were
negatively correlated with the number of inflorescences and nutritional
traits except for the chlorophyll content, while the number of
inflorescences was significantly positively correlated with the plant
height (Figure S5). The mean Fst for the four groups was 0.1519
(0.14949-0.15074), and the overall Qst of seven traits was higher than
the mean Fst, including the corolla diameter, petal length, angle, spur
length, pistil length, inflorescence number and leaf area, indicating
that these traits showed local adaptations (Table 1).

\textbf{Identification of DEGs and Gene Modules Analysis}

Because more floral characteristics showed local adaptation, combined
with genomic and phenotypic variation, transcriptome analysis of the
flower buds of ten populations was performed. A total of 1,323,834,776
raw reads (456 Gb) were obtained from 30 samples from ten populations.
The ten populations were also divided into four groups (NE, EL, CN and
NW) for differential expression analysis. The greatest total numbers of
DEGs were found between NE vs. EL, EL vs. NW and CN vs. NW (928, 1716
and 855 DEGs, respectively). Fewer DEGs were found for EL vs. CN, NE vs.
CN and NE vs. NW (428, 666 and 513 DEGs, respectively). Kyoto
Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were
conducted for the DEGs identified between different groups. The DEGs
identified from EL vs. CN and NE vs. CN were not significantly enriched
in any KEGG pathway. The DEGs identified from other groups were related
to glutathione metabolism (ko00480), cytochrome P450 (ko00199),
cyanoamino acid metabolism (ko00460), metabolism of xenobiotics by
cytochrome P450 (ko00980), drug metabolism -- cytochrome P450 (ko00983),
carotenoid biosynthesis (ko00906), starch and sucrose metabolism
(ko00500), amino sugar and nucleotide sugar metabolism (ko00520),
photosynthesis -- antenna proteins (ko00196), galactose metabolism
(ko00052), CD molecules (ko04090), phenylpropanoid biosynthesis
(ko00940), pentose and glucuronate interconversions (ko00040), and
fructose and mannose metabolism (ko00051) (Figure S6 A - D).

Because the greatest number of DEGs and Fst was found between EL and NW,
we used 23,418 genes expressed in both EL and NW for WGCNA to identify
gene modules related to adaptive floral characteristics (including the
corolla diameter, petal length, angle, spur length, pistil length and
inflorescence number). The analysis identified 33 gene modules, the
number of genes in each module ranged from 58 (darkolivegreen) to 3229
(turquoise), and 24 genes were not classified into any modules. Among
the six adaptive floral characteristics, only spur length presented a
significant correlation with four modules --- namely, magenta (640
genes), darkgray (257 genes), saddlebrown (134 genes) and orange (201
genes). No gene modules were associated with other traits (Figure 5).
KEGG enrichment analysis showed that the genes in the magenta module
were related to DNA replication proteins (ko03032), cytoskeletal
proteins (ko04812), the cell cycle (ko04110), protein kinases (ko01001),
DNA replication (ko03030), homologous recombination (ko03440), and gap
junctions (ko04540). The genes in the darkgray module were related to
ribosomes (ko03008) and spliceosomes (ko03040) (Figure S6 E, F). The
genes in the saddle brown and orange modules were not significantly
enriched in any KEGG pathway.

\textbf{Selected Sites at the Speciation and Driven by Environment}

Based on the Bayesian method applied in BAYESCAN, we identified 20,796
outlier SNPs from the 151,577 filtered SNPs according to a 0.05
threshold for the q-value among the four groups. These outlier SNPs
might putatively be considered to have been selected among the four
groups. The outlier SNPs were located in 5954 genes, of which 677 genes
were differentially expressed in different groups. Eleven genes were
enriched in cytochrome P450 (ko00199, CYP450) (Figure S7),
including\emph{Aqcoe2G002300} (CYP98A), \emph{Aqcoe2G223900}
(CYP707A),\emph{Aqcoe4G255000} (CYP71), \emph{Aqcoe4G282900}
(CYP76A),\emph{Aqcoe6G141100} (CYP707A), \emph{Aqcoe4G0100000}
(CYP79D),\emph{Aqcoe5G393400} (CYP71AN), \emph{Aqcoe5G448000}
(CYP90D),\emph{Aqcoe7G425500} (CYP711A), \emph{Aqcoe5G393900} (CYP71),
and\emph{Aqcoe4G241600} (CYP81Q).

Additionally, 1730 outlier SNPs were identified by BAYENV2 that were
associated with seven environmental variables. There were 1,027 and 862
SNPs associated with Bio2 and Bio8, respectively. Considering the
intersection of the SNPs identified by BAYESCAN and BAYENV2, 347 SNPs
located in the intersection might be presumed to be under environmental
selection and were located in 157 genes (Table S4).

\textbf{Key Genes Regulating the Adaptive Evolution of Floral Size}

Among the 157 genes affected by environmental selection, 21 genes were
DEGs. Most of these genes were located on chromosomes 3, 5 and 7, among
which the greatest number of candidate genes was distributed on
chromosome 7 (Figure 6A and Figure S8). The results of QTL mapping also
revealed the existence of genes related to controlling organ morphology
on chromosome 7 of \emph{Aquilegia} (Edwards et al., 2021). Six genes
whose variations were affected by environmental selection showed
nonsynonymous mutations: \emph{Aqcoe1G340100} (Ankyrin repeat family
protein), \emph{Aqcoe1G488100} (ascorbate peroxidase
1),\emph{Aqcoe3G242600} (glutamate dehydrogenase 2),
\emph{Aqcoe5G182600}(cytokinin oxidase 5), \emph{Aqcoe5G459400}
(OBP3-responsive gene 1), and \emph{Aqcoe7G135600} (pleiotropic drug
resistance 11) (Figure 6A). Except for \emph{Aqcoe1G488100} , which was
related to Bio8, the remaining genes were related to Bio2. These six
genes were essential for plant growth and development. Additionally,
\emph{Aqcoe5G459400} and\emph{Aqcoe7G135600} were located in the orange
module and saddle brown module related to spur length, respectively. The
expression of\emph{Aqcoe5G459400} was significantly negatively
correlated with spur length (\emph{p} = 0.007) (Figure 6C), while that
of \emph{Aqcoe7G35600}showed no significant correlation with spur length
(\emph{p} = 0.3481). The selected SNP in \emph{Aqcoe5G459400} showed two
mutations: the NE group was mostly C, the EL group was mostly A, and the
number of A and C in the CN and NW groups was not much different (Figure
6B). Therefore,\emph{Aqcoe5G459400} may be a candidate gene for \emph{A.
viridiflora}spur length to adapt to environmental evolution. WGCNA
showed that the genes that interacting with \emph{Aqcoe5G459400}
were\emph{Aqcoe2G056600} , \emph{Aqcoe2G172400} , \emph{Aqcoe3G439700}
,\emph{Aqcoe7G244300} and \emph{Aqcoe7G395200} . Additionally, a binding
site exists for the ABI3/VP1 transcription factor located 1.8 kb
upstream of the \emph{Aqcoe5G459400} gene.

\textbf{DISCUSSION}

\textbf{Demographic History among \emph{Aquilegia viridiflora}}

The lower eudicot genus \emph{Aquilegia} , which represents a
phylogenetic midpoint between model species such as
\emph{Arabidopsis}and \emph{Oryza} and has recently undergone adaptive
radiation evolution, has become an ideal model to study adaptive
speciation. Although \emph{A. viridiflora} has a wide distribution and
obvious morphological differences, phylogenetic analysis based on the
genome showed that 20 populations constituted a fully supported
monophyletic clade in our study. Both the ML tree and PCA showed four
groups within\emph{A. viridiflora} . Additionally, the optimum K value
according to the ADMIXTURE result was 4 (Figure 2 and Figure S2A), which
verified the existence of four groups. The degree of genome
differentiation was low (Fst values were all less than 0.25), indicating
that \emph{A. viridiflora} is a ``species on the speciation way''.
Additionally, the gene flow from warm and humid regions (EL and CN) to
cold and arid regions (NE) was greater than that in the opposite
direction, which was in accord with studies different populations of
\emph{Silene ciliata} ,\emph{Pseudotsuga menziesii} and \emph{Eperua
falcata} (Brousseau, Fine, Dreyer, Vendramin, \& Scotti, 2021; George et
al., 2021; Morente-L\selectlanguage{ngerman}ópez et al., 2021). This gene flow pattern was
beneficial to increase the genetic diversity of \emph{A. viridiflora} in
cold and arid environments to enhance its local adaptability and
accelerate speciation. The demographic history analysis showed that
\emph{A. viridiflora}differentiated into two lineages, eastern and
western, at 6.45 Mya. Subsequently, the eastern lineage differentiated
into two groups, NE and EL, at 5.97 Mya, and the western lineage
differentiated into CN and NW groups, at 5.54 Mya (Figure 3B). The
divergence times of these groups were all in the late Miocene. At this
time, a global cooling event (late Miocene coding, LMC) caused by the
uplift of the northeastern part of the Qinghai-Tibet Plateau and a
general climate pattern of monsoons in eastern China and arid inland
areas in northwestern China formed (Chen et al., 2019; Steinthorsdottir
et al., 2021). Environmental changes in the late Miocene also drove the
lineage differentiation of the\emph{Psammobates tentorius} complex and
angraecoids in Africa (Farminhão et al., 2021; Zhao, Heideman, Bester,
Jordaan, \& Hofmeyr, 2020). Therefore, temperature may be a critical
driving force for the differentiation of \emph{A. viridiflora} . The
group in northwestern China presented a high level of nucleotide
polymorphism and a low level of linkage disequilibrium, deducing that
northwestern China may be the origin region of \emph{A. viridiflora} .
These groups differentiated in the late Miocene, when one branch
migrated to the North China region and the other branch migrated to the
northeastern region. The group in the northeastern region then diverged
and migrated southward, forming a group in the East Shandong South
Liaoning area.

\textbf{Effect of Cytochrome P450 on the Early Stage of Speciation}

The selected and differentially expressed genes identified in the four
groups of \emph{A. viridiflora} were enriched in cytochrome P450
(ko00199, CYP450) (Figure S7). CYP450 is the largest family of enzyme
proteins in plants and is involved in fatty acid metabolism, antioxidant
biosynthesis, plant defense, secondary metabolite biosynthesis, hormone
regulation and xenobiotic metabolism (Pandian, Sathishraj,
Djanaguiraman, Prasad, \& Jugulam, 2020). These pathways play a key role
in the biosynthesis of secondary metabolites, antioxidants, and hormones
in higher plants. In our study, the 11 genes enriched in CYP450 belong
to eight subfamilies --- CYP71, CYP76, CYP79, CYP81, CYP90, CYP98,
CYP707, and CYP711. For example, CYP98 is located at a key position in
the biosynthesis of phenylpropanoids, and the diversification of
phenylpropanoids provides a driving force for the development and
adaptation of terrestrial plants in new environments (A. Alber \&
Ehlting, 2012; A. V. Alber et al., 2019). ABA 8'-hydroxylase protein is
encoded by the \emph{CYP707} gene and is the key step in ABA catabolism.
It can regulate the ABA level, thereby regulating the physiological
state of plants in response to environmental signals (Janda et al.,
2021; Zheng, Huang, Xian, Wang, \& Liao, 2012). Other subfamilies play a
prominent role in the plant stress response (Irmisch, Zeltner, Handrick,
Gershenzon, \& Köllner, 2015; Johnson et al., 2014; M. Wang et al.,
2020). Additionally, in response to environmental changes, the gene
expression levels of these subfamilies have also changed accordingly
(Baron, Schroeder, \& Stasolla, 2012; Chan, 2012; Tao et al., 2017).
Therefore, the differentiation of the CYP450 superfamily may have driven
the speciation of \emph{A. viridiflora} , and its diversification
provides a genetic basis to adapt to different environments. Because the
CYP450 superfamily is widely present in angiosperms, its important role
in the early speciation of angiosperms warrants further exploration.

\textbf{Phenotype- and gene-driven Speciation under Environment-mediated
Natural Selection}

The differentiation of floral characteristics is considered the core of
prezygotic isolation in plants and plays an important role in speciation
(Celep et al., 2020; Xu et al., 2011). In the present study, the four
identified groups of \emph{A. viridiflora} were located in different
climate types and had continuous phenotypic variation. The areas where
NE and NW were located showed large mean diurnal temperature ranges,
while the areas where EL and CN were located showed small mean diurnal
temperature ranges. Populations located in environments with a wide mean
diurnal temperature range exhibit longer spurs and pistils than those
located in environments with a narrow mean diurnal temperature range.
The differences in the spur and pistil lengths may promote the formation
of reproductive isolation (Brothers \& Delph, 2017; Kramer, 2009; Sun,
Schlüter, Gross, \& Schiestl, 2015). Additionally, plant growth is
greatly affected by the environment; therefore, adaptive adjustments to
different environments are necessary (Anderson \& Song, 2020). Among the
adaptable traits, flower size (corolla diameter, petal length, and spur
length) was significantly negatively correlated with the inflorescence
number and leaf area (Table 1 and Figure S5). NE and NW had larger and
fewer flowers and smaller leaves, while CN and EL had smaller flowers, a
greater number of flowers and larger leaves. This phenomenon has also
been reported in \emph{Pedicularis siphonantha} populations (Dai,
Amboka, Kadiori, Wang, \& Yang, 2017). The wide mean diurnal temperature
range limits the activities of pollinators. To ensure the production of
offspring, \emph{A. viridiflora} uses more energy for reproductive
growth and less energy for vegetative growth. Additionally, the
populations in this area (NW and NE) used the limited energy available
for reproductive growth to form a small number of flowers, thus
enhancing their ability to attract pollinators. By contrast, abundant
pollinators were found in areas (EL and CN) with a narrow mean diurnal
temperature range, and \emph{A. viridiflora} populations in such areas
use more energy for vegetative growth to ensure their long-term
survival. The groups in these areas use their reproductive growth energy
to form more flowers and produce more seeds.

According to Wu's ``Genic view of the process of speciation'', the key
genes or sites of speciation determine the adaptive differentiation of
the population to the environment in which it is located (Wang, He, Shi,
\& Wu, 2020; Wu, 2001). Among selected and differentially expressed
genes related to the environment, \emph{Aqcoe5G182600} (Cytokinin
oxidase 5) might have played an important role in the early speciation
of \emph{A. viridiflora} . Cytokinin oxidase 5 and its orthologous genes
play an important role in plant grain yield.\emph{Os}
\emph{CKX2}knockout results in increases in the numbers of both panicle
branches and grains per plant under stress conditions in \emph{Oryza
sativa}(Ashikari et al., 2005). Accordingly, \emph{Aqcoe5G182600} may
regulate the activity of flower primordia and control the number of
inflorescences by controlling the levels of cytokinin. Furthermore, five
candidate genes were related to resistance to environmental
stress,\emph{Aqcoe1G340100} (Seong et al., 2007),
\emph{Aqcoe1G488100}(Koussevitzky et al., 2008), \emph{Aqcoe3G242600}
(Magadlela et al., 2019), \emph{Aqcoe5G182600} (S. Li et al., 2019)
and\emph{Aqcoe7G135600} (He et al., 2019), indicating that the
differences in the growth environment drive the genetic differentiation
of different groups of \emph{A. viridiflora} .

The appearance of nectar spurs is a critical morphological
characteristic of highly diversified \emph{Aquilegia} . Nosil et al.
proposed that a `speciation gene' is defined as any gene contributing to
the evolution of reproductive isolation (Nosil \& Schluter, 2011).
Because \emph{Aqcoe5G459400}(OBP3-responsive gene 1) was shown to be
related to nectar spurs by integrating the results of WGCNA, and its
expression level was negatively correlated with the spur length, as well
as the selected SNP in \emph{Aqcoe5G459400} was differentiated and
almost fixed in NE and EL, \emph{Aqcoe5G459400} could be considered a
key gene in the early stage of \emph{A. viridiflora} speciation.
Previous studies have shown that OBP3-responsive gene 3 (\emph{ORG3} )
controls the size of petals by controlling the number of cells
(Omidbakhshfard et al., 2018). Additionally, \emph{A. ecalcarata} has
short or no spurs because cell division ceased and cell differentiation
began earlier (Ballerini, Kramer, \& Hodges, 2019). Therefore,
\emph{Aqcoe5G459400} may control the spur size by controlling the number
of cells in \emph{A. viridiflora} . Additionally, the ABI3/VP1
transcription factor located upstream of \emph{Aqcoe5G459400} and its
orthologs are key genes involved in regulating abscisic acid (ABA),
which not only plays a role in the dormancy of seeds and buds but also
affects the flowering time of plants (Riboni, Robustelli Test, Galbiati,
Tonelli, \& Conti, 2016; Shu et al., 2018; Shu et al., 2016; Shu et al.,
2013). In summary, a cascade reaction between ABI3/VP1 and \emph{ORG3}
likely exists, controlling the spur size in \emph{A. viridiflora} . The
expression level of\emph{Aqcoe5G459400} was lower in the lineages from
regions with a wide mean diurnal temperature range, causing \emph{A.
viridiflora} to produce larger flowers to attract more pollinators and
increasing the chance of successful reproduction. By contrast, its
expression level is higher in the lineages from regions with a narrow
mean diurnal temperature range, and \emph{A. viridiflora} produces
smaller flowers in such areas to increase the number of seeds produced.
Thus, the differentiation of\emph{Aqcoe5G49500} enables \emph{A.
viridiflora} to develop different adaptive evolution mechanisms in
response to different environments and promotes speciation. During
speciation of \emph{A. viridiflora} , except for the differentiation of
\emph{Aqcoe5G49500} , the gene flow between populations may increase at
other sites because of the genetic variation mediated by natural
selection, which accelerates the divergence of\emph{A. viridiflora} into
different species (X. Wang et al., 2020). Additionally,
\emph{Aqcoe2G056600} , \emph{Aqcoe2G172400} ,\emph{Aqcoe3G439700} ,
\emph{Aqcoe7G244300} and \emph{Aqcoe7G395200}interacting with
\emph{Aqcoe5G459400} in the module are related to histone methylation
and ion transport across membranes and therefore may improve the ability
of plants to resist insect pests and adverse stresses (Fedoreyeva,
Vanyushin, \& Baranova, 2020; Liao et al., 2018; Philippe, Ralph,
Külheim, Jancsik, \& Bohlmann, 2010; Rodríguez-Celma, Chou, Kobayashi,
Long, \& Balk, 2019; Surya, 2020).

\textbf{AUTHOR CONTRIBUTIONS}

X. HX. and W. HY designed the study and evaluated the results. W. HY and
Z. W. collected the materials. Z. W., Z. TJ. and F. XX. participated in
data analysis. Z. W. and W. HY. prepared the manuscript. All the authors
read and approved the final manuscript.

\textbf{ACKNOWLEDGMENTS}

The research was supported by the National Natural Science Foundation of
China (32070244) and ``the Fundamental Research Funds for the Central
Universities''. We acknowledge Mingzhou Sun for his help in material
collection.

\textbf{DATA AVAILABILITY STATEMENT}

Raw sequence data are available from the National Center for
Biotechnology Information (NCBI) Sequence Read Archive (SRA) under the
submission PRJNA756091 for genome resequencing and PRJNA666554 for
transcriptome sequencing. Phenotypic and climate data have been archived
in dryads (https://doi.org/10.5061/dryad.sqv9s4n4m).

\textbf{ORCID}

\emph{Hongxing Xiao} https://orcid.org/0000-0002-6040-5443

\textbf{REFERENCES}

Alber, A., \& Ehlting, J. (2012). Cytochrome P450s in lignin
biosynthesis. \emph{Advances in Botanical Research, 61} , 113-143.

Alber, A. V., Renault, H., Basilio\selectlanguage{english}-Lopes, A., Bassard, J. E., Liu, Z.,
Ullmann, P., . . . Werck\selectlanguage{english}-Reichhart, D. (2019). Evolution of coumaroyl
conjugate 3\selectlanguage{english}-hydroxylases in land plants: lignin biosynthesis and
defense. \emph{The Plant Journal, 99} (5), 924-936.

Alexander, D. H., Novembre, J., \& Lange, K. (2009). Fast model-based
estimation of ancestry in unrelated individuals. \emph{Genome research,
19} (9), 1655-1664.

Anderson, J. T., \& Song, B. H. (2020). Plant adaptation to climate
change---Where are we? \emph{Journal of Systematics and Evolution, 58}
(5), 533-545.

Andrew, S. (2010). A quality control tool for high throughput sequence
data. \emph{Fast QC, 532} .

Ashikari, M., Sakakibara, H., Lin, S., Yamamoto, T., Takashi, T.,
Nishimura, A., . . . Matsuoka, M. (2005). Cytokinin oxidase regulates
rice grain production. \emph{science, 309} (5735), 741-745.

Bailey, T. L., \& Elkan, C. (1994). Fitting a mixture model by
expectation maximization to discover motifs in bipolymers.

Ballerini, E. S., Kramer, E. M., \& Hodges, S. A. (2019). Comparative
transcriptomics of early petal development across four diverse species
of \emph{Aquilegia} reveal few genes consistently associated with nectar
spur development. \emph{BMC Genomics, 20} (1), 668.

Baron, K. N., Schroeder, D. F., \& Stasolla, C. (2012). Transcriptional
response of abscisic acid (ABA) metabolism and transport to cold and
heat stress applied at the reproductive stage of development in
Arabidopsis thaliana. \emph{Plant Science, 188} , 48-59.

Bates, D., Machler, M., Bolker, B., \& Walker, S. (2014). Fitting linear
mixed-effects models using lme4. \emph{arXiv preprint arXiv:1406.5823} .

Bolger, A. M., Lohse, M., \& Usadel, B. (2014). Trimmomatic: a flexible
trimmer for Illumina sequence data. \emph{Bioinformatics, 30} (15),
2114-2120.

Brothers, A. N., \& Delph, L. F. (2017). Divergence in style length and
pollen size leads to a postmating\selectlanguage{english}-prezygotic reproductive barrier among
populations of Silene latifolia. \emph{Evolution, 71} (6), 1532-1540.

Brousseau, L., Fine, P. V., Dreyer, E., Vendramin, G. G., \& Scotti, I.
(2021). Genomic and phenotypic divergence unveil microgeographic
adaptation in the Amazonian hyperdominant tree \emph{Eperua
falcata}Aubl.(Fabaceae). \emph{Molecular ecology, 30} (5), 1136-1154.

Browning, S. R., \& Browning, B. L. (2007). Rapid and accurate haplotype
phasing and missing-data inference for whole-genome association studies
by use of localized haplotype clustering. \emph{The American journal of
human genetics, 81} (5), 1084-1097.

Cahill, A. E., Aiello\selectlanguage{english}-Lammens, M. E., Caitlin Fisher\selectlanguage{english}-Reid, M., Hua, X.,
Karanewsky, C. J., Ryu, H. Y., . . . Wiens, J. J. (2014). Causes of
warm\selectlanguage{english}-edge range limits: systematic review, proximate factors and
implications for climate change. \emph{Journal of Biogeography, 41} (3),
429-442.

Celep, F., Atalay, Z., Dikmen, F., Do\selectlanguage{romanian}ǧan, M., Sytsma, K. J., \&
Cla\selectlanguage{ngerman}ßen-Bockhoff, R. (2020). Pollination ecology, specialization, and
genetic isolation in sympatric bee-pollinated Salvia
(Lamiaceae).\emph{International Journal of Plant Sciences, 181} (8),
800-811.

Chan, Z. (2012). Expression profiling of ABA pathway transcripts
indicates crosstalk between abiotic and biotic stress responses in
Arabidopsis. \emph{Genomics, 100} (2), 110-115.

Chen, C., Bai, Y., Fang, X., Guo, H., Meng, Q., Zhang, W., . . .
Murodov, A. (2019). A Late Miocene Terrestrial Temperature History for
the Northeastern Tibetan Plateau's Period of Tectonic
Expansion.\emph{Geophysical Research Letters, 46} (14), 8375-8386.

Cingolani, P., Platts, A., Wang, L. L., Coon, M., Nguyen, T., Wang, L.,
. . . Ruden, D. M. (2012). A program for annotating and predicting the
effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome
of Drosophila melanogaster strain w1118; iso-2; iso-3. \emph{Fly, 6}
(2), 80-92.

Coop, G., Witonsky, D., Di Rienzo, A., \& Pritchard, J. K. (2010). Using
environmental correlations to identify loci underlying local
adaptation.\emph{Genetics, 185} (4), 1411-1423.

Coyne, J. A. (2016). Theodosius Dobzhansky on hybrid sterility and
speciation. \emph{Genetics, 202} (1), 5-7.

Coyne, J. A., \& Orr, H. A. (2004). \emph{Speciation} (Vol. 37): Sinauer
Associates Sunderland, MA.

Dai, W.-K., Amboka, G. M., Kadiori, E. L., Wang, Q.-F., \& Yang, C.-F.
(2017). Phenotypic plasticity of floral traits and pollination adaption
in an alpine plant \emph{Pedicularis siphonantha} D. Don when
transplanted from higher to lower elevation in Eastern
Himalaya.\emph{Journal of Mountain Science, 14} (10), 1995-2002.

Dalgleish, H. J., Koons, D. N., Hooten, M. B., Moffet, C. A., \& Adler,
P. B. (2011). Climate influences the demography of three dominant
sagebrush steppe plants. \emph{Ecology, 92} (1), 75-85.

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E.,
DePristo, M. A., . . . Sherry, S. T. (2011). The variant call format and
VCFtools. \emph{Bioinformatics, 27} (15), 2156-2158.

Darwin, C. (1859). The origin of species. Harvard classics. \emph{PF
Collier \& Son, New York} , 1909-1914.

Des Marais, D. L., \& Rausher, M. D. (2010). Parallel evolution at
multiple levels in the origin of hummingbird pollinated flowers in
Ipomoea. \emph{Evolution: International Journal of Organic Evolution,
64} (7), 2044-2054.

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha,
S., . . . Gingeras, T. R. (2013). STAR: ultrafast universal RNA-seq
aligner. \emph{Bioinformatics, 29} (1), 15-21.

Dobzhansky, T. (1982). \emph{Genetics and the Origin of Species} :
Columbia university press.

Doyle, J. J., \& Doyle, J. L. (1987). \emph{A rapid DNA isolation
procedure for small quantities of fresh leaf tissue} . Retrieved from

Edwards, M. B., Choi, G. P., Derieg, N. J., Min, Y., Diana, A. C.,
Hodges, S. A., . . . Ballerini, E. S. (2021). Genetic architecture of
floral traits in bee-and hummingbird-pollinated sister species
of\emph{Aquilegia} (columbine). \emph{bioRxiv} .

Erst, A., Shaulo, D., \& Schmakov, A. (2013). \emph{Aquilegia
kamelinii}(Ranunculaceae)--a new species from North Asia.
\emph{Turczaninowia, 16} (3), 8-10.

Erst, A. S., Pendry, C. A., Erst, T. V., Ikeda, H., Xiang, K., \& Wang,
W. (2020). Two new taxa and one new record of
\emph{Aquilegia}(Ranunculaceae) from India and Pakistan.
\emph{Phytotaxa, 439} (2), 108-118.

Erst, A. S., Wang, W., Yu, S. X., Xiang, K., Wang, J., Shaulo, D. N., .
. . Nobis, M. (2017). Two new species and four new records
of\emph{Aquilegia} (Ranunculaceae) from China. \emph{Phytotaxa, 316}
(2), 121-137.

Farminhão, J. N., Verlynde, S., Kaymak, E., Droissart, V.,
Simo-Droissart, M., Collobert, G., . . . Stévart, T. (2021). Rapid
radiation of angraecoids (Orchidaceae, Angraecinae) in tropical Africa
characterised by multiple karyotypic shifts under major environmental
instability. \emph{Molecular Phylogenetics and Evolution, 159} , 107105.

Feder, J. L., Opp, S. B., Wlazlo, B., Reynolds, K., Go, W., \& Spisak,
S. (1994). Host fidelity is an effective premating barrier between
sympatric races of the apple maggot fly. \emph{Proceedings of the
National Academy of Sciences, 91} (17), 7990-7994.

Fedoreyeva, L. I., Vanyushin, B. F., \& Baranova, E. N. (2020). Peptide
AEDL alters chromatin conformation via histone binding. \emph{AIMS
Biophysics, 7} (1), 1-16.

Filiault, D. L., Ballerini, E. S., Mandakova, T., Aköz, G., Derieg, N.
J., Schmutz, J., . . . Hayes, R. D. (2018). The \emph{Aquilegia} genome
provides insight into adaptive radiation and reveals an extraordinarily
polymorphic chromosome with a unique history. \emph{Elife, 7} , e36426.

Foll, M., \& Gaggiotti, O. (2008). A genome-scan method to identify
selected loci appropriate for both dominant and codominant markers: a
Bayesian perspective. \emph{Genetics, 180} (2), 977-993.

Futuyma, D. J. (1998). Evolutionary biology. In: Sinauer Associates,
Sunderland, MA.

George, J. P., Schueler, S., Grabner, M., Karanitsch\selectlanguage{english}-Ackerl, S., Mayer,
K., Stierschneider, M., . . . van Loo, M. (2021). Looking for the needle
in a downsized haystack: Whole\selectlanguage{english}-exome sequencing unravels genomic signals
of climatic adaptation in Douglas\selectlanguage{english}-fir (\emph{Pseudotsuga menziesii}
).\emph{Ecology and Evolution} .

Gilbert, K. J., \& Whitlock, M. C. (2015). QST--FST comparisons with
unbalanced half\selectlanguage{english}-sib designs. \emph{Molecular ecology resources, 15} (2),
262-267.

Gronau, I., Hubisz, M. J., Gulko, B., Danko, C. G., \& Siepel, A.
(2011). Bayesian inference of ancient human demography from individual
genome sequences. \emph{Nature genetics, 43} (10), 1031.

Gupta, S., Stamatoyannopoulos, J. A., Bailey, T. L., \& Noble, W. S.
(2007). Quantifying similarity between motifs. \emph{Genome Biology, 8}
(2), 1-9.

Haerty, W., \& Singh, R. S. (2006). Gene regulation divergence is a
major contributor to the evolution of Dobzhansky--Muller
incompatibilities between species of Drosophila. \emph{Molecular biology
and evolution, 23} (9), 1707-1714.

He, Y., Xu, J., Wang, X., He, X., Wang, Y., Zhou, J., . . . Meng, X.
(2019). The \emph{Arabidopsis} pleiotropic drug resistance transporters
PEN3 and PDR12 mediate camalexin secretion for resistance to Botrytis
cinerea. \emph{The Plant Cell, 31} (9), 2206-2222.

Hodges, S. A., Whittall, J. B., Fulton, M., \& Yang, J. Y. (2002).
Genetics of floral traits influencing reproductive isolation between
Aquilegia formosa and Aquilegia pubescens. \emph{the american
naturalist, 159} (S3), S51-S60.

Irmisch, S., Zeltner, P., Handrick, V., Gershenzon, J., \& Kollner, T.
G. (2015). The maize cytochrome P450 CYP79A61 produces
phenylacetaldoxime and indole-3-acetaldoxime in heterologous systems and
might contribute to plant defense and auxin formation. \emph{BMC plant
biology, 15} (1), 1-14.

Janda, T., Tajti, J., Hamow, K. A., Mar\selectlanguage{polish}č\selectlanguage{english}ek, T., Ivanovska, B., Szalai,
G., . . . Dark\selectlanguage{ngerman}ó, É. (2021). Acclimation of photosynthetic processes and
metabolic responses to elevated temperatures in
cereals.\emph{Physiologia plantarum, 171} (2), 217-231.

Johnson, S. M., Lim, F.-L., Finkler, A., Fromm, H., Slabas, A. R., \&
Knight, M. R. (2014). Transcriptomic analysis of Sorghum bicolor
responding to combined heat and drought stress. \emph{BMC Genomics, 15}
(1), 1-19.

Kim, E., \& Donohue, K. (2013). Local adaptation and plasticity
of\emph{Erysimum capitatum} to altitude: its implications for responses
to climate change. \emph{Journal of Ecology, 101} (3), 796-805.

Koussevitzky, S., Suzuki, N., Huntington, S., Armijo, L., Sha, W.,
Cortes, D., . . . Mittler, R. (2008). Ascorbate peroxidase 1 plays a key
role in the response of \emph{Arabidopsis thaliana} to stress
combination. \emph{Journal of Biological Chemistry, 283} (49),
34197-34203.

Kramer, E. M. (2009). \emph{Aquilegia} : a new model for plant
development, ecology, and evolution. \emph{Annual review of plant
biology, 60} , 261-277.

Kumar, S., Stecher, G., Li, M., Knyaz, C., \& Tamura, K. (2018). MEGA X:
molecular evolutionary genetics analysis across computing
platforms.\emph{Molecular biology and evolution, 35} (6), 1547-1549.

Kuriya, S., Hattori, M., Nagano, Y., \& Itino, T. (2015). Altitudinal
flower size variation correlates with local pollinator size in a
bumblebee\selectlanguage{english}-pollinated herb, Prunella vulgaris L.(Lamiaceae).\emph{Journal
of evolutionary biology, 28} (10), 1761-1769.

Langfelder, P., \& Horvath, S. (2008). WGCNA: an R package for weighted
correlation network analysis. \emph{BMC bioinformatics, 9} (1), 1-13.

Li, B., \& Dewey, C. N. (2011). RSEM: accurate transcript quantification
from RNA-Seq data with or without a reference genome. \emph{BMC
bioinformatics, 12} (1), 323.

Li, H., \& Durbin, R. (2009). Fast and accurate short read alignment
with Burrows--Wheeler transform. \emph{Bioinformatics, 25} (14),
1754-1760.

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., .
. . Durbin, R. (2009). The sequence alignment/map format and
SAMtools.\emph{Bioinformatics, 25} (16), 2078-2079.

Li, M., Wang, H., Ding, N., Lu, T., Huang, Y., Xiao, H., . . . Li, L.
(2019). Rapid divergence followed by adaptation to contrasting
ecological niches of two closely related columbine
species\emph{Aquilegia japonica} and \emph{A. oxysepala} . \emph{Genome
biology and evolution, 11} (3), 919-930.

Li, S., An, Y., Hailati, S., Zhang, J., Cao, Y., Liu, Y., . . . Yang, P.
(2019). Overexpression of the cytokinin oxidase/dehydrogenase (CKX)
from\emph{Medicago sativa} enhanced salt stress tolerance
of\emph{Arabidopsis} . \emph{Journal of Plant Biology, 62} (5), 374-386.

Liao, Q., Zhou, T., Yao, J. Y., Han, Q. F., Song, H. X., Guan, C. Y., .
. . Singh, A. K. (2018). Genome-scale characterization of the vacuole
nitrate transporter Chloride Channel (CLC) genes and their
transcriptional responses to diverse nutrient stresses in allotetraploid
rapeseed. \emph{PLoS One, 13} (12).

Lischer, H. E., \& Excoffier, L. (2012). PGDSpider: an automated data
conversion tool for connecting population genetics and genomics
programs. \emph{Bioinformatics, 28} (2), 298-299.

Love, M., Anders, S., \& Huber, W. (2014). Differential analysis of
count data--the DESeq2 package. \emph{Genome Biol, 15} (550), 10.1186.

Lu, T., Li, M.-R., Ding, N., Wang, Z.-H., Lan, L.-Z., Gao, X., \& Li,
L.-F. (2019). Genetic and epigenetic mechanisms underpinning the
adaptive radiation of \emph{Aquilegia} species. \emph{bioRxiv} , 782821.

Luo, Y., Erst, A. S., Yang, C.-X., Deng, J.-P., \& Li, L.
(2018).\emph{Aquilegia yangii} (Ranunculaceae), a new species from
western China. \emph{Phytotaxa, 348} (4), 289-296.

Magadlela, A., Morcillo, R. J. L., Kleinert, A., Venter, M., Steenkamp,
E., \& Valentine, A. (2019). Glutamate dehydrogenase is essential in the
acclimation of V\emph{irgilia divaricata} , a legume indigenous to the
nutrient-poor Mediterranean-type ecosystems of the Cape
Fynbos.\emph{Journal of plant physiology, 243} , 153053.

Mayr, E. (1947). Ecological factors in speciation. \emph{Evolution} ,
263-288.

McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K.,
Kernytsky, A., . . . Daly, M. (2010). The Genome Analysis Toolkit: a
MapReduce framework for analyzing next-generation DNA sequencing
data.\emph{Genome research, 20} (9), 1297-1303.

Morente\selectlanguage{english}-Lopez, J., Lara\selectlanguage{english}-Romero, C., Garcia\selectlanguage{english}-Fernandez, A., Rubio Teso, M.
L., Prieto\selectlanguage{english}-Benitez, S., \& Iriondo, J. M. (2021). Gene flow effects on
populations inhabiting marginal areas: origin matters. \emph{Journal of
Ecology, 109} (1), 139-153.

Mulcare, D. M. (2004). NGS Toolkit, Part 8: The National Geodetic Survey
NADCON Tool. \emph{Professional Surveyor Magazine} .

Munz, P. A. (1946). \emph{Aquilegia: the cultivated and the wild
Columbines} : Bailey Hortorium.

Nguyen, L.-T., Schmidt, H. A., Von Haeseler, A., \& Minh, B. Q. (2015).
IQ-TREE: a fast and effective stochastic algorithm for estimating
maximum-likelihood phylogenies. \emph{Molecular biology and evolution,
32} (1), 268-274.

Nosil, P., \& Schluter, D. (2011). The genes underlying the process of
speciation. \emph{Trends in Ecology \& Evolution, 26} (4), 160-167.

Omidbakhshfard, M. A., Fujikura, U., Olas, J. J., Xue, G.-P., Balazadeh,
S., \& Mueller-Roeber, B. (2018). GROWTH-REGULATING FACTOR 9 negatively
regulates arabidopsis leaf growth by controlling ORG3 and restricting
cell proliferation in leaf primordia. \emph{PLoS genetics, 14} (7),
e1007484.

Pandian, B. A., Sathishraj, R., Djanaguiraman, M., Prasad, P., \&
Jugulam, M. (2020). Role of cytochrome P450 enzymes in plant stress
response. \emph{Antioxidants, 9} (5), 454.

Patterson, N., Moorjani, P., Luo, Y., Mallick, S., Rohland, N., Zhan,
Y., . . . Reich, D. (2012). Ancient admixture in human
history.\emph{Genetics, 192} (3), 1065-1093.

Philippe, R. N., Ralph, S. G., Kulheim, C., Jancsik, S. I., \& Bohlmann,
J. (2010). Poplar defense against insects: genome analysis, full-length
cDNA cloning, and transcriptome and protein analysis of the poplar
Kunitz-type protease inhibitor family. \emph{New Phytologist, 184} (4),
865-884.

Presgraves, D. C. (2010). The molecular evolutionary basis of species
formation. \emph{Nature reviews genetics, 11} (3), 175-180.

Price, A. L., Patterson, N. J., Plenge, R. M., Weinblatt, M. E.,
Shadick, N. A., \& Reich, D. (2006). Principal components analysis
corrects for stratification in genome-wide association
studies.\emph{Nature genetics, 38} (8), 904-909.

Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A.,
Bender, D., . . . Daly, M. J. (2007). PLINK: a tool set for whole-genome
association and population-based linkage analyses. \emph{The American
journal of human genetics, 81} (3), 559-575.

Quattrocchio, F., Wing, J., van der Woude, K., Souer, E., de Vetten, N.,
Mol, J., \& Koes, R. (1999). Molecular analysis of the anthocyanin2 gene
of petunia and its role in the evolution of flower color. \emph{The
Plant Cell, 11} (8), 1433-1444.

Raman, S. (1998). \emph{Endless forms: species and speciation} : Oxford
University Press on Demand.

Riboni, M., Robustelli Test, A., Galbiati, M., Tonelli, C., \& Conti, L.
(2016). ABA-dependent control of GIGANTEA signalling enables drought
escape via up-regulation of FLOWERING LOCUS T in \emph{Arabidopsis
thaliana} . \emph{Journal of experimental botany, 67} (22), 6309-6322.

Rieseberg, L. H., \& Blackman, B. K. (2010). Speciation genes in
plants.\emph{Annals of botany, 106} (3), 439-455.

Rodriguez-Celma, J., Chou, H., Kobayashi, T., Long, T. A., \& Balk, J.
(2019). Hemerythrin E3 Ubiquitin Ligases as Negative Regulators of Iron
Homeostasis in Plants. \emph{Other, 10} .

Schwinn, K., Venail, J., Shang, Y., Mackay, S., Alm, V., Butelli, E., .
. . Martin, C. (2006). A small family of MYB-regulatory genes controls
floral pigmentation intensity and patterning in the genus
Antirrhinum.\emph{The Plant Cell, 18} (4), 831-851.

Seong, E.-S., Choi, D.-I., Cho, H.-S., Lim, C.-K., Cho, H.-J., \& Wang,
M.-H. (2007). Characterization of a stress-responsive ankyrin
repeat-containing zinc finger protein of Capsicum annuum
(CaKR1).\emph{BMB Reports, 40} (6), 952-958.

Shu, K., Chen, F., Zhou, W., Luo, X., Dai, Y., Shuai, H., \& Yang, W.
(2018). ABI4 regulates the floral transition independently of ABI5 and
ABI3. \emph{Molecular Biology Reports, 45} (6), 2727-2731.

Shu, K., Chen, Q., Wu, Y., Liu, R., Zhang, H., Wang, P., . . . Liu, C.
(2016). ABI 4 mediates antagonistic effects of abscisic acid and
gibberellins at transcript and protein levels. \emph{The Plant Journal,
85} (3), 348-361.

Shu, K., Zhang, H., Wang, S., Chen, M., Wu, Y., Tang, S., . . . Xie, Q.
(2013). ABI4 regulates primary seed dormancy by regulating the
biogenesis of abscisic acid and gibberellins in Arabidopsis. \emph{PLoS
Genet, 9} (6), e1003577.

Smith, B. J. (2007). boa: an R package for MCMC output convergence
assessment and posterior inference. \emph{Journal of statistical
software, 21} (11), 1-37.

Steinthorsdottir, M., Coxall, H., De Boer, A., Huber, M., Barbolini, N.,
Bradshaw, C., . . . Henderiks, J. (2021). The Miocene: The future of the
past. \emph{Paleoceanography and Paleoclimatology, 36} (4),
e2020PA004037.

Sun, M., Schluter, P. M., Gross, K., \& Schiestl, F. P. (2015). Floral
isolation is the major reproductive barrier between a pair of rewarding
orchid sister species. \emph{Journal of evolutionary biology, 28} (1),
117-129.

Surya, R. (2020). In-silico characterization of \emph{At5g18130} gene
in\emph{Arabidopsis thaliana} with emphasis on its expression patterns
and functional aspects. \emph{Journal of Pharmacognosy and
Phytochemistry, 9} (5), 987-995.

Tao, X., Wang, M.-X., Dai, Y., Wang, Y., Fan, Y.-F., Mao, P., \& Ma,
X.-R. (2017). Identification and expression profile of CYPome in
perennial ryegrass and tall fescue in response to temperature
stress.\emph{Frontiers in plant science, 8} , 1519.

Turbek, S. P., Browne, M., Di Giacomo, A. S., Kopuchian, C., Hochachka,
W. M., Estalles, C., . . . Lovette, I. J. (2021). Rapid speciation via
the evolution of pre-mating isolation in the Ibera
Seedeater.\emph{science, 371} (6536).

Turelli, M., Barton, N. H., \& Coyne, J. A. (2001). Theory and
speciation. \emph{Trends in Ecology \& Evolution, 16} (7), 330-343.

Wang, M., Yuan, J., Qin, L., Shi, W., Xia, G., \& Liu, S. (2020). Ta CYP
81D5, one member in a wheat cytochrome P450 gene cluster, confers
salinity tolerance via reactive oxygen species scavenging. \emph{Plant
biotechnology journal, 18} (3), 791-804.

Wang, X., He, Z., Shi, S., \& Wu, C.-I. (2020). Genes and speciation: is
it time to abandon the biological species concept? \emph{National
Science Review, 7} (8), 1387-1397.

Whittall, J. B. (2005). \emph{Ecological speciation and convergent
evolution in the North American columbine radiation (Aquilegia,
Ranunculaceae)} : University of California, Santa Barbara.

Wu, C. I. (2001). The genic view of the process of
speciation.\emph{Journal of evolutionary biology, 14} (6), 851-865.

Xu, S., Schluter, P. M., Scopece, G., Breitkopf, H., Gross, K.,
Cozzolino, S., \& Schiestl, F. P. (2011). Floral isolation is the main
reproductive barrier among closely related sexually deceptive
orchids.\emph{Evolution: International Journal of Organic Evolution, 65}
(9), 2606-2620.

Yu, G., Wang, L.-G., Han, Y., \& He, Q.-Y. (2012). clusterProfiler: an R
package for comparing biological themes among gene clusters.\emph{Omics:
a journal of integrative biology, 16} (5), 284-287.

Zhang, C., Dong, S.-S., Xu, J.-Y., He, W.-M., \& Yang, T.-L. (2019).
PopLDdecay: a fast and effective tool for linkage disequilibrium decay
analysis based on variant call format files. \emph{Bioinformatics, 35}
(10), 1786-1788.

Zhang, W., Wang, H., Dong, J., Zhang, T., \& Xiao, H. (2021).
Comparative chloroplast genomes and phylogenetic analysis of
Aquilegia.\emph{Applications in plant sciences, 9} (3), e11412.

Zhao, Z., Heideman, N., Bester, P., Jordaan, A., \& Hofmeyr, M. D.
(2020). Climatic and topographic changes since the Miocene influenced
the diversification and biogeography of the tent tortoise (Psammobates
tentorius) species complex in Southern Africa. \emph{BMC evolutionary
biology, 20} (1), 1-33.

Zheng, Y., Huang, Y., Xian, W., Wang, J., \& Liao, H. (2012).
Identification and expression analysis of the Glycine max CYP707A gene
family in response to drought and salt stresses. \emph{Annals of botany,
110} (3), 743-756.

\section*{FIGURE LEGENDS}

{\label{figure-legends}}

Figure 1. Geographical distributions of sampled \emph{Aquilegia
viridiflora} . Populations in the NE, EL, CN and NW groups have
different floral traits. The scale in the figure is 1:20,000,000.

Figure 2. Phylogenetic relationships of \emph{Aquilegia viridiflora}
with genetic structure. Phylogenetic relationships were inferred
according to the maximum likelihood (ML) method, and the genetic
structure showed ADMIXTURE proportions of genetic clusters for each
individual of the four lineages at the best K value (K = 4). Pale
turquoise, moccasin, light green and salmon represent the NE, EL, CN and
NW lineages, respectively.

Figure 3. Demographic history of\emph{Aquilegia viridiflora} . (A)
Results for the D-statistics. The solid red dot represents that the
absolute value of Z-scores is greater than 3, and the hollow red dot
represents that the absolute value of Z-scores is less than 3. (B)
Ancestral population sizes, population divergence times, and migration
rates were assessed by G-PhoCS.

Figure 4. Phenotype of\emph{Aquilegia viridiflora} . (A) Corolla
diameter. (B) Petal length (not including spur). (C) Angle between
petals and spurs. (D) Spur length. (E) Pistil length. (F) Stamen length.
(G) Calyx length. (H) Number of inflorescences. (I) Leaf area. (J)
Chlorophyll content. (K) Leaf perimeter. (L) Plant height.

Figure 5. WGCNA module and trait associations. Gene modules that are
significantly related to a phenotype are indicated in bold.

Figure 6. Key genes for the adaptive evolution of \emph{Aquilegia
viridiflora} in early speciation. (A) Candidate genes for adaptive
evolution that are differentially expressed in EL and NW. The yellow
dots indicate that the SNP in the sliding window was selected by the
environment; red and blue dots represent other regions in the genome.
(B) Gene structure and mutation of \emph{Aqcoe5G459400} . The 13 blue
rectangles represent the CDS of\emph{Aqcoe5G459400} , the red asterisks
represent selected mutations related to the environment in the first
CDS, and the red triangles represent the ABI3/VP1 transcription factor
upstream of\emph{Aqcoe5G459400} . (C) Expression level of
\emph{Aqcoe5G459400} in flowers and the spur length in the different
\emph{A. viridiflora}groups. Individuals representing the spur length of
the NE, EL, CN and NW groups are from the AE, TS, YM and NW populations,
respectively.

\textbf{Hosted file}

\verb`Table.docx`  available at \url{https://authorea.com/users/500580/articles/712244-mechanism-of-rapid-speciation-and-the-role-of-environmental-factors-a-case-study-of-aquilegia-viridiflora}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-2/Figure-2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-3/Figure-3}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-4/Figure-4}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-5/Figure-5}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-6/Figure-6}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

