\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[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{The Role of \emph{In Situ} Diversification in Island Endemism: Adaptive
Introgression in \emph{Scutellaria} in Taiwan}



\author[1]{Pei-Wei Sun}%
\author[1]{Jui-Tse Chang}%
\author[1]{Min-Xin Luo}%
\author[1]{Chien-Ti Chao}%
\author[2]{Fang Du}%
\author[1]{Pei-Chun Liao}%
\affil[1]{National Taiwan Normal University}%
\affil[2]{Beijing Forestry University}%


\vspace{-1em}



  
  \date{June 18, 2024}


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





\selectlanguage{english}
\begin{abstract}
Island habitats are natural laboratories that offer unique opportunities
to study speciation mechanisms. Recent work indicates that both
anagenetic and cladogenetic speciation processes, driven by allopatric
and sympatric modes, contribute to island species diversity. However,
clear evidence of local adaptation of endemic plant species on islands
requires in-depth studies, which are scarce. This study underscores the
importance of local adaptation in maintaining species entity by
examining how adaptive introgression, hybridization, and local
adaptation contribute to genetic variation in island species.
Specifically, multilocus genome scanning of 51 nuclear genes was used to
investigate the evolutionary relationships of the Scutellaria species
complex on Taiwan Island and assess the role of \emph{in situ}
diversification in generating high endemism and genetic diversity.
Recurrent interspecies introgressions were detected by phylogenetic
networks and ABBA-BABA-based analysis, suggesting ongoing or recent
speciation processes. Approximate Bayesian computation identified hybrid
speciation in S. taiwanensis and S. hsiehii, with evidence of
hybridization between more than two parental species.
Genotype-environment association studies revealed that the influence of
climate, particularly precipitation- and temperature-related factors,
contributed to adaptive genetic divergence between species.
Additionally, adaptive introgressions related to environmental pressures
that may have facilitated the colonization of new island habitats were
identified. This research illustrates how hybridization, introgression,
and adaptation shaped the evolutionary histories and divergence of this
island-endemic plant species complex and sheds light on the multifaceted
mechanisms of speciation on semi-isolated islands.%
\end{abstract}%



\sloppy


\textbf{Introduction}

Introgression is vital for plant evolution (Arnold, 1997; Rieseberg \&
Wendel, 1993; Villa-Mach\selectlanguage{ngerman}ío et al., 2023). The outcome of introgression
depends on the taxa, environmental conditions, and degree of
differentiation (Abbott et al., 2013). Introgression can break species
boundaries and result in genetic homogenization, but it can also
introduce genetic variations that may lead to divergence and promote
speciation along with adaptation.

Introgression between closely related species is most likely to occur
under secondary contact, which has diverse consequences for speciation.
Genetic recombination caused by hybridization can reshape the new
species and isolate it from its parents; for example, in grasshoppers, a
hybrid swarm facilitated the persistence of reproductive isolation from
the parental species (Noguerales \& Ortego, 2022). Interbreeding can
also provide a reservoir of standing genetic variation that produces a
wide range of adaptive outcomes (Barrett \& Schluter, 2008), such as
adaptive introgression across species boundaries, in which beneficial
genetic elements are transferred between species through hybridization
(Suarez-Gonzalez et al., 2018). Detecting the levels, directions, and
drivers of introgression is essential for understanding the speciation
mechanisms of related species on islands. However, incomplete lineage
sorting (ILS), in which alleles from a common ancestor, i.e., ancestral
polymorphisms, are randomly retained, could leave a signature of genetic
admixture similar to that of introgression (Buckley et al., 2006; Du,
Petit, \& Liu, 2009; Funk \& Omland, 2003). Therefore, these two
nonexclusive mechanisms must be distinguished and quantified to unravel
the speciation pattern of a species complex (Funk \& Omland, 2003;
Hedrick, 2013).

In some environments, such as island habitats, local adaptation is
crucial for maintaining species integrity (Algar \& Mahler, 2016; Huang
et al., 2016a). Islands have intrigued taxonomists, evolutionists, and
ecologists since Darwin because, as closed or semi-closed restricted
areas, they offer more opportunities for introgression,
facilitating\emph{in situ} diversification (Ellepola et al., 2022; Ye,
Yang, \& Tian, 2023). When isolated populations encounter diverse and
unique island environments, they face distinct selection pressures that
lead to adaptive divergence. This process involves the accumulation of
genetic variations that enhance the ability of these species to survive
and thrive in their specific habitats (Exposito-Alonso, 2023; Lasky,
Josephs, \& Morris, 2023). For rapidly diverging species with small
population sizes, local adaptation takes on increased importance because
it acts in conjunction with non-adaptive drivers, such as random genetic
drift (Wiens et al., 2022), to strengthen the genetic differentiation
between populations (Rosche et al., 2022). Two processes of
within-island species diversification can be hypothesized: (1) \emph{ex
situ} origination from adjacent areas (e.g., Chiang \& Schaal, 2006) and
(2) \emph{in situ} speciation within the island (e.g., Cowie, 1995). The
former is similar to the anagenetic speciation process involving
colonization of an island, whereas the latter describes cladogenetic
speciation within an island (Stuessy et al., 2006). Emerson and Patino
(2018) suggested the terms ''regional allopatric'' and ''regional
sympatric speciation'' as alternatives to anagenetic and cladogenetic
speciation, respectively. Regardless of whether the speciation mode is
\emph{ex situ} origin or \emph{in situ} , island species live and evolve
within a confined area, which increases contact opportunities and the
possibility of hybridization and introgression. These interbreeding
opportunities accelerate evolution within the island and, in turn,
facilitate \emph{in situ}diversification.

Adaptive introgression may lead to the introduction of novel genetic
combinations that contribute to adaptive traits and the emergence of new
species (Harrison \& Larson, 2014). The transfer of these advantageous
genetic elements across species boundaries can promote the emergence of
individuals with wide adaptability to the island's varying ecological
niches (Grant \& Grant, 2019). The environmental heterogeneity of the
island provides numerous selection pressures that drive the evolution of
these hybrid populations. As a result, the genetic swarm formed through
hybridization enhances the resilience and adaptability of island
species, enabling them to persist and diverge rapidly in isolated
habitats (Grant \& Grant, 2019). In some cases, hybrid speciation may
occur, with the hybrid swarm evolving into a distinct and reproductively
isolated species adapted to the island's specific conditions (e.g.,
Argyranthemum {[}Asteraceae{]} on the Canary Island, Brochmann, Borgen,
\& Stabbetorp, 2000)). Overall, the interplay between local adaptation
and adaptive introgression enriches the genetic diversity and ecological
complexity of island species.

Continental islands, such as Taiwan Island, repeatedly connected and
disconnected from the mainland due to sea level changes during glacial
oscillations (Lin, 1963; Yu, 2003). Continent-island connectivity leads
to the formation of new species from the continental ancestor and
prompts the colonization of islands and founder derivatives from the
continental progenitor, while within-island environmental heterogeneity
facilitates \emph{in situ} diversification (Emerson \& Patino,
2018).\emph{Scutellaria} , commonly known as skullcap, exhibits a high
level of endemism to Taiwan. Among the seven known species, five
(\emph{S. taiwanensis} , \emph{S. hsiehii} , \emph{S. tashiroi} ,
\emph{S. playfairii} , and \emph{S. austrotaiwanensis} ) are endemic,
while the other two (\emph{S. barbata} and \emph{S. indica} ) have wider
distributions across East, South, and Southeast Asia. Among the endemic
species, \emph{S. tashiroi} is distributed in eastern Taiwan and the
Orchid Islet southeast of Taiwan, \emph{S. hsiehii} is confined to
west-central Taiwan, and \emph{S. taiwanensis} , \emph{S.
austrotaiwanensis} , and \emph{S. playfairii} are mainly distributed in
the south. On the basis of three chloroplast DNA fragments and two
low-copy nuclear genes, Chiang, Huang, and Liao (2012) proposed that
endemic \emph{S. tashiroi} , \emph{S. playfairii} , and \emph{S.
austrotaiwanensis} descended from \emph{S. indica} after colonizing
Taiwan and formed the ''\emph{indica} group''. These three
Taiwan-endemic species are thus cladogenetic or regional sympatric
species. By contrast, \emph{S. taiwanensis} and \emph{S. barbata}
(including its synonym \emph{S. taipeiensis} (Chao et al., 2020))
colonized Taiwan independently during the Middle Pleistocene, i.e., 0.4
Mya (Chiang, Huang, \& Liao, 2012), and are so-called anagenetic or
regional allopatric species. The multiple colonization events suggest a
combination of \emph{in situ} speciation and \emph{ex situ} origin of
Taiwan's \emph{Scutellaria} species.

The demographic history of \emph{Scutellaria} on Taiwan Island reveals
frequent southward and northward dispersal and vicariance events
(Chiang, Huang, \& Liao, 2012), and their historical sympatric
distribution has allowed interbreeding (Huang et al., 2017). Their
morphological similarities (Chao et al., 2020; Hsieh, 2013; Hsieh \&
Huang, 1995, 1997) and relatively short speciation history (Chiang,
Huang, \& Liao, 2012) suggest that this group of species has only
recently completed or is still undergoing the speciation process. Their
adaptation to diverse environments and repeated secondary contacts
allowing hybridization and introgression have likely played critical
roles in maintaining their species and genetic diversities. To better
understand how island-endemic species adapt to their surroundings, we
performed multilocus genome scanning using the chromosome-scale genome
of \emph{S. baicalensis} as a reference to examine the evolutionary
relationships among \emph{Scutellaria} species in Taiwan. We analyzed
the occurrence and significance of introgression events between closely
related island species and explored the potential role of these events
in promoting hybrid speciation. Additionally, we investigated the impact
of local adaptation on the maintenance of species boundaries between
rapidly diverging island species.



\textbf{Materials and methods}

\textbf{Sampling}

In 2011--2015, we sampled 27 populations throughout the distribution
range of the seven \emph{Scutellaria} species in Taiwan. Four of these
populations exhibited sympatry: two where both \emph{S. indica}
and\emph{S. barbata} were present, one where both \emph{S.
austrotaiwanensis} and \emph{S. taiwanensis} were present, and one where
both \emph{S. playfairii} and \emph{S. tashiroi} were present (Fig. 1;
Table S1). Species identification of the collected samples was based on
Hsieh (2013), Hsieh and Huang (1995), Hsieh and Huang (1997), and Chao
et al. (2020). Fresh leaves from 4--46 individuals per population were
collected and dried in silica gel. The latitude, longitude, and altitude
of each sampling site were recorded (Table S1). A total of 504
individuals from 27 populations were sequenced in this study, making
this the most comprehensive sampling of \emph{Scutellaria} on Taiwan
Island.

\textbf{Environmental variables}

To characterize the climate conditions at the sampling sites and
identify putative drivers of environment-associated selection, we
downloaded monthly data on 19 bioclimatic variables and average solar
radiation from WorldClim2 (Fick \& Hijmans, 2017). A global aridity
index (Aridity) was added to our analysis of evapotranspiration
processes and rainfall deficits underlying vegetative growth (Trabucco
\& Zomer, 2018). We excluded highly correlated variables with variance
inflation factors (VIF) \textgreater{} 10 and Pearson correlation
coefficients (\emph{r} ) \textgreater{} \textbar{}0.7\textbar{}, leaving
eight climate factors for inclusion in the subsequent analyses: BIO2,
BIO4, BIO5, BIO6, BIO17, BIO18, Solar9, and Aridity.

\textbf{Screening of low-copy-number nuclear genes}

Low-copy-nuclear (LCN) genes are ideal molecular markers for inferring
genetic relationships between plant lineages and identifying potential
mechanisms of selection driven by environmental factors (Du et al.,
2020; Levin, Whelan, \& Miller, 2009; Li et al., 2017). We screened LCN
genes using the genome of \emph{S. baicalensis} and its associated
annotation file as a reference (Zhao et al., 2019). The genes were
screened using two thresholds: (1) gene length \textgreater{} 500 bp and
(2) physical distance \textgreater{} 100 kb. Primer3 (Untergasser et
al., 2012) was used to design primers with a CG content of 40\%--60\%,
length of 18--25 nt, and annealing temperature of 55--65 °C. A total of
132 primer pairs were designed and tested in polymerase chain reaction
(PCR) analysis of 14 individuals representing seven
\emph{Scutellaria}species (Table S2). Amplicons that produced a single
bright band were assumed to be LCN genes and were verified by Sanger
sequencing. Finally, 52 LCN genes that were successfully amplified in
all seven species were retained for subsequent sequencing (Table S6).

\textbf{DNA extraction, library construction, and sequencing}

Genomic DNA was extracted from silica-dried leaves using a modified CTAB
protocol (Porebski, Bailey, \& Baum, 1997). The 52 LCN genes were
amplified by PCR using the protocol of Jiang et al. (2022), and the DNA
library was constructed following de Sousa et al. (2014). After
obtaining fragments of 300--400 bp, PCR was performed on the indexed
samples, and the concentrations of the amplified library were measured
using a NanoDrop 2000c (Thermo Fisher Science, USA). Gene enrichment and
hybridization reactions were performed following the protocol of de
Sousa et al. (2014), and the pooled samples were high-throughput
sequenced using the Illumina MiSeq platform to generate 150-bp
paired-end reads.

\textbf{Read processing, SNP calling, and filtering}

Adapter sequences and low-quality bases were removed from the raw reads
using Trimmomatic (Bolger, Lohse, \& Usadel, 2014). The clean reads were
then mapped to the reference genome of \emph{S. baicalensis} using the
mem algorithm of BWA (Li, 2013). Aligned reads with a mapping quality
score below 20 were discarded. Duplicates were removed with Picard
Toolkit (\emph{http://broadinstitute.github.io/picard/}, accessed in
2019). We used SAMtools (Li et al., 2009) to generate consensus
sequences around target sequences. All genes were aligned in MAFFT
(Katoh \& Standley, 2013), and alignments were trimmed by TrimAL
(Capella-Gutiérrez, Silla-Martínez, \& Gabaldón, 2009). One LCN gene
that was not mapped to any of the nine reference chromosomes was
discarded in the following analyses and left 51 LCN genes.

Finally, single nucleotide polymorphisms (SNPs) were called with
BCFtools (Li, 2011) with default parameters. VCFtools (Danecek et al.,
2011) was used to retain biallelic SNPs with a missing rate \textless{}
0.5, minor allele frequency \textgreater{} 0.01, SNP site with a quality
score \textgreater{} 10, and a mean minimum depth \textgreater{} 3.

\textbf{Genetic diversity and differentiation}

The expected heterozygosity (\emph{H} \textsubscript{exp}), observed
heterozygosity (\emph{H} \textsubscript{obs}),\emph{F}
\textsubscript{IS}, number of polymorphic sites, and \selectlanguage{greek}\emph{π}\selectlanguage{english}for each
population and species were estimated using Arlequin (Excoffier \&
Lischer, 2010). The number of private alleles was calculated by the
package poppr (Kamvar, Tabima, \& Gr\selectlanguage{ngerman}ünwald, 2014) in R 3.6.1(R Core
Team, 2013).

Next, to evaluate the genetic differentiation and ancestry of
the\emph{Scutellaria} species, principal component analysis (PCA) was
performed using adegenet (Jombart, 2008) in R 3.6.1 (R Core Team, 2013).
All SNP sites were retained, and missing data were imputed with mean
values for each SNP site. We further estimated individual assignments
using StrAuto (Chhatre \& Emerson, 2017) with 100,000 MCMC chains,
25,000 steps as burn-in, and ten iterations from \emph{K} = 1 to\emph{K}
= 20. The results were uploaded to CLUMPAK to generate consensus plots
and evaluate the optimal \emph{K} using the \selectlanguage{greek}Δ\selectlanguage{english}\emph{K}method (Kopelman et
al., 2015). Finally, the phylogenetic relationships among the species
network were reconstructed using the NeighborNet program with
\textbackslash{}soutthe concatenated supermatrix of all retained genes
as input implemented with SplitsTree 6 (Huson \& Bryant, 2006).

\textbf{Introgression, incomplete lineage sorting (ILS) and hybrid
speciation}

TreeMix (Pickrell \& Pritchard, 2012) was used to reconstruct the
evolutionary relationships among species using a maximum likelihood tree
and accounting for potential introgression events. We constructed
consensus trees from \emph{m} = 1 to \emph{m} = 10 with 100 bootstrap
replicates; 10 independent runs were executed for each value of\emph{m}
. We estimated the optimal number of migration events using the R
package OptM (Fitak, 2021) by assessing changes in likelihood (\selectlanguage{greek}Δ\selectlanguage{english}\emph{m}
) across incremental values of \emph{m} . The final ML tree with the
best \emph{m} and bootstrap values was visualized according to scripts
written by Zecca, Labra, and Grassi (2020).

To further detect the footprints of introgression, we attempted to
identify hybrid nodes in the phylogenetic networks based on the species
tree (see details in the Supplemental Materials and Methods) using
PhyloNetworks (Sol\selectlanguage{ngerman}ís-Lemus, Bastide, \& Ané, 2017). To reduce the
computational burden, the analyses were run using a subset of 70
individuals randomly selected from each clade and population. Quartet
concordance factors (CF) were calculated from SNP matrices using the R
function SNPs2CF (Olave \& Meyer, 2020). The best number for hybrid
nodes was searched from \emph{hmax} = 1 to \emph{hmax} = 5, and the
scenario with the lowest negative log pseudo-likelihood was evaluated
using bootstrap analysis with the bootsnaq function in PhyloNetworks
with 10 runs per replicate and 100 bootstrap iterations.

We then aimed to detect the specific introgression signal by estimating
Peterson's \emph{D} (ABBA-BABA) and related statistics in Dsuite
(Malinsky, Matschiner, \& Svardal, 2021). \emph{D} and
\emph{f4}statistics were estimated using all possible combinations of
trios within the phylogenetic relationships of the species tree. We used
the\emph{f} -branch test (\emph{fb} ) in Dsuite to allocate gene flows
between specific branches of the species tree, addressing false
positives from \emph{D} -statistics comparisons and ancient
introgression. The \emph{fb} test was also performed to detect the
introgression signal in each gene.

Because of the difficulty of distinguishing recent introgression from
ILS, we used QuIBL (Edelman et al., 2019), a tree-based method that
compares models with ILS and introgression with models with only ILS to
obtain localized information on gene flow. We implemented BIC values to
select the ILS-only model or alternative ILS and introgression model for
each tested triplet. One individual with the highest coverage from each
species was kept, and the sequences of the 52 LCN genes were used as the
input for QuIBL. The percentages of gene sequences and SNPs showing
evidence for introgression were also calculated.

We built and examined two independent sets of evolutionary scenarios
revealed by PhyloNetworks and StrAuto to detect potential introgression
events. First, because the \emph{S. taiwanensis} (TAI) lineage was
genetically admixed with \emph{S. playfairii} (PLA) and \emph{S.
indica}(IND), we investigated the evolutionary relationship of these
three species by focusing on the speciation scenario of TAI (i.e., the
TAI-IND-PLA set)\emph{.} Because the three species are not closely
phylogenetically related and TAI contains genetic components that are
absent from the other two species, we considered the effect of gene flow
from a ghost (unsampled or extinct) lineage to IND or PLA. Regarding
possible\emph{ex situ} sources, the ghost lineage can also be taken as
the hypothetical progenitor of TAI. We proposed three speciation
scenarios (Figure S1 A-C): Recentghost (M1), recent and independent
hybridization between IND and PLA with the ghost lineage; Ancientghost
(M2), ancestral hybridization with the ghost lineage; and Noghost (M3),
ancestral polymorphism. Second, we investigated another genetically
admixed species, \emph{S. hsiehii} (HSI), from \emph{S. indica} lineage
1 (IND1) and \emph{S. austrotaiwanensis} (AUS) (i.e., the HSI-AUS-IND
set). We simulated comparative scenarios of the speciation process
between\emph{S. hsiehii} (HSI), \emph{S. indica} lineage 1 (IND1), and
\emph{S. austrotaiwanensis} (AUS). Three models were designed (Figure S1
D-F): Hybrid (M1), HSI originated from hybridization between IND1 and
AUS; AUSdiv (M2), the speciation of HSI from AUS; and INDdiv (M3), the
speciation of HSI from IND1. A ghost lineage was also added to the
outgroups in all simulated models. All models were simulated and
optimized using the composite likelihoods implemented in fatsimcoal2
(Excoffier et al., 2021), an Approximate Bayesian computation (ABC)
approach. Observed site frequency spectra (SFSs) were generated from
easySFS using the SNP data set. Because the computation of fastsimcoal2
is demanding, observed SFSs were projected to reduce the sample size
while obtaining 90\%--95\% of the maximum segregating sites.

The best model in each evolutionary scenario was selected according to
three criteria. First, the Akaike information criterion (AIC), which
represents the difference between the observed and expected SFSs, was
used to rank the optimized models. Second, the 2D-SFS fit obtained from
the optimized parameters in simulations was visually compared to the
observed SFS. Third, the parametric-bootstrap approach was used to
calculate 100 likelihoods and the 95\% CI for the parameters estimated
by the best models and to select the best-fitting models. The priors and
tested scenarios are given in Supplementary file. The time of divergence
between \emph{S. indica} and the ancestral lineage containing \emph{S.
austrotaiwanensis} and \emph{S. playfairii} was calibrated to 200 kyr
(Chiang, Huang, \& Liao, 2012), with a generation time of one year.

\textbf{Outlier detection and} \textbf{functional annotation of adaptive
genes}

The genetic signal of adaptation was detected using two methods: (1) an
approach that identifies SNPs with significantly high or low genetic
differentiation compared to the overall population structure (\emph{F}
\textsubscript{ST} outliers) implemented in BayeScan (Foll, 2012) and
(2) a landscape genomic method that uncovers the molecular mechanism of
adaptation based on genotype-environment association (GEA) analysis
applied in BayEnv2 (Feng \& Du, 2022; Günther \& Coop, 2013). In brief,
BayeScan was performed with a parameter of posterior probability = 100
and detected loci that showed higher genetic differentiation than the
q-value threshold \textless{} 0.05 (Foll, 2012). BayEnv2 was used to
estimate the covariance of allele frequencies between populations, and
the model was then used as a null hypothesis for detecting SNP sites
associated with environmental variables (Günther \& Coop, 2013).

To detect loci strongly associated with the most important axes for
rearranging genetic components according to environmental factors, the
RDA-based method was used (Capblancq et al., 2018). To minimize false
positives, demography and geography were accounted for in the RDA-based
method by implementing two separate RDA-based strategies as described in
Capblancq et al. (2023): the first kept climate variables as predictors
without conditioning other factors (RDA-raw), whereas the second kept
climate variables as predictors while using the first two genetic PC
axes from the genetic PCA and the first three distance-based Moran's
eigenvector maps (dbMEMs) as conditioning variables (RDA-correct). Loci
shared by at least two of the four scanning methods were retained as
loci putatively under selection by the environment.

A hierarchical approach was used for functional annotation of the
outlier loci. First, KOBAS-I (Bu et al., 2021) was used to annotate the
trimmed nucleotide sequences of all retained genes using the reference
genome of \emph{Arabidopsis thaliana} with a threshold of e-value
\textless{} 10\textsuperscript{-5}. Second, the complete amino acid
sequences of the orthologous genes retrieved from the genome annotation
of \emph{S. baicalensis} for all unmapped genes in the first round were
retrieved and functionally annotated using KOBAS-I in the same manner as
in the first round.

\textbf{Drivers of adaptive genetic variation}

To evaluate the relative contributions of demographic history,
environment, and geography to shaping the genetic variation
of\emph{Scutellaria} in Taiwan, partial redundancy analysis (pRDA) was
performed. First, the geographic coordinates were converted into dbMEMs,
and the first three axes were retained as proxies of geographic
structure. Second, PCA was performed using the allele frequencies among
populations, and PC1 and PC2 were retained as proxies of demographic
history. Variation partitioning was carried out using the allele
frequencies of each population as responses and either of the groups of
predictors, including (1) geographic structure, (2) demographic history,
and (3) the eight environmental factors (Fig. 3C) we selected by VIF as
explanatory variables, while conditioning the remaining predictors. The
variation explained by the full model and each testable portion was
assessed and tested for significance. Significant predictors were
selected by forward selection. The pRDA analyses, dbMEM, and forward
selection were conducted using the R package vegan (Oksanen et al.,
2007).

\textbf{Results}

\textbf{Genetic diversity and differentiation}

After quality filtering, 224,000 reads (SD = 36,000) were recovered per
individual, and a mean coverage of 127 x (SD = 20 x) was obtained for
the target regions. Genes that did not map to any of the nine
chromosomes of the reference genome were discarded in the following
analyses. A total of 489 biallelic SNPs (average 9.4 (SD = 7.4) SNP
sites per gene) were recovered from the enriched genes after filtering
with VCFtools (Table S6).

The genetic diversity varied between populations and species (Table S3).
The species-specific \emph{H} \textsubscript{exp} ranged from 0.18 to
0.3 and was highest in \emph{S. taiwanensis} (0.3) and \emph{S.
hsiehii}(0.26). Moreover, private alleles were most abundant in the
non-endemic species \emph{S. barbata,} with 61 private alleles, followed
by another non-endemic species, \emph{S. indica} , with 16 private
alleles.

In \emph{S. barbata} and \emph{S. indica} , the intraspecies\emph{F}
\textsubscript{ST} values were larger than the\emph{F}
\textsubscript{ST} values between species (Fig. S3A).\emph{F}
\textsubscript{IS} was negative in all populations (Table S3). The PCA
indicated considerable intraspecific differentiation in \emph{S.
barbata} , \emph{S. indica} , and \emph{S. taiwanensis} . \emph{S.
tashiroi} and \emph{S. playfairii} formed distinct clusters within
species in the first three axes (Fig. 1B). Strong genetic structure was
observed in \emph{S. barbata} , \emph{S. indica} , \emph{S. taiwanensis}
.\emph{S. tashiroi} , and \emph{S. playfairii,} and the two best K
values revealed substantial admixture between most species (Fig. 1A, C;
Fig. S2). \emph{S. taiwanensis} was grouped with some populations of
\emph{S. indica} and \emph{S. tashiroi} . A similar phenomenon was
observed in\emph{S. hsiehii} , which contained genetic components from
\emph{S. austrotaiwanensis} and some \emph{S. indica} populations (Fig.
1 C). The FH and CC populations of \emph{S. barbata} contained genetic
components from \emph{S. tashiroi} and an unknown source. The AWD and
INDWL populations of \emph{S. indica} included genetic components
from\emph{S. tashiroi} and \emph{S. hsiehii} . The HDD population of
\emph{S. indica} and the PLAWLLD population of \emph{S. playfairii} also
contained genetic components from \emph{S. tashiroi} . Conversely, the
TASWLLD population of \emph{S. tashiroi} included genetic components
from \emph{S. playfairii} . On the basis of the results of the STRUCTURE
analysis, we conclude that hybridization and introgression occurred
frequently between Taiwanese \emph{Scutellaria} species.

Consistent with these findings, the distance-based phylogenetic network
revealed a web-like and reticulated evolutionary history; not all
species were well resolved as monophyletic groups (Fig. 2). \emph{S.
barbata, S. indica, S. tashiroi,} and \emph{S. taiwanensis} each formed
two distinct phylogenetic groups. Within \emph{S. tashiroi} , TAS1 was
resolved as a sister to PLA, whereas TAS2 was resolved as a sister to
IND1 (Fig. 2). \emph{S. indica} exhibited similar relationships with
other species, with IND1 and IND2 resolved as sisters to TAI2 and TAS1,
respectively.

\textbf{Contributions of ILS and introgression to genetic variation}

QuIBL analysis indicated evidence of introgression at 32\% (34/105) of
trios and 46\% (227/489) of SNPs (Table S4). In addition,
Patterson's\emph{D} and \emph{f4} inferred multiple significant
interspecies introgression events. The most frequently introgressed
species were\emph{S. barbata} and \emph{S. taiwanensis} , which
introgressed with 5 and 4 species, respectively (Fig. S4). By contrast,
\emph{fb} statistics revealed only one introgression event, which
occurred between \emph{S. barbata} and the ancestor of the group
containing \emph{S. indica} ,~\emph{S. austrotaiwanensis} ,~\emph{S.
tashiroi} , \emph{S. hsiehii,} and \emph{S. playfairii} (Fig. 3B). The
\emph{fb} statistics showed significant signals of gene flow in
\emph{LYC} , \emph{OPCL1} , and\emph{PHYB} (Fig. 3A; Table 1). These
genes were associated with two {[}BIO2 (mean diurnal range) and BIO18
(precipitation of warmest quarter){]}, six {[}BIO2, BIO4 (temperature
seasonality), BIO5 (maximum temperature of warmest month), BIO6 (minimum
temperature of coldest month), BIO17 (temperature annual range) and
BIO18{]}, and one (BIO2) environmental variables related to
precipitation and temperature, respectively (Table 1). A recent
admixture between \emph{S. tashiroi}and \emph{S. hsiehii} was inferred
by TreeMix (Fig. 3B; Fig. S5). Finally, the results of PhyloNetworks
analysis indicated that the data were best fit with the lowest BIC by
the model containing one reticulation event, and all bootstrap
replications supported the introgression between the ancestor of
\emph{S. playfairii} and \emph{S. taiwanensis} (Fig. 3B).

\textbf{Evolutionary history of genetic admixture in \emph{S.
taiwanensis} and \emph{S. hsiehii}}

The results of fastsimcoal2 analysis indicated that the Recentghost and
Hybrid models best explained the evolution of \emph{S. taiwanensis}
and\emph{S. hsiehii} , respectively (Fig. 4). For the first set,
TAI-IND-PLA, the AIC and bootstrapping likelihood comparison of all
models revealed that recent introgression from \emph{S. playfairii}
and\emph{S. indica} to \emph{S. taiwanensis} was the optimal model (Fig.
S7; Fig. 4A). This model also had the highest congruence with the
observed SFS (Figs. S8-S10). The best-fitting model suggested that a
ghost lineage diverged from the common ancestor at approximately 700
kya. This divergence was followed by more recent independent
introgression events with \emph{S. indica} and \emph{S. playfairii} at
7.62 kya and 4.77 kya, respectively, with \textasciitilde{}60\% of the
genetic component inherited from introgression (Fig. 4A; Fig. S11; Table
S7). The time of divergence from the ghost lineage (700 kya) aligns with
the coalescent time estimated in the Stairway plot (Fig. S12;
Supplemental Materials and Methods), whereas the introgression time
corresponds to the onset time of population decline (Fig. S12). Our
approach supports that TAI originated from multiple introgressions of
more than one \emph{Scutellaria} species.

For the second set, HSI-AUS-IND, we again constructed and compared three
speciation scenarios to investigate the evolutionary history of \emph{S.
hsiehii} . Although the AIC values did not informatively distinguish
between hybrid speciation (Hybrid model) and divergence from \emph{S.
austrotaiwanensis} (AUSdiv model) (Fig. 4B)\emph{,} the bootstrapping
likelihoods strongly supported the model corresponding to hybridization
between \emph{S. austrotaiwanensis} and \emph{S. indica} (Fig. S13-16;
Table S8). However, the genetic components provided by each species in
the estimated hybridization model differed (Fig. S4C; Table S8): more
than 90\% of the genetic components in \emph{S. hsiehii} were inherited
from \emph{S. austrotaiwanensis} approximately 9.07 kya according to the
results of STRUCTURE, indicating another substantial but potential
hybrid speciation event in the studied species.

\textbf{Contributions of climate to genetic divergence}

Eight climate factors were included in the analysis of the contribution
of climate to genetic divergence (Fig. 3C). A total of 272 SNPs were
scanned by BayeScan and three GEA methods as significant outliers (Fig.
3D; Fig. S6). We observed a substantial intersection between the results
of RDA and BayEnv regardless of whether the former accounted for
population structure; 31 and 41 of the SNPs detected by BayEnv (i.e.,
24\% of each set) were detected by RDA.raw and RDA.correct,
respectively. We also observed a sizable intersection (i.e., 33 SNPs;
20\%) between the RDA.raw and RDA.correct sets. The 89 SNPs that were
detected by at least two methods were considered putative adaptive SNPs
(Fig. S3D). Forward selection also identified genetic structure
(PC1-PC2) and BIO18 as the only significant predictors explaining the
genetic allele frequencies of adaptive loci (Table S5). The 89 adaptive
SNPs belonged to 34 of the 51 genes enriched in this study. Many of
these are protein-coding genes that are known to be involved in the
adaptation of other species. A complete list of genes, potential
functions, and associated climate variables is given in Table S6.

Climate, demography, and geography together explained 71\% of the
variance by pRDA (Table 2; Fig. S6). When two of the three factors were
controlled, demography was the only significant predictor, explaining
52.1\% of the total variance, whereas pure climate and pure geography
explained 1\% and 2\% of the total variance, respectively. These results
suggest that the breaking of species boundaries by genetic groupings
(e.g., Fig. 2B) due to hybridization may play a stronger role in
enhancing the genetic diversity of these species within Taiwan Island
than differentiation caused by geographic or climate barriers. Notably,
16.3\% of the variance was confounded among the predictors (Table 2),
suggesting some degree of collinearity among genetic structure, climate,
and geography.

\textbf{Discussion}

The present examination of 51 nuclear genes shows that frequent
hybridization and introgression events were important for \emph{S.
taiwanensis} origination and \emph{in situ} speciation of
the\emph{indica} group through increased genetic diversity. Standing
variation is crucial during the early stages of plant speciation and
governs the repeatability of speciation processes (Bock et al., 2023).
Incorporating standing variation through hybridization helps hybrid
organisms survive and reach new adaptive peaks. Hybrids that acquire new
adaptive genes and subsequent mutations that function later in the
speciation process (Bock et al., 2023) can widen their ecological
potential, adapt to new habitats, and reduce vulnerability to climate
change (Brauer et al., 2023).

Chiang, Huang, and Liao (2012) established that \emph{S. barbata} is the
most divergent species (and a heterotypic synonym of \emph{S.
taipeiensis} (Chao et al., 2020; Hsiung et al., 2017)) in the phylogeny,
followed by \emph{S. taiwanensis} and the ''\emph{indica} group'' (which
includes \emph{S. indica} , \emph{S. tashiroi} , \emph{S. playfairii}
,\emph{S. austrotaiwanensis} (Chiang, Huang, \& Liao, 2012), and the new
species \emph{S. hsiehii} (Hsieh, 2013)). Hybridization has occurred not
only within the \emph{indica} group but also between the
\emph{indica}group and \emph{S. barbata} and \emph{S. taiwanensis} ,
implying a lack of reproductive isolation between species. We focused on
the two endemic species with genetic admixture, \emph{S. taiwanensis}
and \emph{S. hsiehii,} which both have hybrid origins. Hybridization of
the progenitor of \emph{S. taiwanensis} or an extinct species (i.e.,
ghost) with other \emph{in situ} diverged species brought \emph{S.
taiwanensis}closer to the \emph{indica} group in the phylogeny. The
speciation time for sympatric speciation (i.e., the \emph{indica} group)
was relatively recent, leading to weak reproductive barriers with
frequent introgressions. These frequent hybridizations and
introgressions widened the ecological potential of \emph{Scutellaria}
species to adapt to the highly heterogeneous habitats of Taiwan Island
by increasing genomic complexity, consistent with the hybrid model,
which showed that hybrid genome instability enhances adaptation to wild
environments (Porretta \& Canestrelli, 2023; Stelkens \& Bendixsen,
2022). In this study, we used the \emph{Scutellaria} species limited to
Taiwan Island as a model system to explore the relationships of frequent
hybridization with origin and adaptation.

\textbf{Genetic evidence of adaptative variation}

To identify adaptive genetic variants, we used multivariate approaches,
which perform better than univariate procedures while maintaining a low
false-positive rate (Forester et al., 2018; Sang et al., 2022). However,
we observed significant overlap between the climate and demographic
factors, indicating high multicollinearity (Capblancq et al., 2023). To
balance potentially missing weak-effect sites and minimize
false-positive rates, we adopted a strategy in which sites detected by
more than two approaches were considered outliers (Ahrens et al., 2018;
Forester et al., 2018). We deliberately avoided relying solely on sites
detected by all methods simultaneously. This approach was chosen to
increase the true-positive rate while circumventing limitations inherent
in the least potent methods.

Chiang, Huang, and Liao (2012) proposed repeated dispersal and
vicariance within the current geographic distribution but did not
address how local adaptation contributes to species divergence. The
detection of adaptive SNPs points to local adaptation within these
species. Notably, partial RDA revealed 15.2\% overlap in variation
between climate factors and genetic structure, indicating that climate
adaptation has played a crucial role in maintaining species divergence.
Among the genes implicated in climate-related adaptation, \emph{CBF2}and
\emph{PHYB} function in the response and acclimation to cold
in\emph{Arabidopsis thaliana} (Jiang et al., 2020; Zhao et al.,
2016).\emph{ABA1} and \emph{HSA32} play essential roles in responding to
heat and temperature stimuli in plants (Charng et al., 2006; Suzuki et
al., 2016), and \emph{YUC6} and \emph{CRY2} have been linked to drought
recovery in \emph{Arabidopsis thaliana} (Cha et al., 2015; Mao et al.,
2005). \emph{PKT3} and \emph{MYB4} are associated with UV-B and
high-light-intensity responses in plants (He et al., 2016; Zhao et al.,
2007). The adaptive divergence of these climate-related genes suggests
that local climate has been a significant driver of the rapid speciation
and maintenance of species divergence in Taiwanese
\emph{Scutellaria}species.

Moreover, previous studies have highlighted the importance of positive
selection in the diversification of Taiwanese \emph{Scutellaria} . For
example, Huang et al. (2015) identified interspecific differences in the
expression of growth and environmental regulation genes,
specifically\emph{R2R3-MYB} , suggesting that these genes are associated
with positive selection and functional gene divergence in
Taiwanese\emph{Scutellaria} species. Positive selection of the
anthocyanin biosynthesis gene \emph{UFGT} has also been linked to
accelerated floral diversity on Taiwan Island (Huang et al., 2016a).
Moreover, adaptive subfunctionalization of the anthocyanin gene
\emph{DFR} was instrumental in anthocyanin diversification (Huang et
al., 2016b), which is ecologically relevant in terms of attracting
pollinators, repelling herbivores and parasites, and providing abiotic
stress resistance (e.g., UV resistance) (Lev-Yadun \& Gould, 2009).

In summary, the genetic variations of Taiwanese
\emph{Scutellaria}species have been shaped by a combination of factors.
Hybridization and ghost introgression provide greater genetic variation
by increasing standing genetic variation, and further local adaptation
driven by climate-related genes and adaptive introgression has likely
contributed to the ability of these species to thrive in specific
environments. Furthermore, the positive selection of functional genes
has played a crucial role in driving their diversification and
ecological characteristics. These findings shed light on the intricate
mechanisms that shaped the evolutionary history of \emph{Scutellaria} on
Taiwan Island in response to ecological selection pressures.

\emph{S. indica} \textbf{and \emph{S. playfairii} as sources of
hybridization of \emph{S. taiwanensis}}

The fastsimcoal2 simulation indicated that \emph{S. taiwanensis} arose
from hybridizations of a ghost lineage with \emph{S. indica} and\emph{S.
playfairii} approximately 7.62 and 4.77 kya\emph{,}respectively,
suggesting multiple origins of extant \emph{S. taiwanensis} . Older
specimens of \emph{S. taiwanensis} are often mistaken morphologically
for \emph{S. indica} , and the habitat of\emph{S. taiwanensis} on shady
slopes is similar to that of \emph{S. playfairii} (Hsieh \& Huang,
1995). Despite this multiple-source hybridization, \emph{S. taiwanensis}
is considered a single species since the two hybrid lineages (TAI1 and
TAI2) are monophyletic in multispecies coalescent (ASTRAL) inferences.
Hybridization events may displace certain gene trees, leading to gene
tree incongruence or topological inconsistencies between multispecies
coalescent and concatenated analyses (Degnan \& Rosenberg, 2009), as
supported by the low \emph{F} \textsubscript{ST} with other
\emph{indica} group members (Fig. S3A). Approximately 60\% of the
genetic components of \emph{S. taiwanensis} are from \emph{indica} group
members.

\emph{S. taiwanensis} is scattered in the shallow mountains of southern
Taiwan and partially overlaps with \emph{S. playfairii} , enabling
hybridization between the ancestor of \emph{S. taiwanensis} (ghost
lineage) and \emph{S. playfairii} . Surprisingly, \emph{S.
taiwanensis}hybridized with \emph{S. indica} , a species that has no
distributional overlap with \emph{S. taiwanensis} and is distributed in
north-central Taiwan. In addition to long-distance dispersal, we cannot
exclude the possibility that an ancestor or close relative of \emph{S.
indica}spread to the south (Chiang, Huang, \& Liao, 2012) and hybridized
with\emph{S. taiwanensis} after the glacial climate began. Although the
middle and low altitudes of Taiwan Island did not experience glacial
coverage, the arid and cold glacial climate forced many organisms to
hide in multiple relict refugia during the Pleistocene glacials. In
Taiwan, glacial refugia are commonly found in mountainous areas in the
south or southeast (Huang, Hwang, \& Lin, 2002; Lee et al., 2006; Wu et
al., 2006) and Mount Sylvia (Xue-shan) in the northeast (Cheng, Hwang,
\& Lin, 2005; Huang, Hwang, \& Lin, 2002; Lee et al., 2006; Wu et al.,
2006). The warm and humid southern mountains of Taiwan served as a
glacial refugium for numerous species, including \emph{Castanopsis
carlesii} (Cheng, Hwang, \& Lin, 2005), \emph{Quercus glauca} (Huang,
Hwang, \& Lin, 2002; Lee et al., 2006), \emph{Cinnamomum kanehirae}(Liao
et al., 2010), and \emph{Machilus kusanoi} (Wu et al., 2006), resulting
in high genetic diversity.

During the refuge period, hybridization between closely related species
may have occurred. A previous study showed that the Pleistocene glacial
oscillations accelerated \emph{Scutellaria} species interbreeding in
southern Taiwan and that current environmental heterogeneity is
responsible for shaping the genetic structure of
\emph{Scutellaria}populations (Huang et al., 2017). Due to this genetic
legacy, glacial refugia in southern Taiwan contributed to genetic
diversity in extant populations. Genetic mixing in refugia in the south
may have been delayed until the postglacial period (ca. 7-4 kya from
ABC). This inference is supported by the fact that demography explains
more than half of the adaptive genetic variation, whereas the current
climate explains only a small share. This result suggests that
demographic history, such as changes in population size and gene flow,
was crucial in shaping the genetic structure of the current
population/species.

The spread of beneficial genes may broaden the ecological niche or
enable hybridized species to colonize novel habitats (Suarez-Gonzalez et
al., 2018). We identified two adaptive genes involved in the
introgression between \emph{S. taiwanensis} and each of the potential
hybridized species: \emph{PHYB} (between \emph{S. taiwanensis}
and\emph{S. indica} ) and \emph{OPCL1} (between \emph{S. taiwanensis}
and\emph{S. playfairii} ). \emph{PHYB} was associated with the mean
diurnal range (BIO2), and \emph{S. indica} was present within a broad
range of BIO2, suggesting that this species can acclimate to a diverse
distribution of habitats and tolerate large temperature variations
within a year. \emph{PHYB} has been shown to promote plant adaptation to
cold environments (Jiang et al., 2020; Zhao et al., 2016). The
introgression of \emph{PHYB} between \emph{S. indica} and \emph{S.
taiwanensis} may have enabled \emph{S. taiwanensis} to colonize and
thrive in southern Taiwan. \emph{OPCL1} was associated with several
temperature and precipitation-related factors, including the mean
diurnal range (BIO2), temperature seasonality (BIO4), maximum
temperature of the warmest month (BIO5), minimum temperature of the
coldest month (BIO6), temperature annual range (BIO17), and
precipitation of the warmest quarter (BIO18). Among these environmental
variables, variations in BIO5 were larger in the habitats of \emph{S.
playfairii} than in the habitats of \emph{S. taiwanensis} . The
transmission of alleles of \emph{OPCL1} may further increase the
adaptability of \emph{S. taiwanensis} to lower temperatures in warm
seasons compared to most \emph{Scutellaria} species in Taiwan except
for\emph{S. hsiehii} . However, the niches of \emph{S. taiwanensis}
differed from those of either potential hybridized species and were
characterized by lower BIO4, higher BIO6, and higher BIO18, indicating a
shift toward temperature- and precipitation-related traits associated
with tropical climates in southern Taiwan. Alternatively, these
differences could indicate relics of ancestral traits of ghost lineages
that are not related to the hybridization events. Overall, we infer that
climate drivers of the introgression of adaptive alleles combined with
demographic history associated with climate oscillation generated the
unique genetic components and high genetic diversity of \emph{S.
taiwanensis} .

\textbf{Origin of \emph{S. hsiehii} through asymmetric hybridization
between \emph{S. austrotaiwanensis} and \emph{S. indica}}

The most plausible explanation for the origin of \emph{S. hsiehii} is an
asymmetric hybridization event between \emph{S. austrotaiwanensis}
and\emph{S. indica} . In this scenario, \emph{S.
austrotaiwanensis}contributed a significant proportion (91.2\%) of the
genetic makeup of\emph{S. hsiehii} . The strong asymmetry of the genetic
contributions from \emph{S. austrotaiwanensis} and \emph{S. indica}
supports the second-best hypothesis, namely, that \emph{S. hsiehii} is
derived from\emph{S. austrotaiwanensis} , despite the relatively small
likelihood difference with the best hybridization scenario.

Although the contribution of \emph{S. indica} is minor, the STRUCTURE
results indicated that it is genetically inseparable from the formation
of \emph{S. hsiehii} . Additionally, the distribution range of \emph{S.
indica} extends across a wide range of altitudes in central Taiwan, up
to 2400 m (Hsieh \& Huang, 1995), which includes the area where \emph{S.
hsiehii} is found (specifically at 1600 m in central Taiwan (Hsieh,
2013)). This observation is further evidence that \emph{S. indica} may
be a parental species.

An analysis of Grinnellian niches showed that \emph{S. hsiehii} had
different temperature preferences than \emph{S. indica} and \emph{S.
austrotaiwanensis} . Specifically, \emph{S. hsiehii} was associated with
a higher mean diurnal range (BIO2) and lower maximum temperature of the
warmest month (BIO5) and minimum temperature of the coldest month (BIO6)
than its parental species (Fig. 3c). This finding suggests a shift in
the temperature-related adaptive traits of \emph{S. hsiehii} compared to
its parental species following the hybridization event.

Although we estimate that 91.2\% of the genetic composition of \emph{S.
hsiehii} originated from \emph{S. austrotaiwanensis} , its distribution
does not align with that of \emph{S. austrotaiwanensis} , which is
typically found in mountainous areas at altitudes of 400-500 m in
southern Taiwan (Hsieh \& Huang, 1997). This discrepancy indicates a
shift in the distribution of \emph{S. hsiehii} to higher elevations
compared to its parent species \emph{S. austrotaiwanensis} . A previous
study of \emph{S. hsiehii} (Hsieh, 2013) compared this species with
another species of the \emph{indica} group, \emph{S. playfairii} , and
stated that the floral parts of the two species were morphologically
similar; however, this study did not compare \emph{S. hsiehii}
with\emph{S. austrotaiwanensis} . We found that \emph{S. hsiehii} most
closely resembled \emph{S. austrotaiwanensis} , although the long-ovate
leaved form of \emph{S. hsiehii} differs greatly from those of the other
Taiwanese \emph{Scutellaria} species (oval to obovate leaves).
Although\emph{S. austrotaiwanensis} is primarily found in the south and
is geographically isolated from \emph{S. hsiehii} in central Taiwan,
similar to \emph{S. taiwanensis} , it is possible that this
\emph{indica}group species also retreated to the southern glacial
refugium, where the hybridization event may have occurred during or
after the glacial period. However, the substantiation of these
conclusions requires further evidence.

\textbf{Adaptive introgression from \emph{S. tashiroi} to \emph{S.
taiwanensis}, \emph{S. indica}, and \emph{S. barbata}}

Genetic introgression may diminish the species boundary but may transfer
adaptive genes into another species. The multiple introgression events
among Taiwanese \emph{Scutellaria} species indicate frequent genetic
replacement. Specifically, we found evidence of ancient gene flow
between \emph{S. hsiehii} and \emph{S. tashiroi} , \emph{S.
playfairii}and \emph{S. taiwanensis} , and the ancestors of \emph{S.
barbata} and the \emph{indica} group (Fig. 3B). Although some of the
analyses were incongruent, these findings provide insights into
occasional interbreeding that occurred during the evolutionary history
of these species.

The results revealed adaptive introgression involving three specific
genes (\emph{LYC} , \emph{OPCL1} , and \emph{PHYB} ) responsible for
antioxidant biosynthesis, pathogen-resistance processes, and
light-promotion responses (Fig. 3A, Table S6). These genes are linked to
the response to environmental stress, which could indirectly affect
species survival, growth, and overall fitness. Notably, the adaptive
introgression of these genes was correlated with BIO2 (mean diurnal
range of temperature) and BIO18 (precipitation of the warmest quarter)
(Table 1). These correlations suggest that variations in temperature and
rainfall impacted populations across altitudes within a narrow latitude
range of the island, creating diverse stress environments. Introducing
genetic material through introgression may have benefited species
survival under these varying stress conditions. Moreover, the rainfall
during the warmest season after spring aligns with the growth season
after germination of these \emph{Scutellaria} species in Taiwan, which
may influence their growth rate. However, the irregularity of BIO2 and
BIO18 between the source and sink species, along with other associated
bioclimatic variables (BIO4, BIO5, BIO6, and BIO17) (Fig. 3C), implies
that these variables are associated with other genetic elements through
hitchhiking.

This study shows that stress-related genes were adapted and transferred
between species in response to environmental changes, significantly
impacting their survival, growth, and fitness. The correlations of
introgression with temperature and rainfall factors likely aided species
survival across altitude variations. The exchange of stress-related
genes between species probably enabled them to tolerate or resist
adverse conditions, leading to a convergence in adaptive responses to
environmental cues and light/shade tolerance.

\textbf{Conclusions}

This study provides comprehensive insights into the complex evolutionary
dynamics of Taiwanese \emph{Scutellaria} species. A combination of
factors, including hybridization, introgression, and climate-driven
selection, influenced their origins and adaptations. The intricate
interplay between \emph{ex situ} and \emph{in situ} speciation
processes, driven by allopatric and sympatric modes, shaped the
diversity within the island. Frequent hybridization events contributed
to genetic diversity and facilitated \emph{in situ} speciation,
broadening ecological potential and enhancing adaptation to diverse
habitats. Local climate heterogeneity also played a pivotal role,
allowing genetically distinct yet ecologically similar species to
hybridize and adapt to changing environments. The hybridization events
occasionally blurred species boundaries but simultaneously transferred
adaptive genes, further aiding survival and adaptation. The genetic
composition of these species reflects a complex interplay of historical
demographic events and climate factors. By investigating the genetic
basis of adaptation, hybridization, and introgression, this study
highlights the complex mechanisms driving the evolution
of\emph{Scutellaria} species within the unique ecological context of
Taiwan.

\textbf{Acknowledgments}

The authors express their gratitude to Dr. Jianhuo Fang for his
assistance with the NGS sequencing platform and library preparation.
Special thanks to Dr. Bing-Hong Huang for his invaluable contributions
to sample collection and experimental work. The authors also thank Dawn
Schmidt for her English editing.

\textbf{References}

Abbott, R., Albach, D., Ansell, S., Arntzen, J. W., Baird, S. J.,
Bierne, N., . . . Buggs, R. (2013). Hybridization and
speciation.\emph{Journal of Evolutionary Biology, 26} (2), 229-246.

Ahrens, C. W., Rymer, P. D., Stow, A., Bragg, J., Dillon, S., Umbers, K.
D., \& Dudaniec, R. Y. (2018). The search for loci under selection:
trends, biases and progress. \emph{Molecular ecology, 27} (6),
1342-1356.

Algar, A. C., \& Mahler, D. L. (2016). Area, climate heterogeneity, and
the response of climate niches to ecological opportunity in island
radiations of \emph{Anolis} lizards. \emph{Global Ecology and
Biogeography, 25} (7), 781-791.
doi:\emph{https://doi.org/10.1111/geb.12327}

Arnold, M. L. (1997). \emph{Natural hybridization and evolution} :
Oxford University Press.

Barrett, R. D. H., \& Schluter, D. (2008). Adaptation from standing
genetic variation. \emph{Trends in Ecology \& Evolution, 23} (1), 38-44.
doi:\emph{https://doi.org/10.1016/j.tree.2007.09.008}

Bock, D. G., Cai, Z., Elphinstone, C., González-Segovia, E.,
Hirabayashi, K., Huang, K., . . . Rieseberg, L. H. (2023). Genomics of
plant speciation. \emph{Plant Communications} , 100599.
doi:\emph{https://doi.org/10.1016/j.xplc.2023.100599}

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

Brauer, C. J., Sandoval-Castillo, J., Gates, K., Hammer, M. P., Unmack,
P. J., Bernatchez, L., \& Beheregaray, L. B. (2023). Natural
hybridization reduces vulnerability to climate change. \emph{Nature
Climate Change, 13} (3), 282-289. doi:10.1038/s41558-022-01585-1

Brochmann, C., Borgen, L., \& Stabbetorp, O. E. (2000). Multiple diploid
hybrid speciation of the Canary Island endemic \emph{Argyranthemum
sundingii} (Asteraceae). \emph{Plant Systematics and Evolution, 220}
(1), 77-92. doi:10.1007/BF00985372

Bu, D., Luo, H., Huo, P., Wang, Z., Zhang, S., He, Z., . . . Guo, J.
(2021). KOBAS-i: intelligent prioritization and exploratory
visualization of biological functions for gene enrichment
analysis.\emph{Nucleic acids research, 49} (W1), W317-W325.

Buckley, T. R., Cordeiro, M., Marshall, D. C., \& Simon, C. (2006).
Differentiating between hypotheses of lineage sorting and introgression
in New Zealand alpine cicadas (Maoricicada Dugdale). \emph{Systematic
biology, 55} (3), 411-425.

Capblancq, T., Lachmuth, S., Fitzpatrick, M. C., \& Keller, S. R.
(2023). From common gardens to candidate genes: exploring local
adaptation to climate in red spruce. \emph{New Phytologist, 237} (5),
1590-1605.

Capblancq, T., Luu, K., Blum, M. G., \& Bazin, E. (2018). Evaluation of
redundancy analysis to identify signatures of local
adaptation.\emph{Molecular ecology resources, 18} (6), 1223-1233.

Capella-Gutiérrez, S., Silla-Martínez, J. M., \& Gabaldón, T. (2009).
trimAl: a tool for automated alignment trimming in large-scale
phylogenetic analyses. \emph{Bioinformatics, 25} (15), 1972-1973.

Cha, J.-Y., Kim, W.-Y., Kang, S. B., Kim, J. I., Baek, D., Jung, I. J.,
. . . Nakajima, M. (2015). A novel thiol-reductase activity of
Arabidopsis YUC6 confers drought tolerance independently of auxin
biosynthesis. \emph{Nature communications, 6} (1), 8041.

Chao, C.-T., Huang, B.-H., Chang, J.-T., \& Liao, P.-C. (2020).
Taxonomic notes on \emph{Scutellaria taipeiensis} (Lamiaceae) from
morphological and molecular data. \emph{PhytoKeys, 140} , 33-45.

Charng, Y.-y., Liu, H.-c., Liu, N.-y., Hsu, F.-c., \& Ko, S.-s. (2006).
Arabidopsis Hsa32, a novel heat shock protein, is essential for acquired
thermotolerance during long recovery after acclimation. \emph{Plant
Physiology, 140} (4), 1297-1305.

Cheng, Y. P., Hwang, S. Y., \& Lin, T. P. (2005). Potential refugia in
Taiwan revealed by the phylogeographical study of \emph{Castanopsis
carlesii} Hayata (Fagaceae). \emph{Molecular Ecology, 14} (7),
2075-2085. doi:10.1111/j.1365-294X.2005.02567.x

Chhatre, V. E., \& Emerson, K. J. (2017). StrAuto: automation and
parallelization of STRUCTURE analysis. \emph{BMC bioinformatics, 18}
(1), 1-5.

Chiang, T.-Y., \& Schaal, B. A. (2006). Phylogeography of plants in
Taiwan and the Ryukyu Archipelago. \emph{TAXON, 55} (1), 31-41.

Chiang, Y.-C., Huang, B.-H., \& Liao, P.-C. (2012). Diversification,
biogeographic pattern, and demographic history of Taiwanese Scutellaria
species inferred from nuclear and chloroplast DNA. \emph{PLoS One, 7}
(11), e50844.

Cowie, R. H. (1995). Variation in species diversity and shell shape in
Hawaiian land snails: \emph{In situ} speciation and ecological
relationships. \emph{Evolution, 49} (6), 1191-1202.
doi:10.1111/j.1558-5646.1995.tb04446.x

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.

de Sousa, F., Bertrand, Y. J., Nylinder, S., Oxelman, B., Eriksson, J.
S., \& Pfeil, B. E. (2014). Phylogenetic properties of 50 nuclear loci
in Medicago (Leguminosae) generated using multiplexed sequence capture
and next-generation sequencing. \emph{PLoS One, 9} (10), e109704.

Degnan, J. H., \& Rosenberg, N. A. (2009). Gene tree discordance,
phylogenetic inference and the multispecies coalescent. \emph{Trends in
Ecology \& Evolution, 24} (6), 332-340.
doi:\emph{https://doi.org/10.1016/j.tree.2009.01.009}

Du, F. K., Petit, R. J., \& Liu, J. Q. (2009). More introgression with
less gene flow: chloroplast vs. mitochondrial DNA in the Picea asperata
complex in China, and comparison with other conifers. \emph{Molecular
ecology, 18} (7), 1396-1407.

Du, F. K., Wang, T., Wang, Y., Ueno, S., \& de Lafontaine, G. (2020).
Contrasted patterns of local adaptation to climate change across the
range of an evergreen oak, Quercus aquifolioides. \emph{Evolutionary
Applications, 13} (9), 2377-2391.

Edelman, N. B., Frandsen, P. B., Miyagi, M., Clavijo, B., Davey, J.,
Dikow, R. B., . . . Neafsey, D. E. (2019). Genomic architecture and
introgression shape a butterfly radiation. \emph{Science, 366} (6465),
594-599.

Ellepola, G., Pie, M. R., Pethiyagoda, R., Hanken, J., \& Meegaskumbura,
M. (2022). The role of climate and islands in species diversification
and reproductive-mode evolution of Old World tree
frogs.\emph{Communications Biology, 5} (1), 347.

Emerson, B. C., \& Patino, J. (2018). Anagenesis, cladogenesis, and
speciation on islands. \emph{Trends in Ecology \& Evolution, 33} (7),
488-491. doi:10.1016/j.tree.2018.04.006

Excoffier, L., \& Lischer, H. E. (2010). Arlequin suite ver 3.5: a new
series of programs to perform population genetics analyses under Linux
and Windows. \emph{Molecular ecology resources, 10} (3), 564-567.

Excoffier, L., Marchi, N., Marques, D. A., Matthey-Doret, R., Gouy, A.,
\& Sousa, V. C. (2021). fastsimcoal2: demographic inference under
complex evolutionary scenarios. \emph{Bioinformatics, 37} (24),
4882-4885.

Exposito-Alonso, M. (2023). Understanding local plant extinctions before
it is too late: bridging evolutionary genomics with global
ecology.\emph{New Phytologist, 237} (6), 2005-2011.
doi:\emph{https://doi.org/10.1111/nph.18718}

Feng, L., \& Du, F. K. (2022). Landscape genomics in tree conservation
under a changing environment. \emph{Frontiers in Plant Science, 13} ,
822217.

Fick, S. E., \& Hijmans, R. J. (2017). WorldClim 2: new 1\selectlanguage{english}-km spatial
resolution climate surfaces for global land areas. \emph{International
journal of climatology, 37} (12), 4302-4315.

Fitak, R. R. (2021). OptM: estimating the optimal number of migration
edges on population trees using Treemix. \emph{Biology Methods and
Protocols, 6} (1), bpab017.

Foll, M. (2012). BayeScan v2. 1 user manual. \emph{Ecology, 20} (10),
1450-1462.

Forester, B. R., Lasky, J. R., Wagner, H. H., \& Urban, D. L. (2018).
Comparing methods for detecting multilocus adaptation with multivariate
genotype--environment associations. \emph{Molecular ecology, 27} (9),
2215-2233.

Funk, D. J., \& Omland, K. E. (2003). Species-level paraphyly and
polyphyly: frequency, causes, and consequences, with insights from
animal mitochondrial DNA. \emph{Annual Review of Ecology, Evolution, and
Systematics, 34} (1), 397-423.

Grant, P. R., \& Grant, B. R. (2019). Hybridization increases population
variation during adaptive radiation. \emph{Proceedings of the National
Academy of Sciences, 116} (46), 23216-23224. doi:10.1073/pnas.1913534116

Gunther, T., \& Coop, G. (2013). Robust identification of local
adaptation from allele frequencies. \emph{Genetics, 195} (1), 205-220.

Harrison, R. G., \& Larson, E. L. (2014). Hybridization, introgression,
and the nature of species boundaries. \emph{Journal of Heredity, 105}
(S1), 795-809. doi:10.1093/jhered/esu033

He, D., Damaris, R. N., Fu, J., Tu, J., Fu, T., Xi, C., . . . Yang, P.
(2016). Differential molecular responses of rapeseed cotyledons to light
and dark reveal metabolic adaptations toward autotrophy
establishment.\emph{Frontiers in Plant Science, 7} , 988.

Hedrick, P. W. (2013). Adaptive introgression in animals: examples and
comparison to new mutation and standing variation as sources of adaptive
variation. \emph{Molecular ecology, 22} (18), 4606-4618.

Hsieh, T.-H. (2013). \emph{Scutellaria hsiehii} (Lamiaceae), a New
Species from Taiwan. \emph{58} (4), 242-245.

Hsieh, T.-H., \& Huang, T.-C. (1995). Notes on the flora of Taiwan
(20)-- \emph{Scutellaria} (Lamiaceae) in Taiwan. \emph{Taiwania, 40}
(1), 35-56.

Hsieh, T.-H., \& Huang, T.-C. (1997). Notes on the Flora of Taiwan
(29)\selectlanguage{english}-\emph{Scutellaria austrotaiwanensis} Hsieh \& Huang \emph{sp.
nov.} (Lamiaceae) from Taiwan. \emph{42} (2), 109-116.

Hsiung, H.-Y., Huang, B.-H., Chang, J.-T., Huang, Y.-M., Huang, C.-W.,
\& Liao, P.-C. (2017). Local climate heterogeneity shapes population
genetic structure of two undifferentiated insular
\emph{Scutellaria}species. \emph{Frontiers in Plant Science, 8} , 159.

Huang, B.-H., Chen, Y.-W., Huang, C.-L., Gao, J., \& Liao, P.-C.
(2016a). Diversifying selection of the anthocyanin biosynthetic
downstream gene \emph{UFGT} accelerates floral diversity of
island\emph{Scutellaria} species. \emph{BMC Evolutionary Biology, 16}
(1), 191. doi:10.1186/s12862-016-0759-0

Huang, B.-H., Chen, Y.-W., Huang, C.-L., Gao, J., \& Liao, P.-C.
(2016b). Imbalanced positive selection maintains the functional
divergence of duplicated \emph{DIHYDROKAEMPFEROL 4-REDUCTASE}
genes.\emph{Scientific Reports, 6} (1), 39031. doi:10.1038/srep39031

Huang, B.-H., Huang, C.-W., Huang, C.-L., \& Liao, P.-C. (2017).
Continuation of the genetic divergence of ecological speciation by
spatial environmental heterogeneity in island endemic
plants.\emph{Scientific Reports, 7} (1), 5465.
doi:10.1038/s41598-017-05900-1

Huang, B.-H., Pang, E., Chen, Y.-W., Cao, H., Ruan, Y., \& Liao, P.-C.
(2015). Positive selection and functional divergence of R2R3-MYB
paralogous genes expressed in inflorescence buds of
\emph{Scutellaria}species (Labiatae). \emph{International Journal of
Molecular Sciences, 16} (3), 5900-5921. Retrieved from
doi:10.3390/ijms16035900

Huang, S. S. F., Hwang, S.-Y., \& Lin, T.-P. (2002). Spatial pattern of
chloroplast DNA variation of \emph{Cyclobalanopsis glauca} in Taiwan and
East Asia. \emph{Molecular Ecology, 11} (11), 2349-2358.
doi:10.1046/j.1365-294X.2002.01624.x

Huson, D. H., \& Bryant, D. (2006). Application of phylogenetic networks
in evolutionary studies. \emph{Molecular biology and evolution, 23} (2),
254-267.

Jiang, B., Shi, Y., Peng, Y., Jia, Y., Yan, Y., Dong, X., . . . Gong, Z.
(2020). Cold-induced CBF--PIF3 interaction enhances freezing tolerance
by stabilizing the phyB thermosensor in Arabidopsis. \emph{Molecular
plant, 13} (6), 894-906.

Jiang, L., Bao, Q., He, W., Fan, D. M., Cheng, S. M., Lopez\selectlanguage{english}-Pujol, J., .
. . Gedik, A. (2022). Phylogeny and biogeography of Fagus (Fagaceae)
based on 28 nuclear single/low\selectlanguage{english}-copy loci. \emph{Journal of Systematics
and Evolution, 60} (4), 759-772.

Jombart, T. (2008). adegenet: a R package for the multivariate analysis
of genetic markers. \emph{Bioinformatics, 24} (11), 1403-1405.

Kamvar, Z. N., Tabima, J. F., \& Grunwald, N. J. (2014). Poppr: an R
package for genetic analysis of populations with clonal, partially
clonal, and/or sexual reproduction. \emph{PeerJ, 2} , e281.

Katoh, K., \& Standley, D. M. (2013). MAFFT multiple sequence alignment
software version 7: improvements in performance and
usability.\emph{Molecular biology and evolution, 30} (4), 772-780.

Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., \&
Mayrose, I. (2015). Clumpak: a program for identifying clustering modes
and packaging population structure inferences across K. \emph{Molecular
ecology resources, 15} (5), 1179-1191.

Lasky, J. R., Josephs, E. B., \& Morris, G. P. (2023).
Genotype--environment associations to reveal the molecular basis of
environmental adaptation. \emph{The Plant Cell, 35} (1), 125-138.
doi:10.1093/plcell/koac267

Lee, Y., Jr., Hwang, S.-Y., Ho, K.-C., \& Lin, T.-P. (2006). Source
populations of \emph{Quercus glauca} in the Last Glacial Age in Taiwan
revealed by Nuclear microsatellite markers. \emph{Journal of Heredity,
97} (3), 261-269. doi:10.1093/jhered/esj030

Lev-Yadun, S., \& Gould, K. S. (2009). Role of Anthocyanins in Plant
Defence. In C. Winefield, K. Davies, \& K. Gould
(Eds.),\emph{Anthocyanins: Biosynthesis, Functions, and Applications}
(pp. 22-28). New York, NY: Springer New York.

Levin, R. A., Whelan, A., \& Miller, J. S. (2009). The utility of
nuclear conserved ortholog set II (COSII) genomic regions for
species-level phylogenetic inference in Lycium
(Solanaceae).\emph{Molecular Phylogenetics and Evolution, 53} (3),
881-890.

Li, H. (2011). A statistical framework for SNP calling, mutation
discovery, association mapping and population genetical parameter
estimation from sequencing data. \emph{Bioinformatics, 27} (21),
2987-2993.

Li, H. (2013). Aligning sequence reads, clone sequences and assembly
contigs with BWA-MEM. \emph{arXiv preprint arXiv:1303.3997} .

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

Li, Z., De La Torre, A. R., Sterck, L., Canovas, F. M., Avila, C.,
Merino, I., . . . Van de Peer, Y. (2017). Single-copy genes as molecular
markers for phylogenomic studies in seed plants. \emph{Genome Biology
and Evolution, 9} (5), 1130-1147.

Liao, P.-C., Kuo, D.-C., Lin, C.-C., Ho, K.-C., Lin, T.-P., \& Hwang,
S.-Y. (2010). Historical spatial range expansion and a very recent
bottleneck of \emph{Cinnamomum kanehirae} Hay.(Lauraceae) in Taiwan
inferred from nuclear genes. \emph{BMC Evolutionary Biology, 10} (1),
124.

Lin, C. C. (1963). Geology and ecology of Taiwan prehistory. \emph{Asian
Perspectives, 7} (1/2), 203-213.

Malinsky, M., Matschiner, M., \& Svardal, H. (2021). Dsuite\selectlanguage{english}-Fast
D\selectlanguage{english}-statistics and related admixture evidence from VCF
files.\emph{Molecular ecology resources, 21} (2), 584-595.

Mao, J., Zhang, Y.-C., Sang, Y., Li, Q.-H., \& Yang, H.-Q. (2005). A
role for Arabidopsis cryptochromes and COP1 in the regulation of
stomatal opening. \emph{Proceedings of the National Academy of Sciences,
102} (34), 12270-12275.

Noguerales, V., \& Ortego, J. (2022). Genomic evidence of speciation by
fusion in a recent radiation of grasshoppers. \emph{Evolution, 76} (11),
2618-2633. doi:10.1111/evo.14508

Oksanen, J., Kindt, R., Legendre, P., O'Hara, B., Stevens, M. H. H.,
Oksanen, M. J., \& Suggests, M. (2007). The vegan
package.\emph{Community ecology package, 10} (631-637), 719.

Olave, M., \& Meyer, A. (2020). Implementing large genomic single
nucleotide polymorphism data sets in phylogenetic network
reconstructions: a case study of particularly rapid radiations of
cichlid fish. \emph{Systematic biology, 69} (5), 848-862.

Pickrell, J., \& Pritchard, J. (2012). Inference of population splits
and mixtures from genome-wide allele frequency data. \emph{Nature
Precedings} , 1-1.

Porebski, S., Bailey, L. G., \& Baum, B. R. (1997). Modification of a
CTAB DNA extraction protocol for plants containing high polysaccharide
and polyphenol components. \emph{Plant molecular biology reporter, 15} ,
8-15.

Porretta, D., \& Canestrelli, D. (2023). The ecological importance of
hybridization. \emph{Trends in Ecology \& Evolution} .
doi:10.1016/j.tree.2023.07.003

R Core Team, R. (2013). R: A language and environment for statistical
computing.

Rieseberg, L. H., \& Wendel, J. F. (1993). Introgression and its
consequences in plants. \emph{Hybrid zones and the evolutionary process,
70} , 109.

Rosche, C., Baasch, A., Runge, K., Brade, P., Trager, S., Parisod, C.,
\& Hensen, I. (2022). Tracking population genetic signatures of local
extinction with herbarium specimens. \emph{Annals of Botany, 129} (7),
857-868. doi:10.1093/aob/mcac061

Sang, Y., Long, Z., Dan, X., Feng, J., Shi, T., Jia, C., . . . Zhang, H.
(2022). Genomic insights into local adaptation and future
climate-induced vulnerability of a keystone forest tree in East
Asia.\emph{Nature communications, 13} (1), 6541.

Solis-Lemus, C., Bastide, P., \& Ane, C. (2017). PhyloNetworks: a
package for phylogenetic networks. \emph{Molecular biology and
evolution, 34} (12), 3292-3298.

Stelkens, R., \& Bendixsen, D. P. (2022). The evolutionary and
ecological potential of yeast hybrids. \emph{Current Opinion in Genetics
\& Development, 76} , 101958.
doi:\emph{https://doi.org/10.1016/j.gde.2022.101958}

Stuessy, T. F., Jakubowsky, G., Gomez, R. S., Pfosser, M., Schluter, P.
M., Fer, T., . . . Kato, H. (2006). Anagenetic evolution in island
plants. \emph{Journal of Biogeography, 33} (7), 1259-1265.
doi:10.1111/j.1365-2699.2006.01504.x

Suarez-Gonzalez, A., Hefer, C. A., Lexer, C., Douglas, C. J., \& Cronk,
Q. C. B. (2018). Introgression from \emph{Populus balsamifera} underlies
adaptively significant variation and range boundaries
in\emph{P.~trichocarpa} . \emph{New Phytologist, 217} (1), 416-427.
doi:\emph{https://doi.org/10.1111/nph.14779}

Suzuki, N., Bassil, E., Hamilton, J. S., Inupakutika, M. A., Zandalinas,
S. I., Tripathy, D., . . . Kumazaki, A. (2016). ABA is required for
plant acclimation to a combination of salt and heat stress. \emph{PLoS
One, 11} (1), e0147625.

Trabucco, A., \& Zomer, R. J. (2018). Global aridity index and potential
evapotranspiration (ET0) climate database v2. \emph{CGIAR Consort Spat
Inf, 10} , m9.

Untergasser, A., Cutcutache, I., Koressaar, T., Ye, J., Faircloth, B.
C., Remm, M., \& Rozen, S. G. (2012). Primer3---new capabilities and
interfaces. \emph{Nucleic acids research, 40} (15), e115-e115.

Villa\selectlanguage{english}-Machio, I., Heuertz, M., Alvarez, I., \& Nieto Feliner, G. (2023).
Demography\selectlanguage{english}-driven and adaptive introgression in a hybrid zone of the
Armeria syngameon. \emph{Molecular ecology} .

Wiens, B. J., Combe, F. J., Dickerson, B., Divine, L. M., Padula, V. M.,
Sage, G. K., . . . Hope, A. G. (2022). Genetic drift drives rapid
speciation of an Arctic insular endemic shrew (\emph{Sorex
pribilofensis} ). \emph{Molecular Ecology, 31} (20), 5231-5248.
doi:\emph{https://doi.org/10.1111/mec.16658}

Wu, S.-H., Hwang, C.-Y., Lin, T.-P., Chung, J.-D., Cheng, Y.-P., \&
Hwang, S.-Y. (2006). Contrasting phylogeographical patterns of two
closely related species, \emph{Machilus thunbergii} and \emph{Machilus
kusanoi} (Lauraceae), in Taiwan. \emph{Journal of Biogeography, 33} (5),
936-947. doi:10.1111/j.1365-2699.2006.01431.x

Ye, J. W., Yang, Z. Z., \& Tian, B. (2023). Tempo\selectlanguage{english}-spatial evolution of
seed plant endemism in Taiwan island. \emph{Journal of Biogeography} .

Yu, H.-S. (2003). Geological Characteristics and Distribution of
Submarine Physiographic Features in the Taiwan Region. \emph{Marine
Georesources \& Geotechnology, 21} (3-4), 139-153. doi:10.1080/713773391

Zecca, G., Labra, M., \& Grassi, F. (2020). Untangling the evolution of
American wild grapes: Admixed species and how to find
them.\emph{Frontiers in Plant Science, 10} , 1814.

Zhao, C., Zhang, Z., Xie, S., Si, T., Li, Y., \& Zhu, J.-K. (2016).
Mutational evidence for the critical role of CBF transcription factors
in cold acclimation in Arabidopsis. \emph{Plant Physiology, 171} (4),
2744-2759.

Zhao, J., Zhang, W., Zhao, Y., Gong, X., Guo, L., Zhu, G., . . . Guo, Y.
(2007). SAD2, an importin \selectlanguage{greek}β-\selectlanguage{english}like protein, is required for UV-B response
in Arabidopsis by mediating MYB4 nuclear trafficking. \emph{The Plant
Cell, 19} (11), 3805-3818.

Zhao, Q., Yang, J., Cui, M.-Y., Liu, J., Fang, Y., Yan, M., . . .
Yidiresi, R. (2019). The reference genome sequence of Scutellaria
baicalensis provides insights into the evolution of wogonin
biosynthesis. \emph{Molecular plant, 12} (7), 935-950.

\textbf{Data Accessibility}

The following primary datasets associated with this paper have been
deposited on Figshare under the DOI: 10.6084/m9.figshare.24314641:

1. Variant Call Format (VCF) files encompassing all individuals
subjected to analysis in this study.

2. Parameters and scenarios employed in all Approximate Bayesian
Computation (ABC) simulations.

3. Environmental data corresponding to each population.

\textbf{Funding}

This work was supported by the National Science and Technology Council
under project ID NSTC 112-2621-B-003-001-MY3.

\textbf{Competing Interests}

The authors have no relevant financial or non-financial interests to
disclose.

\textbf{Author contribution}

F.K.D. and P.C.L. designed and supervised the project. P.W.S. performed
laboratory experiments and statistical analyses. P.W.S., P.C.L., and
D.F.K. drafted the manuscript. J.T.C., M.X.L., and C.T.C. assisted with
data interpretation, provided valuable feedback, and revised the
manuscript. All authors have read and approved the final version of the
manuscript.

\textbf{Tables and Figures}

Table 1. Genes exhibiting signals of introgression and scanned by GEA
procedures.\selectlanguage{english}
\begin{longtable}[]{@{}lllll@{}}
\toprule
Pop1 & Pop2 & Gene ID & Gene name & Variable\tabularnewline
TAI & (PLA,(HSI,AUS)) & Sb01g09340 & LYC & BIO2,BIO18\tabularnewline
TAI & PLA & Sb05g13050 & OPCL1 &
BIO2,BIO4,BIO5,BIO6,BIO17,BIO18\tabularnewline
TAS & BAR & Sb05g13050 & OPCL1 &
BIO2,BIO4,BIO5,BIO6,BIO17,BIO18\tabularnewline
TAS & TAI & Sb05g13050 & OPCL1 &
BIO2,BIO4,BIO5,BIO6,BIO17,BIO18\tabularnewline
TAS & IND & Sb05g13050 & OPCL1 &
BIO2,BIO4,BIO5,BIO6,BIO17,BIO18\tabularnewline
TAS & BAR & Sb09g00180 & PHYB & BIO2\tabularnewline
TAS & TAI & Sb09g00180 & PHYB & BIO2\tabularnewline
TAS & IND & Sb09g00180 & PHYB & BIO2\tabularnewline
TAI & IND & Sb09g00180 & PHYB & BIO2\tabularnewline
TAI & HSI & Sb09g00180 & PHYB & BIO2\tabularnewline
\bottomrule
\end{longtable}

Table 2. The influences of climate, geography, and demography on genetic
variation as decomposed by pRDA. *p\textless{}0.05, **p\textless{}0.01,
***p\textless{}0.001.\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
Partial RDA models & P (\textgreater{}F) & Proportion of explainable
variation & Proportion of total variation\tabularnewline
Full model & 0.001 *** & 1 & 0.711\tabularnewline
Pure climate & 0.383 & 0.011 & 0.008\tabularnewline
Pure geography & 0.22 & 0.027 & 0.019\tabularnewline
Pure demography & 0.001 *** & 0.733 & 0.521\tabularnewline
Joint climate/geography/demography & & & 0.163\tabularnewline
Total unexplained & & & 0.289\tabularnewline
\bottomrule
\end{longtable}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-1-neworder/fig-1-neworder}
\caption{{Sampling locations and genetic differentiation among the investigated
individuals. (A) Sampling locations with population assignment with K =
9. We collected 26 populations encompassing 7 species across Taiwan,
including 4 sympatric populations and a \emph{S. tashiroi} population
collected from Orchid Islet. (B) PCA of all individuals in the space of
the first two principal components and in the space of the second and
the third principal components. (C) Assignments of each of the 504
individuals by the most-supported models with K = 2 and K = 9 genetic
clusters. TAI, \emph{S. taiwanensis} ; IND,\emph{S. indica} ; PLA,
\emph{S. playfairii} ; AUS, \emph{S. austrotaiwanensis} ; HSI, \emph{S.
hsiehii} ; TAS, \emph{S. tashiroi} ; BAR, \emph{S. barbata~}.
{\label{396037}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-2-new-0913/fig-2-new-0913}
\caption{{Distance-based phylogenetic network implemented in SplitsTree 6. TAI,
IND, BAR, and TAS are each separated into two intraspecies lineages that
exhibit reticulated relationships with other species. TAI,\emph{S.
taiwanensis} ; IND, \emph{S. indica} ; PLA, \emph{S. playfairii} ; AUS,
\emph{S. austrotaiwanensis} ; HSI, \emph{S. hsiehii} ; TAS, \emph{S.
tashiroi} ; BAR, \emph{S. barbata~}.
{\label{758243}}%
}}
\end{center}
\end{figure}

~\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-3/fig-3}
\caption{{Tests for introgression and putative adaptive loci. (A) Genes with a
significant footprint of introgression as estimated by the\emph{f}
-branch statistic. Darker colors indicate greater evidence of gene flow
between species. The dotted lines in the phylogeny represent the
ancestral lineage. (B) Tests for introgression using the all-biallelic
SNP data set. The arrows represent the direction and the species with
gene flow estimated by TreeMix (orange dotted line), Dsuite (brown
dotted line), and PhyloNetworks (blue lines). The blue numbers represent
the proportion of genes inherited via gene flow from each parent. (C)
Environmental differences between the sampling locations and species as
measured by the eight retained factors: mean diurnal range (BIO2),
temperature seasonality (BIO4), maximum temperature of the warmest month
(BIO5), minimum temperature of the coldest month (BIO6), precipitation
of the driest quarter (BIO17), precipitation of the warmest quarter
(BIO18), solar radiation in September (Solar9), and global aridity index
(Aridity). The y-axis represents the value of the associated
environmental factor. (D) Outlier SNP sites scanned with each deployed
method. The numbers represent the number of scanned and overlapping SNP
sites from each procedure.
{\label{308529}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fac/fac}
\caption{{Simulated demographic histories inferred from\emph{fastsimcoal2} . (A)
Illustration of the best demographic model (Recentghost) revealed by
fastsimcoal2 and the AIC values for TAI-IND-PLA. The model suggests the
occurrence of two introgression events with different species that
affected the genetic component of TAI during the last 8 kyr. We also
included a ghost lineage representing individuals not covered by our
sampling scheme. (B) Illustration of the best demographic model (hybrid
speciation) revealed by fastsimcoal2 and the AIC values for HSI-IND-AUS.
The effective population size (Ne) with 95\% CI of each population is
depicted on the plot. The number on the arrow represents the proportion
of inherited gene flow. TAI, \emph{S. taiwanensis} ; IND, \emph{S.
indica} ; PLA, \emph{S. playfairii} ; AUS,\emph{S. austrotaiwanensis} ;
HSI, \emph{S. hsiehii} ; TAS, \emph{S. tashiroi} ; BAR, \emph{S.
barbata~}.
{\label{284653}}%
}}
\end{center}
\end{figure}

~

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

