ORIGINAL RESEARCH article
Sec. Plant Systematics and Evolution
Volume 13 - 2022 | https://doi.org/10.3389/fpls.2022.824158
Divergence With Gene Flow and Contrasting Population Size Blur the Species Boundary in Cycas Sect. Asiorientales, as Inferred From Morphology and RAD-Seq Data
- 1School of Life Science, National Taiwan Normal University, Taipei, Taiwan
- 2Botanic Garden, Field Science Center for Northern Biosphere, Hokkaido University, Sapporo, Japan
- 3Department of Anthropology, Smithsonian Institution, National Museum of Natural History, Washington, DC, United States
The divergence process of incipient species is fascinating but elusive by incomplete lineage sorting or gene flow. Species delimitation is also challenging among those morphologically similar allopatric species, especially when lacking comprehensive data. Cycas sect. Asiorientales, comprised of C. taitungensis and C. revoluta in the Ryukyu Archipelago and Taiwan, diverged recently with continuous gene flow, resulting in a reciprocal paraphyletic relationship. Their previous evolutionary inferences are questioned from few genetic markers, incomplete sampling, and incomprehensive morphological comparison by a long-term taxonomic misconception. By whole range sampling, this study tests the geographic mode of speciation in the two species of Asiorientales by approximate Bayesian computation (ABC) using genome-wide single nucleotide polymorphisms (SNPs). The individual tree was reconstructed to delimit the species and track the gene-flow trajectory. With the comparison of diagnostic morphological traits and genetic data, the allopatric speciation was rejected. Alternatively, continuous but spatially heterogeneous gene flow driven by transoceanic vegetative dispersal and pollen flow with contrasting population sizes blurred their species boundary. On the basis of morphological, genetic, and evolutionary evidence, we synonymized these two Cycas species. This study highlights not only the importance of the Kuroshio Current to species evolution but also the disadvantage of using species with geographically structured genealogies as conservation units.
Incipient species are critical for evolutionary biologists to study speciation, but they also challenge taxonomy due to gene flow or ancestral polymorphism. The former and contrasting population size lead to larger intraspecific than interspecific variations, a phenomenon called the species-definition anomaly zone (Jiao and Yang, 2021). The latter results in inconsistent species delimitation and gene tree, called the species-tree anomaly zone (Degnan and Rosenberg, 2006). Because speciation is a continuous process, various species concepts that impose discrete (sub)species status also confuses taxonomic status. Currently, the unified species concept (USC), one of 34 species concepts (Zachos, 2018), integrates species as separately evolving metapopulations and considers other species concepts as operational criteria to mitigate the problem (but also see Freudenstein et al. (2017), for the “phenophyletic view” of the species concept).
Cycas L. is a taxonomically complex group because of the recent divergence of extant taxa in the Paleocene, morpho-ecological conservatism, and the fragmented geographic distribution (Norstog and Nicholls, 1997; Nagalingum et al., 2011; Xiao and Moller, 2015; Zheng et al., 2017; Habib et al., 2021; Liu et al., 2021b). These factors make Cycas a sound study system for the geographic mode of speciation. The mainstream geography-of-speciation research has focused on gene flow patterns during speciation. However, this approach is inadequate without considering the distribution and range size to distinguish allopatric from non-allopatric speciation (Hernandez-Hernandez et al., 2021). Further, next-generation sequencing (NGS) techniques such as restriction site-associated DNA (RAD) improve the inferences of speciation processes (Maguilla et al., 2017; Nikolic et al., 2020; Goodall et al., 2021). Estimating allopatric-to-non-allopatric speciation frequencies is important to learn the community assembly process and evolutionary and ecological mechanisms of biodiversity (Skeels and Cardillo, 2019; García-Navas et al., 2020). However, considerations on the geographic mode of speciation are scarce in recent Cycas studies (Liu et al., 2015; Feng et al., 2016, 2021; Wang et al., 2019). Although most systematic studies conducted integrative taxonomy following the USC, few genetic markers were used, making them questionable given the porous and semipermeable features of the genome (Wu, 2001; Harrison, 2012; Wolf and Ellegren, 2017).
Cycas sect. Asiorientales comprises C. taitungensis C. F. Shen et al. (1994) and C. revoluta Thunb. The former is endemic to Taiwan with two small populations; the latter is broadly distributed in southern Kyushu, and the Ryukyu Archipelago, Japan, and a small population in Fujian, China, with large but fragmented populations (Liu et al., 2018; Figure 1). The Ryukyu Archipelago and Taiwan are at the continental margin, and the land configuration has changed drastically to become fragmented. They have been reconnected over time through crustal movements and paleoclimatic sea-level changes (Wang et al., 2014). The Kerama and Tokara Gaps are two major geological gaps dividing the Ryukyu Archipelago into the northern, central, and southern groups (Figure 1). For species with limited dispersibility like C. revoluta and C. taitungensis (Dehgan and Yuen, 1983; Kono and Tobe, 2007), these landscape dynamics may influence their population genetic structure and evolutionary trajectories. Previous studies of C. revoluta have shown that morphogenetic variations correspond to geography (Setoguchi et al., 2009; Kyoda and Setoguchi, 2010). Compared with C. revoluta, C. taitungensis exhibits high genetic diversity despite the small census population size (Chiang et al., 2009; Zheng et al., 2017; Chang et al., 2019). These two sister species are in the stage of incipient speciation, showing a species tree anomaly zone based on different markers (Chiang et al., 2009, 2018; Chang et al., 2019). The high ancestral polymorphism and historical gene flow further complicate the study of their speciation (Chiang et al., 2009, 2018; Smith et al., 2018; Chang et al., 2019). Although the southward founder speciation of C. taitungensis has been inferred (Chang et al., 2019), the few dominant markers, incomplete sampling, and ignorance of different gene flow possibilities would mask or even bias the inference, for example, miscalculating the heterozygotes or mistaking the intraspecific and interspecific variations, especially for those with large genomes (mean size = 26.2 pg/2C in Cycas) (Funk and Omland, 2003; Zonneveld, 2011).
Figure 1. Sampling localities of Cycas revoluta and C. taitungensis for double-digest restriction site-associated DNA sequencing. The gray route is the course of the Kuroshio Current. The red and brown colors denote northern and southern groups according to the morphogenetic structure. CMR: Coastal Mountain Range; N. Ryu, C. Ryu, and S. Ryu refer to the northern, central, and southern Ryukyus, respectively.
Taxonomically, C. taitungensis had long been mistaken for C. taiwaniana Carruth. in section Stangerioides (Shen et al., 1994). Due to the taxonomic misconception, the morphological comparison between C. taitungensis and C. revoluta is lacking or erroneous, even in the protolog of C. taitungensis (Shen et al., 1994). Shen et al. (1994) did not compare characters between the two species but just adopted the problematic comparison of Yamamoto (1928) (see taxonomic status below). Other references have distinguished C. revoluta from C. taitungensis by the following characters in the latter species: (1) shorter, narrower, and keeled leaves, and leaflets; (2) strongly recurved or revolute leaflet margins; (3) more tightly imbricate female cones; (4) smaller and ovate (vs. elliptical) seeds with darker sarcotesta and non-grooved sclerotesta; and (5) ovate to lanceolate megasporophyll terminal lamella (Wang, 1996; De Laubenfels and Adema, 1998; Chen and Stevenson, 1999; Whitelock, 2002; Hill, 2008). However, the compared specimens were mostly from restricted areas and did not take into account comprehensive intraspecific variation.
To understand the speciation process and clarify comprehensive genetic variations, we applied genome-wide double-digest RAD sequencing (ddRAD-seq) alongside broad-range sampling. Morphometrics with a focus on diagnostic traits was also evaluated. We aimed to (1) clarify the population genetic structure across the whole distribution range of sect. Asiorientales, (2) infer the geographic mode of speciation in C. revoluta and C. taitungensis, and (3) reappraise the species status in sect. Asiorientales. For the geographic mode of speciation, the allopatric and non-allopatric speciation hypotheses were tested. The former expected no gene flow during speciation, and the divergence time may correspond to the submergence of land bridges in the Ryukyu Archipelago-Taiwan. In contrast, gene flow can still be detected during species divergence if speciation was not driven by geographic isolation, and the divergence time may be dated in prevalent land bridges in the Ryukyu Archipelago-Taiwan. The gene flow and their divergence time will be estimated using the coalescence-based method, and the hypotheses will be tested via model selection under the Bayesian method.
Materials and Methods
The packages for analyses were all conducted in R version 3.6.1 (R Core Team, 2015).
Population Sampling for ddRAD-Seq
We sampled 26 populations of C. revoluta (n = 209) in the Ryukyu Archipelago and southern Kyushu in Japan, the Fujian Province in China, and two populations of C. taitungensis (n = 30) in Taiwan (Figure 1 and Supplementary Table 1). The numbers of sampling populations and individuals were different between the two species but inevitably biased due to the small census population size and merely two populations of C. taitungensis. Although the removal of siblings would inflate the genetic diversity and effective population size (Hendricks et al., 2018), the frequent vegetative bulbils fallen from C. revoluta may inversely underestimate these indices. Therefore, the individuals separated by at least 5 m were sampled to prevent the replication from asexual individuals. Also, we only selected the populations growing in clusters, and hence, the information from nearby siblings would be retained. Fresh leaflets were dried in silica gel and preserved at 4°C for DNA extraction.
Bioinformatics of Single Nucleotide Polymorphism Calling
Genomic DNA was extracted following a modified protocol from Doyle (1991). Genomic DNA quality and quantity were checked using a Qubit 2.0 fluorometer (Life Technologies, Carlsbad, CA, United States). The reduced-representation ddRAD approach was used to acquire reliable and high-coverage biallelic single nucleotide polymorphism (SNPs). The qualified DNA was double digested by Sbf1 (recognition site: 5′-CCTGCAGG-3′) and Msp1 (5′-GGCC-3′) and ligated to Illumina adapters with individual barcodes and library indices. Following the library construction protocol of Peterson et al. (2012) with size selection targeting 250–500 bp, 150-bp paired-end reads were sequenced in 12 libraries using the Illumina HiSeq X Ten (Illumina, San Diego, CA, United States).
The raw reads were checked by multiQC with a Phred score > 30 before de novo mapping. Stacks 2.2 (Catchen et al., 2013) was used to optimize the parameters in ustacks (M and n from 1 to 6) and cstacks (m from 3 to 6) considering the error rate, locus number, and coverage optimization following the methods of Mastretta-Yanes et al. (2015). The dataset was processed with a de novo pipeline with M = n = 6 and m = 3 without predefined groups. A total of 57 samples (N = 45 and 12 for C. revoluta and C. taitungensis, respectively) were used to construct catalogs in gstacks. Of the resulting SNPs, only one was retained in each RAD locus to decrease the linkage effect, and sites with a missing rate > 0.4 (p = 1, r = 0.4), a minor allele frequency (MAF) < 0.01 (–min–maf = 0.01), and maximum heterozygosity > 0.8 (–max–obs = 0.8) were excluded in the populations program. The PCR duplicates had little effect on genotype calling and heterozygote in high-quality DNA and sequencing data (Supplementary Table 2) and were therefore not removed (Euclide et al., 2020). Unlinked SNPs with a mean allele depth < 5, missing rate > 0.55, and individual missing rate > 0.5 were removed by Vcftools (Danecek et al., 2011).
To understand the population genetic structure of the whole distribution range, the model-free principal coordinate analysis (PCoA) and model-based ancestry coefficient estimations were performed. Nei’s distance matrix among sampling populations was transformed by nei.dist function in the poppr package (Kamvar et al., 2014) and input for PCoA in the cmdscale function from the stats package.
The ancestry coefficient was estimated without a predefined group to clarify the intra-section structure. There was a total of 14 populations in our study (Figure 2), and we assumed at most one population contained its specific structure. Hence, K was set from 1 to 15 in StrAuto 0.3.1 (Chhatre, 2012) and the sparse nonnegative matrix factorization (sNMF) algorithm in the snmf function of the LEA package (Frichot et al., 2015). The former analysis was optimized by 100,000 iterations with 10,000 burn-in, and the consensus was reached with 10 independent runs; the latter was optimized by 1,000 iterations, and the consensus was also reached with 10 independent runs. Evanno’s method in CLUMPAK (Kopelman et al., 2015) and cross-entropy were used to validate the best ancestral proportions (i.e., K) in StrAuto and sNMF, respectively.
Figure 2. Island population genetic structure inferred by (A) the sparse nonnegative matrix factorization (sNMF) algorithm and (B) the StrAuto software. CR and CT denote Cycas revoluta and Cycas taitungensis, respectively. The island names are abbreviated as the first three characters of the full names.
Population clustering exhibited larger intraspecific than interspecific variations, corresponding to the geographic distribution (see Results). Consequently, the boundaries and strengths of genetic barriers were estimated considering the sampling locality connections and Nei’s genetic distances in the adegenet 1.3–1 package (Jombart and Ahmed, 2011). Three connection networks were drawn by chooseCN function to represent the intensive [Delaunay triangulation (DT)], moderate [nearest neighbor (NN)], and low [minimum spanning tree (MP)] connection strength. Particularly for NN, due to the lack of dispersibility for sect. Asiorientales, K was set to 7 to consider the connections of a quarter of the sampling (i.e., N/4 = 28/4). The barriers were searched by 20 optimizations with the optimize.monmonier function and the local differences were set to zero to search possible boundaries. The island-level pairwise Fst values were also computed to quantify the divergence patterns using arlsumstat in the ABCtoolbox (Wegmann et al., 2010).
Speciation Modeling of Cycas revoluta and Cycas taitungensis
To test the hypotheses of the geographic mode of speciation, the approximate Bayesian computation (ABC) was conducted. Because of the different range sizes between the two cycads, four speciation models were proposed to consider gene flow and demography. Scenarios combining gene flow and demographic changes were set to fit the poor transoceanic dispersibility of these two species (Dehgan and Yuen, 1983; Kono and Tobe, 2007). Based on Coyne and Orr (2004), allopatric speciation is defined as complete geographic isolation without gene flow. Due to the lack of fossil records to track back the paleodistribution of the onset of speciation, allopatric, and non-allopatric speciation was determined by the absence and presence of gene flow in the incipient speciation stages, respectively (also see Dong et al., 2020). We also accounted for the post-divergence gene flow by secondary contact. Overall, the non-allopatric speciation models were (1) speciation with continuous gene flow (CG) and (2) primary contact (PC); the allopatric speciation models were (3) speciation with complete isolation (CI) and (4) secondary contact (SC). In parameter settings, the divergence of the two species was set < 9 mya based on Nagalingum et al. (2011) and transferred to 3 × 105 generations considering 30 years per generation. There are several records for the Cycas generation time, from more than 10 to 40 years in the International Union for Conservation of Nature (IUCN; Mankga and Yessoufou, 2017). Because the IUCN generation time records are rough and lack detailed justifications (Zheng et al., 2017), and C. taitungensis has been recorded to grow fast—5 cm annually (Huang et al., 2004)—we have conservatively taken 30 years for the Cycas generation time (Stewart and Globig, 2011; Table 1 and Figure 3; see Supplementary Table 3 and Supplementary Method 1 for prior settings and detailed model description).
Figure 3. The four speciation models testing the demography and gene flow patterns in the approximate Bayesian computation models. The parameters are defined in Table 1.
One million simulations for four divergence models were implemented in fastsimcoal 2 (Excoffier et al., 2013), and 20 summary statistics (i.e., species-specific segregating sites; number of pairwise differences; He and its standard deviation; global mean and standard deviation of segregating sites, number of pairwise differences, Fst, and Fit; interspecific Fst) were computed using arlsumstat in ABCtoolbox (Wegmann et al., 2010). Posterior probability was performed using the best 1,000 simulations after removing highly correlated summary statistics for the model selection by the Bayes factor (Jeffreys, 1998) in the ABCtoolbox (see Supplementary Method 1 for detailed quality control and validations before model selection). The precision of the best model on the observed data was quantified by the goodness-of-fit test. Parameters of the best model were log-transformed and estimated by retained summary statistics by 15 hidden layers and 50 feed-forward neural networks in non-linear regression using the abc function in abc package (Csilléry et al., 2012).
Population and Individual Phylogeny Reconstructions
Since frequent gene flow was detected between these two cycads (see “Results”), we re-delineated the species boundary using the phylogenetic analysis. The transformed Nei’s distance was used to reconstruct the neighbor-joining tree with 1,000 bootstraps by the aboot function in poppr package (Kamvar et al., 2014). Because of apparent genetic admixture, the individual tree was further reconstructed to detail the intraspecific relationship for species delimitation.
To reduce the intensive computation, we selected 110 individuals (6–9 individuals per island) with low individual missing rate (range from 16 to 50%) according to their genetic components. For C. taitungensis, the cryptic lineage (i.e., the green genetic component in STRUCTURE and sNMF) and the others were sampled separately (Supplementary Table 1). Further, we concatenated 907 RAD loci for phylogenetic reconstruction for a more accurate tree topology (Leache et al., 2015) via the –phylip-var-all option in Stacks 2.2 (Catchen et al., 2013). BEAST analysis (Bouckaert et al., 2014) was conducted with empirical site frequency and GTR+GAMMA in the substitution model and the gamma distribution birth rate in the Yule model prior with other settings default (xml file in figshare).1 We performed five independent runs, each conducting 50 million Markov chain Monte Carlo (MCMC) simulations and sampling every 1,000 simulations. The convergence and the posterior probability were then checked by Tracer 1.7 (Rambaut et al., 2018). The independent runs were all converged, and we selected the one with the highest posterior probability and effective sample size (ESS) > 200 as the input for TreeAnnotator (Bouckaert et al., 2014). The target tree was summarized by maximum clade credibility with 40% burn-in.
Tree-Based Species Delimitation by Discovery Methods
We used the individual tree to validate the species boundary using the discovery methods (Carstens et al., 2013) by the single-locus species delimitation, general mixed Yule-coalescent (GMYC) (Fujisawa and Barraclough, 2013), and Poisson tree processes (PTP) (Zhang et al., 2013), based on concatenated RAD loci (Fujisawa and Barraclough, 2013). We implemented the calibrated tree from BEAST2 in multiple threshold GMYC (mGMYC) and the Bayesian version of GMYC (bGMYC) in splits (Ezard et al., 2017) and bGMYC packages (Reid, 2014), respectively. In mGMYC, the target tree was estimated using the gmyc function. In bGMYC, 500 individual subtrees from BEAST were randomly sampled for 50,000 MCMC simulations, a 100 thinning interval, and 10% burn-in by the bgmyc.multiphylo function. The convergence was visualized by the spec.prombat function.
The PTP method considering only the number of substitutions was implemented in multi-rate PTP (mPTP) by standalone mptp (v 0.2.4) software (Kapli et al., 2017) and the Bayesian version of PTP (bPTP) in an online webserver.2. In mPTP, the minimum branch length threshold was detected by –minbr_auto and four independent runs with 10 million MCMC samplings; a 1,000 sampling interval and 10% burn-in were set. The combined likelihood plot and the three starting delimitations (i.e., randomly generated, maximum likelihood, and single exponential delimitation) were used to assess the convergence. For bPTP performed using the online webserver, the species were delimited by a target Bayesian individual tree with 500,000 generations with 100 sampling intervals and 20% burn-in. Both original maximum likelihood (PTP ML) and Bayesian implementations (bPTP) were reported.
In addition to genetic evidence, the morphological variations were also scrutinized to evaluate species status. Eleven vegetative and 10 reproductive diagnostic characters (Wang, 1996; De Laubenfels and Adema, 1998; Chen and Stevenson, 1999; Whitelock, 2002; Walters and Osborne, 2004; Hill, 2008) were selected for Cycas classification (Figure 4 and Table 2). Although the female cone compactness is a diagnostic trait between C. revoluta and C. taitungensis, we did not consider it because of the difficulty standardizing it due to large variation through the developmental stage (Walters and Osborne, 2004).
Figure 4. Morphometrics of leaf, leaflet, seed, and megasporophyll. The number block is the sample size. The ratio is calculated as the length divided by the width. The black and red stars are not significant and significant, respectively, diagnostic traits between Cycas revoluta and Cycas taitungensis. The draft of the measured organ is modified from Li and Keng (1975). Insert. angle, insertion angle; Leaflet num. per side, leaflet number per side; Mega. lateral spine num., megasporophyll lateral spine number; Mega. length, megasporophyll lamina length; Mega. ratio, megasporophyll lamina ratio; Mega. width, megasporophyll lamina width. The part labels a–q are the plant trait corresponding to depiction and Table 2.
The traits were measured in 109 specimens from nine herbaria (74 C. revoluta and 35 C. taitungensis) and 21 C. taitungensis from the coastal mountain range (N = 16) and the National Chung Hsing University campus (N = 5, cultivated individuals from Hongye Nature Reserve) (Supplementary Table 1). Specimen from herbaria HAST, KAG, KYO, PE, TAI, TAIF, TI, TNS, and UPS or their websites were examined. Herbarium acronyms follow Index Herbariorum (Thiers, 2019, continuously updated). Specimens with unclear collection information were excluded.
To measure the recurvation of leaflet margin, an additional 188 samples of leaflets used for DNA extraction (178 C. revoluta and 10 C. taitungensis) were examined (Supplementary Table 1). The recurvation degree of the leaflet margin was categorized into four classes (Figure 5A). The middle part of dried mature leaflets was applied to maintain consistency in the measurements. The Mann–Whitney U-test was used for testing the significant differences of interspecific variations by wilcox.test function of the stats package. Finally, four significantly different interspecific traits relevant to leaflets (i.e., leaflet width, ratio, basal width, and spacing; see “Results”) were measured to construct the island dendrogram. The traits were scaled first and transformed by Euclidean distance using the dist function of the stats package. Among four clustering methods [unweighted pair group method with arithmetic mean (UPGMA), complete and single linkage, and Ward’s method], the best one was chosen by evaluating the agglomerative coefficient. Then the dendrogram was drawn by the agnes function in the cluster package (Maechler et al., 2012).
Figure 5. Comparisons of the revolute extent in the leaflet margin (A) between species and (B) among geological regions. The number above the bars are measured sample size. CR, Cycas revoluta; CT, Cycas taitungensis. The cross section of leaflet is modified from Li and Keng (1975).
Bioinformatics of ddRAD-Seq and Genetic Diversity
After demultiplexing and mapping 194,992,789 pair-end reads from 248 samples, the average contig size was 242.8 bp with a mean overlap size of 28.9 bp. Among 806,955 loci, the mean coverage was 59.8× (minimum = 28.3×). The population program kept 2,055 loci with a mean length of 238.21 bp and a total length of 491,971 bp. After quality control, there were 907 unlinked SNPs among 239 individuals (n = 209 and 30 for C. revoluta and C. taitungensis, respectively) with a missing rate of 35.4%. Among 110 selected individuals for individual phylogenetic reconstruction, there were 188,059 bp with 3,509 SNPs in 907 concatenated RAD loci. The genetic diversity showed fewer private alleles, but higher He and π in C. taitungensis (He = 0.238; π = 36.58; private allele = 108) than C. revoluta (He = 0.138; π = 27.35; private allele = 260). The species divergence was low (Fst = 0.077). The SNP data were deposited in the Figshare repository and are available at https://doi.org/10.6084/m9.figshare.16987954.v1.
Discordance of the Population Structure With Geographic Gaps and Larger Intraspecific Than Interspecific Genetic Divergence in Cycas Sect. Asiorientales
The island population structure revealed K = 4 and 7 as the best clustering in sNMF and StrAuto, respectively (Supplementary Figure 1). Both results separated two groups clearly, north of the Amami Islands (hereafter the northern group) and south of the Okinawa Islands (hereafter the southern group) in the middle of the central Ryukyus (Figure 2). The genetic structure did not correspond to the geographic gaps and the species boundary. The genetic compositions across the Tokara and Kerama Gaps were more homogeneous with red and orange genetic components (Figure 2). At the species level, the genetic compositions of C. taitungensis were similar to C. revoluta of the central and northern Ryukyus except for the green lineage. The inconsistency of the geographic boundary also appeared in the Fujian population with similar genetic compositions with the northern group, which was geographically separated from Fujian by the Okinawa Trough.
The neighbor-joining tree and PCoA consistently exhibited greater intraspecific divergence of the northern and southern groups of C. revoluta than the interspecific divergence (Figure 6). Although C. taitungensis formed a monophyletic group, both northern and southern groups of C. revoluta were more closely related to C. taitungensis than each other (Figure 6A). The population barriers and island pairwise Fst also indicated stronger intraspecific than interspecific variations (Figure 7 and Supplementary Table 4). In three connection networks, the major barrier was congruently assigned as the strongest barrier between the southern and northern groups (pairwise Fst: 0.32–0.61), rather than between C. taitungensis and C. revoluta in the Ryukyus (Fst: 0.056–0.273) and between Fujian and the southern group (Fst: 0.11–0.298). For the weaker interspecific barrier, most C. taitungensis individuals were genetically admixed with C. revoluta, except for six individuals harboring the nearly fixed green genetic component. Consequently, the interspecific barrier assigned by Nei’s distance was mainly ascribed to this distinct cryptic green lineage. Another minor barrier that appears in the Tokara Gap in the MS and DT connection networks was concordant with the different genetic compositions of the Amami Islands in the northern group in sNMF when K = 7 (Figure 2).
Figure 6. Population genetic structure inferred by (A) the neighbor-joining tree and (B) principal coordinate analysis (PCoA) by Nei’s distance. The island names are abbreviated as the first three characters of the full names, except for CMR (Coastal Mountain Range). Within-island populations are listed in Supplementary Table 1.
Figure 7. Population genetic boundaries based on three connection networks. The solid black and orange lines represent the connection network and barriers, respectively, and the dashed black lines are the geological boundaries separating the northern, central, and southern Ryukyus.
Ongoing but Incomplete Speciation After the Divergence
Among the 100,000 simulations, 13 summary statistics were retained for the further model selections and parameter estimation after correlation removal. The confusion matrix in model validation discerned the four models appropriately, and the multidimensional space from summary statistics embraced the observed value in the PCA (Supplementary Figure 2). After model selection using 1 million simulations, the continuous gene flow was the best model with the highest Bayes factor and posterior probability (Table 3). Moreover, the observed data fit well with the simulated ones in the goodness-of-fit test (P = 0.51) (Supplementary Figure 3).
Regarding parameter estimation of the continuous gene flow model (Figure 3, Supplementary Figure 4, and Table 4), sect. Asiorientales diverged from its sister group t3 generations ago [weighted median: 162,546.4; weighted 95% highest posterior density (HPD): 26,558–282,427] with a large effective population size (weighted median: 37,409.2), and the two species diverged t2 generations ago (weighted median: 4,588.2; weighted 95% HPD: 51–86,158). At t1 generations ago (weighted median: 109.7; weighted 95% HPD: 15.5–4,176), bottleneck events occurred in both species with NCR1 (weighted median: 34,430.6) contracted 83% to NCR2 (weighted median: 5,835.7; weighted 95% HPD: 2,002.8–8,471.5), which is similar to a 77% contraction in NCT1 (weighted median: 12,067.5) to NCT2 (weighted median: 2,669.8; weighted 95% HPD: 614.2–4,697.1) that is observed today. Notably, the current larger effective population size of C. revoluta than C. taitungensis contradicts the findings of Chang et al. (2019). Following bottleneck events, the interspecific migration rates also decreased about an order of magnitude from mCT1 (weighted median: 3.1 × 10–4; weighted 95% HPD: 1.6 × 10–7–0.44) to mCT2 (weighted median: 1.8 × 10–5; weighted 95% HPD: 6.8 × 10–8–0.041) and mCR1 (weighted median: 1.4 × 10–5; weighted 95% HPD: 5.7 × 10–8–0.02) to mCR2 (weighted median: 4 × 10–6; weighted 95% HPD: 7 × 10–8–2.5 × 10–4). This reduced interspecific gene flow also indicates ongoing speciation. Nonetheless, the migration from C. taitungensis to C. revoluta is stronger than in the opposite direction.
Northward Serial Divergence of Cycas Sect. Asiorientales
The individual tree revealed polyphyly in C. taitungensis and C. revoluta and could be separated into five clades (A–E) with high posterior support (Figure 8). In terms of geography, there was a propensity of serial divergence in sect. Asiorientales from south to north of the Ryukyu Archipelago. Clade E representing the cryptic green lineage in C. taitungensis diverged first from other lineages. The subtending Clades C and D include all southern Ryukyus, plus the Iheya, Tonaki, and southern Okinawa Islands of the central Ryukyus. Specifically, Clade C was admixed with C. taitungensis and Fujian individuals, so Clades A and B diverged last. Considering Clades A and B, the former included Tonaki Island and the northern Okinawa Islands of the central Ryukyus, and the latter included the northern group.
Figure 8. The Bayesian individual tree reconstructed by concatenated restriction site-associated DNA (RAD) loci and the species delimitation results. The branch width indicates the posterior probability. The color bar in each region corresponds to geological regions (i.e., N., C., and S. Ryukyu; Taiwan; and Fujian). For the species delimitation results, the gray bars are posterior probabilities < 0.6. The island names are abbreviated as the first three characters of the full names. OkiN and OkiS represent the northern and southern Okinawa populations, respectively. The bar is unscaled time with a substitution rate equal to 1.
Inconsistency of Species Delimitation
The species delimitation failed between the two cycads. If C. revoluta and C. taitungensis are two “good” species, we expected a distinguished species boundary. However, the paraphyly of sect. Asiorientales rejected the delimitation of C. revoluta and C. taitungensis based on tree-based delimitation approaches (i.e., GMYC and PTP analogs). Results of GMYC and PTP reflect the intraspecies level (Figure 8). For GMYC, the multi-threshold and Bayesian versions obtained similar patterns. In mGMYC, 14 groups were recognized; Clades A and B clustered together, while four, seven, and two groups were recognized in Clades C, D, and E, respectively. bGMYC detected only nine groups that recovered Clades A, B, and E, while separating two and four groups in Clades C and D, respectively.
For PTP, the results were incongruent with GMYC and between multi-rate (mPTP) and non-multi-rate PTP versions (PTP ML and bPTP) (Figure 8). The Bayesian and likelihood versions were oversplit, with 81 and 60 groups discovered, respectively. The nearly individual-level separation indicated a strong false positive of delimitation (Luo et al., 2018). The mPTP outperformed the former version by accounting for intraspecific variation with different exponential distributions. As expected by Kapli et al. (2017), the mPTP yielded fewer putative species than PTP with Clade E recovered and merging of Clades A–D.
Overlapped Morphological Variation Between Cycas taitungensis and Cycas revoluta
Among 11 vegetative and 10 reproductive traits, there were nine traits (four vegetative and five reproductive) significantly different between the species (Table 2), and six of them were the described diagnostic traits (red stars in Figure 4 and leaflet recurvation in Figure 5). In general, C. taitungensis yielded wider and longer leaflets with larger inter-leaflet spacing and a wider basal width than C. revoluta, and the megasporophyll tended to be ovate and shorter with larger and longer seeds. However, these variations were continuous with C. taitungensis included in C. revoluta.
For the five non-significant reproductive traits, counterexamples were found between the two species for three diagnostic traits, namely, seed shape, sarcotesta color, and sclerotesta groove Supplementary Figures 5–7). The indistinguishable sclerotesta groove affirmed the descriptions from Hill (2008). For the seven non-significant vegetative traits, four diagnostic traits (leaflet recurvation, leaflet length, leaf insertion angle, and leaf length) were not different between the two species (black stars in Figure 4). Even though the epithet of C. revoluta is derived from the recurved leaflet margin, 15% of leaflet margins were slightly or not recurved (class 1) and, hence, were not discernible from C. taitungensis. Most flattened leaflets belonged to the northern group (Figure 5). Alternatively, the typical recurved leaflet margins were predominantly observed in the Fujian and southern groups. Focusing on the Ryukyus and Kyushu, these patterns are consistent with our and previous population genetic clustering (i.e., the northern and southern groups) (Kyoda and Setoguchi, 2010), pairwise Fst, and leaflet morphometric cladogram (Figures 2, 6, Supplementary Figure 8, and Supplementary Table 4; see also Setoguchi et al. (2009)).
Among four clustering methods, Ward’s method explained the most structure (agglomerative coefficient 0.94). The compound traits integrating the four significantly different leaflet traits failed to distinguish the species and showed the disjunct distribution of phenotypically similar populations. The dendrogram revealed a similar topology to the Bayesian individual tree, in which C. taitungensis was primarily clustered with the Fujian, Okinawa, and northern groups (Supplementary Figure 8). In summary, half of the diagnostic characters (5/10) were significantly different between species, but the variation of all interspecific traits overlapped.
Cycads are conserved in morphology and genetics. Comprehensive morpho-genetic data help validate taxonomic status by tracking their speciation history. This study combined the genome-wide ddRAD and morphological traits to illuminate the evolutionary history and reappraise species status of sect. Asiorientales. The ABC results indicated non-allopatric divergence following continuous gene flow. The morpho-genetic evidence unveiled larger intra- than inter-species variations laid the foundations of species reappraisal in sect. Asiorientales.
Continuous but Decreasing Gene Flow Supports Recent Non-allopatric Speciation
Our best ABC model suggests a continuous gene flow with the bottleneck after the divergence between C. revoluta and C. taitungensis. This scenario indicates non-allopatric speciation. The origin of sect. Asiorientales around 4.8 mya, as calculated as 30 × t3 (162,546 generations, Figure 3 and Table 4), matches the timing of the emergence of the Ryukyu Arc as a series of continental islands ca. 5–10 mya (Wang et al., 2014; Government of Japan, 2019), underpinning the vicariance of sect. Asiorientales from its sister sect. Panzhihuaenses on the continent (Xiao and Moller, 2015; Mankga et al., 2020). Incipient speciation occurred in the middle Pleistocene (t2, ca. 0.14 mya) after the synchronous isolation of the Ryukyu Arc into the northern, central, and southern Ryukyus, and Taiwan around 1.55 mya (Osozawa et al., 2013). Since then, isolation between the southern Ryukyus and Taiwan may have persisted due to the Yonaguni Strait (Gallagher et al., 2015; Zheng et al., 2016; Government of Japan, 2019). Nevertheless, the gene flow was continuous, thus rejecting allopatric speciation. The mismatch between the speciation event and the simultaneous landscape fragmentation suggests that continuous gene flow delayed the speciation, prolonged the lineage sorting time, and preserved high ancestral polymorphisms. Consequently, current spatial morphogenetic variations are shaped by both gene flow and ancestral polymorphism.
Mechanisms of Inter-Island Gene Flow After Divergence
After speciation, inter-island gene flow was facilitated by sea-level change. The land bridges prevailed in the Last Glacial Maximum (LGM) among the Ryukyu Arc, facilitating the colonization of species with poor transoceanic dispersal ability (Government of Japan, 2019). Transoceanic dispersal of cycad seeds and pollens has been documented. Cycad’s seed dispersibility varies from less than 1 meter to transcontinental by both abiotic (e.g., gravity, sea current, and rain) and biotic media (Burbidge and Whelan, 1982; Hall and Walter, 2013; Cousins and Witkowski, 2017; Liu et al., 2021a). However, some species could not be transported by sea due to a lack of spongy layer in seed, such as sect. Asiorientales (Dehgan and Yuen, 1983). Therefore, the stepping stone process was supposed in inter-island colonization of sect. Asiorientales via land bridge in LGM. The northward colonization may be dominant based on the calibrated individual tree (Figure 8). The genetic compositions remained homogeneous across the long-lasting Kerama and Tokara geological barriers than the Amami-Okinawa genetic boundary (Ota, 1998), indicating that further ecological factors might engage in the divergence of northern and southern groups.
In addition to stepping stone colonization, overseas gene flow may occur between the two more distantly distributed species after the bottleneck event. The overseas gene flow would be mainly from pollen carried by wind and insects. However, the wind-pollinated distance is short in cycads, and hence insects would be more likely as overseas pollination medium (Kono and Tobe, 2007). The primary pollinator, Carpophilus chalybeus, is broadly distributed in Taiwan and the Ryukyu Archipelago (Kirejtshuk, 2005; Chung and Shao, 2021). With the overlapped coning period between C. revoluta and C. taitungensis in spring to summer, insect pollination would facilitate bidirectional gene flow when the northeast and southwest monsoon are prevalent. Northward moving typhoons may also cause more frequent northward gene flow.
According to the individual tree, the gene flow would be spatially heterogeneous, suggesting other potential mechanisms (Figure 8). Apart from Clade E, all C. taitungensis individuals were closely related to Clades A–C of C. revoluta, which may be caused by more frequent gene flow than between C. taitungensis and Clade D of C. revoluta (Fontaine et al., 2015). Although the transoceanic seed dispersal of sect. Asiorientales is limited (Dehgan and Yuen, 1983; Liu et al., 2021a), other diaspores, such as bulbils, may lead to spatially heterogeneous gene flow. Many land plant transoceanic dispersal cases only focused on seeds, but few on vegetative dispersals (Ridley, 1930; Vargas et al., 2014, 2015; Price and Wagner, 2018). Bulbils are the vectors for fast and genetically homogeneous population colonization (Norstog and Nicholls, 1997; Thomas et al., 2005; Mizuki et al., 2010; Huang and Hsieh, 2014). Successful and profuse vegetative propagation has been reported in many cycad species (Raju and Jonathan, 2010; Raju and Rao, 2011; Cousins and Witkowski, 2017), but the importance of its transoceanic dispersal is less discussed (Ridley, 1930; Keppel, 2001). Bulbil development is commonly observed in near-sea streamside C. taitungensis and seashore C. revoluta (Stopes, 1910; Stevenson, 2020) and would lead to long-distance dispersal by hitchhiking by the Kuroshio Current, although not formally recorded. The Kuroshio Current is strong and warm, persistently flowing northwestward from the eastern Philippines, offshore eastern Taiwan, and along the western Ryukyu Archipelago. It changes its course eastward to the Pacific Ocean at the Tokara Gap (Gallagher et al., 2015; Figure 1). This route could be dated back to ca. 2–1.5 mya after the speciation of C. revoluta and C. taitungensis. The likelihood of the Kuroshio promoting dispersal and landing on islands is low but possible, especially under monsoon support (Kurita and Hikida, 2014; Yang et al., 2018; Kaifu et al., 2020). Over a period of 30 years, Kaifu et al. (2020) recorded that ca. 7% of drifters (8/114) analogous to bulbils reached the northern Ryukyus with higher probability than the central or southern Ryukyus by the Kuroshio, set off from Taiwan within half a month (for the coral species example, see Selmoni et al. (2020)), faster than the adventitious root development (ca. 1 year) in cycad bulbils (Dehgan, 1997; Stevenson, 2020).
The occasional long-distance vegetative dispersal following its faster growth and colonization (compared with seedlings) facilitates a founder-takes-all situation (i.e., new colonizers on a new habitat expand rapidly and block the latecomers; Waters et al., 2013). This spatially heterogeneous gene flow and genetic influence may explain the disjunct distribution and barrier patterns of morphogenetically similar populations of C. taitungensis in Taiwan and C. revoluta in the central and northern Ryukyus (Figure 7). The population genetic structure, pairwise FST, and paraphyletic individual tree demonstrated the cohesiveness of C. taitungensis and the central and northern Ryukyus of C. revoluta and not the southern Ryukyus (Figures 2, 6, 8 and Supplementary Table 4), similarly to the morphology, which exhibits less recurvation, a larger length-to-width ratio, spacing, and a wider leaflet and basal width (Figures 4, 5). Specifically, the typical form of strongly recurved leaflets in C. revoluta was discovered mainly in the southern group. These similarities of disjunct morphogenetic variations correspond to leapfrog patterns, which underlie the disjunct distribution of phenotypically similar populations (Remsen, 1984).
Causes of the Leapfrog Morphogenetic Pattern
There are two main hypotheses for the leapfrog pattern (Remsen, 1984; Cadena et al., 2011; Marquez et al., 2020). The first hypothesis is the shared common ancestor of disjunct populations, which could be attributed to the long-distance Kuroshio hitchhike of bulbils. The spatially heterogeneous gene flow decelerates the lineage sorting and retains a more common ancestral polymorphism between C. taitungensis and the northern group of C. revoluta. The second hypothesis is the evolutionary convergence of leapfrogged populations or adaptive divergence of the intervening populations, which this study cannot validate. The convergence would occur between C. taitungensis and the northern group of C. revoluta. Multiple selective sweeps may contribute to adaptive divergence of the southern group of the C. revoluta. Along the Ryukyu Archipelago, the leapfrog pattern has been discovered in other species (Matsumura et al., 2009), but has not been tested to be driven by natural selection or stochastic process. Consequently, future studies can assess possible environmental and geographic forces influencing the evolution of functional genes to test these hypotheses.
The Species-Definition Anomaly Zone Leads to Problematic Taxonomy
The phenomenon of the greater intraspecific divergence than between species is called the species-definition anomaly zone. It can occur when two species harbor extremely different geographical distributions under gene flow. The subdivided population will inflate the effective population size for the species with a larger range and facilitate the species-definition anomaly zone (Wang and Caballero, 1999; Jiao and Yang, 2021). The effective population size of C. revoluta is two times that of C. taitungensis, in concert with the spatially heterogeneous gene flow, possibly enhancing the impact of the species-definition anomaly zone. The larger intra- than inter-specific variations indicate the problematic taxonomy for C. taitungensis and C. revoluta. The morphogenetic leapfrog pattern suggests that the divergence of C. revoluta and C. taitungensis should be regarded at the population level rather than the species level, considering them as separately evolving metapopulations (based on the unified species concept) (Wiley, 1978; De Queiroz, 2007).
Reduced gene flow and allopatric distribution may have initiated the speciation of these two cycads. However, recent human-induced sympatry by the horticultural introduction of C. revoluta into Taiwan may disrupt C. taitungensis speciation (Chen and Stevenson, 1999) by messing with the gene pool of natural populations (Chiang et al., 2013). The compromise between natural and anthropogenic processes is worth considering in our case of C. revoluta and C. taitungensis.
Taxonomy Status of Cycas taitungensis: A Synonym of Cycas revoluta
Incomprehensive Morphological Comparisons Misclassify the Interspecific Variation
The leapfrog morphogenetic variations and compatible reproduction under the same coning time indicate interconnected evolutionary fates. Our previous model-based analyses also support high ancestral polymorphism with recent gene flow (Chiang et al., 2009, 2018). Genetic studies have consistently demonstrated the cohesion of the two species, but morphological studies are more complicated, especially regarding species diagnosis and morphological comparisons.
The misidentification of the two species is not only because C. taitungensis had long been misidentified as C. taiwaniana, but also because of the problematic morphological comparisons in the C. taitungensis protolog. The earliest collection of C. taitungensis was made by Shun-ichi Sasaki in 1920. His collection was identified as C. taiwaniana native to the continent (southern China) and classified in sect. Stangerioides until Shen et al. (1994) published C. taitungensis as a new species. However, its morphological comparisons are questionable due to the cited reference (Yamamoto, 1928), which erroneously compared C. revoluta with mixed samples of genuine C. taiwaniana and C. taitungensis. Yamamoto (1928) described that C. taiwaniana has longer leaflets, leaves, and petioles, a smaller revolute margin, a larger insertion angle, and more leaflets. For the reproductive organs, the shape of the male cone is subcylindrical, and the megasporophyll is orbicular to ovate for C. taiwaniana with glabrous ovule (oblong or ovate to lanceolate with tomentose ovule for C. revoluta). However, tomentose ovule is a synapomorphy of sect. Asiorientales (Walters and Osborne, 2004), and Yamamoto’s description suggest that the cited specimens had been misidentified. This also implies that the taxonomic studies (Kanehira, 1936; Li and Keng, 1954, 1975; Liu, 1960) following Yamamoto’s work are questionable (Shen et al., 1994).
After the publication of C. taitungensis (Shen et al., 1994), some questionable examined specimens were still from unknown or erroneous distributions (Chen and Stevenson, 1999; Whitelock, 2002). In addition, due to large morphological variations, insufficient specimen measurements also lead to biased diagnoses. Wang (1996) and Hill (2008) first compared specimens correctly, but their comparison may be biased regionally. They distinguished C. taitungensis from C. revoluta based on (1) longer and flatter non-keeled leaf and non-revolute leaflet, (2) darker and larger seeds, and (3) more tightly, cabbage-like female cones. However, Wang (1996) and Hill (2008) only examined C. revoluta from Fujian and one individual from the northern Ryukyus, respectively (Zheng et al., 2017). The bias also appeared in the seed diagnosis. Whitelock (2002) and Hill (2008) illustrated grooved or non-grooved sclerotesta of C. revoluta and C. taitungensis, contradicting other references that described C. revoluta and C. taitungensis as having a smooth and grooved sclerotesta, respectively (De Laubenfels and Adema, 1998; Chen and Stevenson, 1999). In summary, the aforementioned diagnostic characters were insufficient and did not consider all variations. Our study for a more comprehensive comparison revealed continuous variations in these traits.
Why Not Define Subspecies?
In many species, an allopatric or edge population like C. taitungensis of sect. Asiorientales in Taiwan with slight morphogenetic differences would be treated as subspecies. However, the subspecies treatment is not recommended in cycads (Walters and Osborne, 2004) and would also not be appropriate in our case. Subspecies are traditionally geographic variants of polymorphic species that may be in the early stage of speciation (Mayr, 1982). The applications, however, are not standardized, and there are pronounced differences across taxonomic groups (Table 1; Huang and Knowles, 2016). The incongruence derives from anchoring the discrete subspecies through the initial stage of a speciation continuum.
Statistically evaluating and comparing the divergence level from multidimensional axes between species and subspecies would be a more objective way of assigning ranks (Huang and Knowles, 2016). Species units may still be the greatest challenge for cycad systematics due to the lack of intraspecific data (Walters and Osborne, 2004), let alone comparing species and subspecies levels across the speciation continuum. Therefore, the level below species is not recommended. Even taking the broad view of subspecies (de Queiroz, 2020), equating it at the species level but within a more inclusive lineage, the larger intraspecific divergence between the southern and northern groups of C. revoluta than the interspecific divergence between C. revoluta and C. taitungensis also markedly reduces inclusiveness. Alternatively, if the three subspecies were assigned (i.e., Taiwan and the northern and southern groups of the Ryukyu Archipelago), the taxonomic units are still confused and oversplit due to several morphological variants on small islands (Satake, 1981; Whitelock, 2002) and based on our species delimitation results (Figure 8). In summary, synonymization would be the most appropriate taxonomic treatment.
According to the results and discussion, we established the following treatment.
Cycas revoluta Thunb. in Verhandelingen Uitgegeeven Door De Hollandse Maatschappy Der Weetenschappen, Te Harlem 20 (2): 424, 426–427 (1782)
= C. taitungensis C.F. Shen, K.D. Hill, C.H. Tsou and C.J. Chen in Botanical Bulletin of Academia Sinica 35:133–140 (1994). Type: Taiwan. Taitung county, growing on the steep mountain slope in the secondary forest or well-exposed areas. C. H. Tsou 825 (Holotype: HAST!, Isotype: A, BM, IBSC, K, NSW, NY, P, PE!, TAI!) syn. nov.
Cycas revoluta is widely distributed along southern Kagoshima, the Ryukyu Archipelago, and eastern Taiwan. The distribution of natural populations in Fujian needs further verification (Zheng et al., 2017). In Japan, the habitats are mostly near the coastal area on the coral reef. Some populations are under or near forests with an altitude from sea level to 500 m. In Taiwan, two populations are located in the Taitung Hongye Village Taitung Cycas Nature Reserve (台東紅葉村台東蘇鐵自然保留區) along the Luyeh River Valley and the Coastal Range Taitung Cycas Forest Reserve (海岸山脈台東蘇鐵自然保護區), the western and eastern side of the East Rift Valley, respectively. Both populations generally grow on steep and exposed slopes with some in sparse forests at 300–900 m.
The threats and drivers of global biodiversity decline are mostly anthropogenically mediated (Mankga and Yessoufou, 2017), as is sect. Asiorientales. Extensive habitat disturbance and individual poaching have led to population fragmentation (e.g., C. taitungensis in Taiwan) and even possibly local extinction (e.g., C. revoluta in the Fujian population) (Shen et al., 1994; Chen and Stevenson, 1999; Whitelock, 2002; Huang et al., 2004). Secondary succession in lowland forest and pests, such as cycad aulacaspis scale (CAS) Aulacaspis yasumatsui and the cycad blue butterfly Chilades pandava, also reduce the habitat quality and individual health (Bailey and Lai, 2006; Wu et al., 2010). However, these threats are more severe for C. taitungensis than C. revoluta. The latest conservation assessment of the IUCN Red List revealed the improvement of C. revoluta from near threatened (NT) to least concern (LC) (Hill, 2010), while C. taitungensis has gone from vulnerable (VU) to endangered (EN) (Haynes, 2010). The worse conservation status of C. taitungensis is due to the increasing CAS and decreasing habitat quality and range size (Haynes, 2010). The genetic drift will enhance the subsequent reduction of heterozygosity if habitat loss continuous. With fewer private alleles than C. revoluta (108 vs. 260), the conservation status of the Taiwanese population will be inferior.
With the break of the species boundary, even if the Taiwan population is described as a unique evolutionary significant unit (ESU) (Chang et al., 2019), it will not be considered in the IUCN global assessment of species extinction risk. This situation signifies the negative effect of synonymized treatment on conservation evaluation. The absence of the conservation concern after synonym treatment has also appeared in other “globally common, but locally rare” taxa (Thakur et al., 2018). Under the considerations of evolutionary processes in conservation biology (Crandall et al., 2000), species as the only management unit is inappropriate (Casetta et al., 2019). Consequently, IUCN global assessment only at the species level would be insufficient.
Nevertheless, the guidelines for regional or national levels have also been proposed and could alleviate this demerit (IUCN, 2012). Such an index has been followed by the Taiwanese government to assess the conservation status of C. taitungensis and was classified as the “nationally critically endangered (NCR)” for its small area of occupancy (AOO) and fluctuation in the number of reproducible individuals (Yang et al., 2017). After our taxonomic treatment, the conservation of C. revoluta in Taiwan still needs to be considered on a national scale. This index could solve the predicament by considering the edge, higher evolvability, or regionally endangered populations and evaluating the conservation status more carefully. More attention on ESU at the regional or national scale would make the conservation management more comprehensive and focused.
The examined specimens from morphometrics are listed in Supplementary Table 5.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://figshare.com/, https://doi.org/10.6084/m9.figshare.16987990.v1, https://figshare.com/, https://doi.org/10.6084/m9.figshare.16987954.v1.
P-CL conceived, designed the project, and interpreted the data with J-TC. J-TC collected field samples, performed lab experiments and statistical analyses, and wrote the manuscript. C-TC, KN, M-XL, and H-LL critically reviewed and provided constructive suggestions on the first draft. All authors read and approved the final manuscript.
The ddRAD techniques were supported by Jui-Hua Chu from Technology Commons, College of Life Science at National Taiwan University. This project was supported by the Ministry of Science and Technology (MOST 110-2628-B-003-001 and MOST 109-2621-B-003-003-MY3) and subsidized by the National Taiwan Normal University (NTNU).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
We are grateful to Bing-Hong Huang, Nguyen Hung, Ariella Gladstein, and Huan-Yi Hsiung for their suggestions for data analyses. We also express our gratitude to Chien-Fan Chen and his colleagues of the Taiwan Forest Research Institute, and Masaharu Amano of the Okinawa Churashima Research Center for helping to photograph specimens.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.824158/full#supplementary-material
Bailey, R., and Lai, P. (2006). Bionomics of Cybocepahalus nipponicus endrody-younga (coleoptera: cybocephaidae) preying on cycad scale, Aulacaspis yasumatsui takagi (hemiptera: diaspididae). J. Int. Coop. 1, 1–14.
Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C.-H., Xie, D., et al. (2014). BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 10:e1003537. doi: 10.1371/journal.pcbi.1003537
Burbidge, A. H., and Whelan, R. J. (1982). Seed dispersal in a cycad, Macrozamia riedlei. Aust. J. Ecol. 7, 63–67. doi: 10.1111/j.1442-9993.1982.tb01300.x
Cadena, C. D., Cheviron, Z. A., and Funk, W. C. (2011). Testing the molecular and evolutionary causes of a ‘leapfrog’ pattern of geographical variation in coloration. J. Evol. Biol. 24, 402–414. doi: 10.1111/j.1420-9101.2010.02175.x
Carstens, B. C., Pelletier, T. A., Reid, N. M., and Satler, J. D. (2013). How to fail at species delimitation. Mol. Ecol. 22, 4369–4383. doi: 10.1111/mec.12413
Casetta, E., Silva, J. M. D., and Vecchi, D. (2019). From Assessing to Conserving Biodiversity: Conceptual and Practical Challenges. Cham: Springer.
Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A., and Cresko, W. A. (2013). Stacks: an analysis tool set for population genomics. Mol. Ecol. 22, 3124–3140. doi: 10.1111/mec.12354
Chang, J.-T., Huang, B.-H., and Liao, P.-C. (2019). Genetic evidence of the southward founder speciation of Cycas taitungensis from ancestral C. revoluta along the Ryukyu Archipelagos. Conserv. Genet. 20, 1045–1056. doi: 10.1007/s10592-019-01193-1
Chen, C.-J., and Stevenson, D. W. (1999). Flora of China, cycadaceae. Flora China 4, 1–7. doi: 10.1016/j.pld.2018.07.002
Chhatre, V. E. (2012). StrAuto: A Python Program.
Chiang, Y.-C., Huang, B.-H., Chang, C.-W., Wan, Y.-T., Lai, S.-J., Huang, S., et al. (2013). Asymmetric introgression in the horticultural living fossil cycas sect. Asiorientales using a genome-wide scanning approach. Int. J. Mol. Sci 14, 8228–8251. doi: 10.3390/ijms14048228
Chiang, Y.-C., Huang, B.-H., Li, N., Gong, Y.-Q., Shen, H.-H., Huang, S., et al. (2018). “Multilocus genome evidence for a paraphyletic relationship and past interspecific gene flow between species of Cycas section Asiorientales,” in Proceedings of the 9th International Conference on Cycad Biology: Cycad Biology and Conservation, Bronx, NY, 230–249.
Chiang, Y.-C., Hung, K.-H., Moore, S.-J., Ge, X.-J., Huang, S., Hsu, T.-W., et al. (2009). Paraphyly of organelle DNAs in Cycas Sect. Asiorientales due to ancient ancestral polymorphisms. BMC Evol. Biol. 9:161. doi: 10.1186/1471-2148-9-161
Chung, K.-F., and Shao, K. T. (2021). Catalogue of Life in Taiwan. Web Electronic Publication [Online]. Available online at: http://taibnet.sinica.edu.tw (accessed June 5, 2021).
Cousins, S., and Witkowski, E. (2017). African cycad ecology, ethnobotany and conservation: a synthesis. Bot. Rev. 83, 152–194. doi: 10.1007/s12229-017-9183-4
Coyne, J. A., and Orr, H. A. (2004). Speciation. Sunderland, MA: Sinauer Associates, Inc.
Crandall, K. A., Bininda-Emonds, O. R., Mace, G. M., and Wayne, R. K. (2000). Considering evolutionary processes in conservation biology. Trends Ecol. Evol. 15, 290–295. doi: 10.1016/s0169-5347(00)01876-0
Csilléry, K., François, O., and Blum, M. G. (2012). abc: an R package for approximate Bayesian computation (ABC). Methods Ecol. Evol. 3, 475–479. doi: 10.1111/j.2041-210x.2011.00179.x
Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., Depristo, M. A., et al. (2011). The variant call format and VCFtools. Bioinformatics 27, 2156–2158. doi: 10.1093/bioinformatics/btr330
De Laubenfels, D. J., and Adema, F. (1998). A taxonomic revision of the genera Cycas and Epicycas gen. nov.(Cycadaceae). Blumea 43, 351–400.
De Queiroz, K. (2007). Species concepts and species delimitation. Syst. Biol. 56, 879–886. doi: 10.1080/10635150701701083
de Queiroz, K. (2020). An updated concept of subspecies resolves a dispute about the taxonomy of incompletely separated lineages. Herpetol. Rev. 51, 459–461.
Degnan, J. H., and Rosenberg, N. A. (2006). Discordance of species trees with their most likely gene trees. PLoS Genet 2:e68. doi: 10.1371/journal.pgen.0020068
Dehgan, B. (1997). “Propagation and culture of cycads: a practical approach,” in Proceedings of the II International Symposium on Ornamental Palms & Other Monocots From the Tropics 486, Tenerife, 123–132. doi: 10.17660/actahortic.1999.486.16
Dehgan, B., and Yuen, C. K. K. H. (1983). Seed morphology in relation to dispersal, evolution, and propagation of Cycas L. Bot. Gazette 144, 412–418. doi: 10.1086/337391
Dong, F., Hung, C. M., and Yang, X. J. (2020). Secondary contact after allopatric divergence explains avian speciation and high species diversity in the Himalayan-Hengduan Mountains. Mol. Phylogenet. Evol. 143:106671. doi: 10.1016/j.ympev.2019.106671
Doyle, J. (1991). “DNA protocols for plants,” in Molecular Techniques in Taxonomy. NATO ASI Series (Series H: Cell Biology), Vol. 57, eds G. M. Hewitt, A. W. B. Johnston, and J. P. W. Young (Berlin: Springer).
Euclide, P. T., Mckinney, G. J., Bootsma, M., Tarsa, C., Meek, M. H., and Larson, W. A. (2020). Attack of the PCR clones: rates of clonality have little effect on RAD-seq genotype calls. Mol. Ecol. Resour. 20, 66–78. doi: 10.1111/1755-0998.13087
Excoffier, L., Dupanloup, I., Huerta-SañNchez, E., Sousa, V. C., and Foll, M. (2013). Robust demographic inference from genomic and SNP data. PLoS Genet. 9:e1003905. doi: 10.1371/journal.pgen.1003905
Ezard, T., Fujisawa, T., and Barraclough, T. (2017). SPecies’ LImits by Threshold Statistics. R package version [Online]. Available online at: https://R-Forge.R-project.org/projects/splits/ (accessed May 20, 2021).
Feng, X. Y., Wang, X. H., Chiang, Y. C., Jian, S. G., and Gong, X. (2021). Species delimitation with distinct methods based on molecular data to elucidate species boundaries in the Cycas taiwaniana complex (Cycadaceae). Taxon 70, 477–491. doi: 10.1002/tax.12457
Feng, X., Liu, J., and Gong, X. (2016). Species delimitation of the Cycas segmentifida complex (cycadaceae) resolved by phylogenetic and distance analyses of molecular data. Front. Plant Sci. 7:134. doi: 10.3389/fpls.2016.00134
Fontaine, M. C., Pease, J. B., Steele, A., Waterhouse, R. M., Neafsey, D. E., Sharakhov, I. V., et al. (2015). Extensive introgression in a malaria vector species complex revealed by phylogenomics. Science 347:125824. doi: 10.1126/science.1258524
Freudenstein, J. V., Broe, M. B., Folk, R. A., and Sinn, B. T. (2017). Biodiversity and the species concept-lineages are not enough. Syst. Biol. 66, 644–656. doi: 10.1093/sysbio/syw098
Frichot, E., François, O., and O’meara, B. (2015). LEA: an R package for landscape and ecological association studies. Methods Ecol. Evol. 6, 925–929. doi: 10.1111/2041-210x.12382
Fujisawa, T., and Barraclough, T. G. (2013). Delimiting species using single-locus data and the generalized mixed yule coalescent approach: a revised method and evaluation on simulated data sets. Syst. Biol. 62, 707–724. doi: 10.1093/sysbio/syt033
Funk, D. J., and Omland, K. E. (2003). Species-level paraphyly and polyphyly: frequency, causes, and consequences, with insights from animal mitochondrial DNA. Annu. Rev. Ecol. Evol. Syst. 34, 397–423. doi: 10.1016/j.ympev.2009.08.024
Gallagher, S. J., Kitamura, A., Iryu, Y., Itaki, T., Koizumi, I., and Hoiles, P. W. (2015). The Pliocene to recent history of the Kuroshio and Tsushima Currents: a multi-proxy approach. Prog. Earth Planet. Sci. 2:17.
García-Navas, V., Kear, B. P., and Westerman, M. (2020). The geography of speciation in dasyurid marsupials. J. Biogeogr. 47, 2042–2053. doi: 10.1111/jbi.13852
Goodall, J., Westfall, K. M., Magnúsdóttir, H., Pálsson, S., örnólfsdóttir, E. B., and Jónsson, Z. O. (2021). RAD sequencing of common whelk, Buccinum undatum, reveals fine-scale population structuring in Europe and cryptic speciation within the North Atlantic. Ecol. Evol. 11, 2616–2629. doi: 10.1002/ece3.7219
Government of Japan (2019). Nomination of Amami-Oshima Island, Tokunoshima Island, Northern Part of Okinawa Island, and Iriomote Island for Inscription on the World Heritage List. Government of Japan.
Habib, S., Dong, S., Liu, Y., Liao, W., and Zhang, S. (2021). The complete mitochondrial genome of Cycas debaoensis revealed unexpected static evolution in gymnosperm species. PLoS One 16:e0255091. doi: 10.1371/journal.pone.0255091
Hall, J. A., and Walter, G. H. (2013). Seed dispersal of the Australian cycad Macrozamia miquelii (Zamiaceae): are cycads megafauna-dispersed “grove forming” plants? Am. J. Bot. 100, 1127–1136. doi: 10.3732/ajb.1200115
Harrison, R. G. (2012). The language of speciation. Evolution 66, 3643–3657. doi: 10.1111/j.1558-5646.2012.01785.x
Haynes, J. (2010). Cycas taitungensis. IUCN Red List Threat. Species 2010:e.T42067A10642224.
Hendricks, S., Anderson, E. C., Antao, T., Bernatchez, L., Forester, B. R., Garner, B., et al. (2018). Recent advances in conservation and population genomics data analysis. Evol. Appl. 11, 1197–1211. doi: 10.1111/eva.12659
Hernandez-Hernandez, T., Miller, E. C., Roman-Palacios, C., and Wiens, J. J. (2021). Speciation across the Tree of Life. Biol. Rev. Camb. Philos. Soc. 96, 1205–1242.
Hill, K. D. (2008). The genus Cycas (Cycadaceae) in China. Telopea 12, 71–118. doi: 10.7751/telopea20085804
Hill, K. D. (2010). Cycas revoluta. IUCN Red List Threat. Species 2010:e.T42080A10622557.
Huang, C.-T., and Hsieh, C.-F. (2014). Asexual bulbil development and diversification of reproductive strategy between Remusatia vivipara and Remusatia pumila (Araceae). Taiwania 59, 220–230.
Huang, J.-P., and Knowles, L. L. (2016). The species versus subspecies conundrum: quantitative delimitation from integrating multiple data types within a single bayesian approach in Hercules beetles. Syst. Biol. 65, 685–699. doi: 10.1093/sysbio/syv119
Huang, S., Hsieh, H. T., Fang, K., and Chiang, Y. C. (2004). Patterns of genetic variation and demography of Cycas taitungensis in Taiwan. Bot. Rev. 70, 86–92. doi: 10.1663/0006-8101(2004)070[0086:pogvad]2.0.co;2
IUCN (2012). Guidelines for Application of IUCN Red List Criteria at Regional and National Levels: Version 4.0. Cambridge: IUCN.
Jeffreys, H. (1998). The Theory of Probability. Oxford: Oxford University Press.
Jiao, X., and Yang, Z. (2021). Defining species when there is gene flow. Syst. Biol. 70, 108–119. doi: 10.1093/sysbio/syaa052
Jombart, T., and Ahmed, I. (2011). adegenet 1.3-1: new tools for the analysis of genome-wide SNP data. Bioinformatics 27, 3070–3071. doi: 10.1093/bioinformatics/btr521
Kaifu, Y., Kuo, T. H., Kubota, Y., and Jan, S. (2020). Palaeolithic voyage for invisible islands beyond the horizon. Sci. Rep. 10:19785. doi: 10.1038/s41598-020-76831-7
Kamvar, Z. N., Tabima, J. F., and Grünwald, N. J. J. P. (2014). Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ 2:e281. doi: 10.7717/peerj.281
Kanehira, R. (1936). Formosan Trees. Taiwan: Kyushu Imperial University.
Kapli, P., Lutteropp, S., Zhang, J., Kobert, K., Pavlidis, P., Stamatakis, A., et al. (2017). Multi-rate Poisson tree processes for single-locus species delimitation under maximum likelihood and Markov chain Monte Carlo. Bioinformatics 33, 1630–1638. doi: 10.1093/bioinformatics/btx025
Keppel, G. (2001). Notes on the natural history of Cycas seemannii (Cycadaceae). South Pac. J. Nat. Appl. Sci. 19, 35–41. doi: 10.1071/sp01007
Kirejtshuk, A. G. (2005). On the fauna of Nitidulidae (Insecta, Coleoptera) from Taiwan with some taxonomical notes. Ann. Hist. Nat. Musei Natl. Hung. 97, 51–113.
Kono, M., and Tobe, H. (2007). Is Cycas revoluta (Cycadaceae) wind-or insect-pollinated? Am. J. Bot. 94, 847–855. doi: 10.3732/ajb.94.5.847
Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., and Mayrose, I. (2015). Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 15, 1179–1191. doi: 10.1111/1755-0998.12387
Kurita, K., and Hikida, T. (2014). Divergence and long-distance overseas dispersals of island populations of the Ryukyu five-lined skink, Plestiodon marginatus (Scincidae: Squamata), in the Ryukyu Archipelago, Japan, as revealed by mitochondrial DNA phylogeography. Zool. Sci. 31, 187–194. doi: 10.2108/zs130179
Kyoda, S., and Setoguchi, H. (2010). Phylogeography of Cycas revoluta Thunb. (Cycadaceae) on the Ryukyu Islands: very low genetic diversity and geographical structure. Plant Syst. Evol. 288, 177–189. doi: 10.1007/s00606-010-0322-1
Leache, A. D., Banbury, B. L., Felsenstein, J., De Oca, A. N., and Stamatakis, A. (2015). Short tree, long tree, right tree, wrong tree: new acquisition bias corrections for inferring SNP phylogenies. Syst. Biol. 64, 1032–1047. doi: 10.1093/sysbio/syv053
Li, H.-L., and Keng, H. (1954). Icones gymnospermum formosanarum. Taiwania 5, 25–83.
Li, H.-L., and Keng, H. (1975). Flora of Taiwan, Cycadaceae. Taipei: National Science Council of the Republic of China.
Liu, J., Lindstrom, A. J., Marler, T. E., and Gong, X. (2021b). Not that young: combining plastid phylogenomic, plate tectonic and fossil evidence indicates a Palaeogene diversification of Cycadaceae. Ann. Bot. 129, 217–230. doi: 10.1093/aob/mcab118
Liu, J., Lindstrom, A. J., Chen, Y. S., Nathan, R., and Gong, X. (2021a). Congruence between ocean-dispersal modelling and phylogeography explains recent evolutionary history of Cycas species with buoyant seeds. New Phytol. 232, 1863–1875. doi: 10.1111/nph.17663
Liu, J., Zhang, S., Nagalingum, N. S., Chiang, Y.-C., Lindstrom, A. J., and Gong, X. (2018). Phylogeny of the gymnosperm genus Cycas L. (Cycadaceae) as inferred from plastid and nuclear loci based on a large-scale sampling: evolutionary relationships and taxonomical implications. Mol. Phylogenet. Evol. 127, 87–97. doi: 10.1016/j.ympev.2018.05.019
Liu, J., Zhou, W., and Gong, X. (2015). Species delimitation, genetic diversity and population historical dynamics of Cycas diannanensis (Cycadaceae) occurring sympatrically in the Red River region of China. Front. Plant Sci. 6:696. doi: 10.3389/fpls.2015.00696
Liu, T.-S. (1960). Illustrations of Native and Introduced Ligneous Plants of Taiwan. Taikei: National Taiwan University, College of Agriculture.
Luo, A., Ling, C., Ho, S. Y. W., and Zhu, C. D. (2018). Comparison of methods for molecular species delimitation across a range of speciation scenarios. Syst. Biol. 67, 830–846. doi: 10.1093/sysbio/syy011
Maechler, M., Rousseeuw, P., Struyf, A., Hubert, M., and Hornik, K. (2012). Cluster: cluster analysis basics and extensions. R Package Version 1:56
Maguilla, E., Escudero, M., Hipp, A. L., and Luceno, M. (2017). Allopatric speciation despite historical gene flow: divergence and hybridization in Carex furva and C. lucennoiberica (Cyperaceae) inferred from plastid and nuclear RAD-seq data. Mol. Ecol. 26, 5646–5662. doi: 10.1111/mec.14253
Mankga, L. T., and Yessoufou, K. (2017). Factors driving the global decline of cycad diversity. AoB Plants 9:lx022. doi: 10.1093/aobpla/plx022
Mankga, L. T., Yessoufou, K., Mugwena, T., and Chitakira, M. (2020). The cycad genus Cycas may have diversified from indochina and occupied its current ranges through vicariance and dispersal events. Front. Ecol. Evol. 8:44. doi: 10.3389/fevo.2020.00044
Marquez, R., Linderoth, T. P., Mejia-Vargas, D., Nielsen, R., Amezquita, A., and Kronforst, M. R. (2020). Divergence, gene flow, and the origin of leapfrog geographic distributions: the history of colour pattern variation in Phyllobates poison-dart frogs. Mol. Ecol. 29, 3702–3719. doi: 10.1111/mec.15598
Mastretta-Yanes, A., Arrigo, N., Alvarez, N., Jorgensen, T. H., Pinero, D., and Emerson, B. C. (2015). Restriction site-associated DNA sequencing, genotyping error estimation and de novo assembly optimization for population genetic inference. Mol. Ecol. Resour. 15, 28–41. doi: 10.1111/1755-0998.12291
Matsumura, S. I., Yokoyama, J., Fukuda, T., and Maki, M. (2009). Origin of the disjunct distribution of flower colour polymorphism within Limonium wrightii (Plumbaginaceae) in the Ryukyu Archipelago. Biol. J. Linn. Soc. Lond. 97, 709–717. doi: 10.1111/j.1095-8312.2009.01253.x
Mayr, E. (1982). Of what use are subspecies? Auk 99, 593–595.
Mizuki, I., Ishida, K., Tani, N., and Tsumura, Y. (2010). Fine-scale spatial structure of genets and sexes in the dioecious plant Dioscorea japonica, which disperses by both bulbils and seeds. Evol. Ecol. 24, 1399–1415. doi: 10.1007/s10682-010-9396-z
Nagalingum, N., Marshall, C., Quental, T., Rai, H., Little, D., and Mathews, S. (2011). Recent synchronous radiation of a living fossil. Science 334, 796–799. doi: 10.1126/science.1209926
Nikolic, N., Liu, S., Jacobsen, M. W., Jonsson, B., Bernatchez, L., Gagnaire, P. A., et al. (2020). Speciation history of European (Anguilla anguilla) and American eel (A. rostrata), analysed using genomic data. Mol. Ecol. 29, 565–577. doi: 10.1111/mec.15342
Norstog, K., and Nicholls, T. J. (1997). The Biology of the Cycads. Ithaca, NY: Cornell University Press.
Osozawa, S., Su, Z.-H., Oba, Y., Yagi, T., Watanabe, Y., and Wakabayashi, J. (2013). Vicariant speciation due to 1.55?Ma isolation of the Ryukyu islands, Japan, based on geological and GenBank data. Entomol. Sci. 16, 267–277. doi: 10.1111/ens.12037
Ota, H. (1998). Geographic patterns of endemism and speciation in amphibians and reptiles of the Ryukyu Archipelago, Japan, with special reference to their paleogeographical implications. Res. Popul. Ecol. 40, 189–204. doi: 10.1007/bf02763404
Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S., and Hoekstra, H. E. (2012). Double digest RADseq: an inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS One 7:e37135. doi: 10.1371/journal.pone.0037135
Price, J. P., and Wagner, W. L. (2018). Origins of the Hawaiian flora: phylogenies and biogeography reveal patterns of long-distance dispersal. J. Syst. Evol. 56, 600–620. doi: 10.1111/jse.12465
R Core Team (2015). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.
Raju, A. S., and Jonathan, K. H. (2010). Anemophily, accidental cantharophily, seed dispersal and seedling ecology of Cycas sphaerica Roxb.(Cycadaceae), a data-deficient red-listed species of northern Eastern Ghats. Curr. Sci. 99, 1105–1111.
Raju, A., and Rao, N. G. (2011). Taxonomic aspects and coning ecology of Cycas circinalis L.(Cycadales: Cycadaceae), a threatened species of India. J. Threat. Taxa 3, 1425–1431. doi: 10.11609/jott.o2372.1425-31
Rambaut, A., Drummond, A. J., Xie, D., Baele, G., and Suchard, M. A. (2018). Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67:901. doi: 10.1093/sysbio/syy032
Reid, N. (2014). bGMYC: a Bayesian MCMC implementation of the general mixed Yule-coalescent model for species delimitation. R Package Version 1. Available online at: https://nreid.github.io/assets/bGMYC_instructions_14.03.12.txt. doi: 10.1186/1471-2148-12-196 (accessed May, 2021).
Remsen, J. (1984). High incidence of “leapfrog” pattern of geographic variation in Andean birds: implications for the speciation process. Science 224, 171–173. doi: 10.1126/science.224.4645.171
Ridley, H. N. (1930). The Dispersal of Plants Throughout the World. Ashford: L. Reeve & Company, Limited.
Satake, T. (1981). The varieties of Cycas revoluta. Cycad Newsl. 4, 3–10.
Selmoni, O., Rochat, E., Lecellier, G., Berteaux-Lecellier, V., and Joost, S. (2020). Seascape genomics as a new tool to empower coral reef conservation strategies: an example on north-western Pacific Acropora digitifera. Evol. Appl. 13, 1923–1938. doi: 10.1111/eva.12944
Setoguchi, H., Kyoda, S., and Maeda, Y. (2009). Geographical variation in leaflet morphology of Cycas revoluta (Cycadaceae) on the Ryukyu Islands. J. Phytogeogr. Taxon. 57, 1–6. doi: 10.25100/socolen.v46i2.7111
Shen, C.-F., Hill, K. D., Tsou, C.-H., and Chia-Jui, C. (1994). Cycas taitungensis CF Shen, KD Hill, CH Tsou & CJ Chen, sp. nov.(Cycadaceae), a new name for the widely known cycad species endemic in Taiwan. Bot. Stud. 35, 133–140.
Skeels, A., and Cardillo, M. (2019). Reconstructing the geography of speciation from contemporary biodiversity data. Am. Nat. 193, 240–255. doi: 10.1086/701125
Smith, J. F., Ooi, M. T. Y., and Clark, J. L. (2018). Incipient speciation in a neotropical Gesneriaceae: Columnea kucyniakii is nested within C. strigosa. Plant Syst. Evol. 304, 511–519. doi: 10.1007/s00606-018-1502-7
Stevenson, D. W. (2020). Observations on vegetative branching in cycads. Int. J. Plant Sci. 181, 564–580. doi: 10.1086/708812
Stewart, P., and Globig, S. (2011). Vascular Plants and Paleobotany. Boca Raton, FL: CRC Press.
Stopes, M. C. (1910). Adventitious budding and branching in Cycas. New Phytol. 9, 235–241. doi: 10.1111/j.1469-8137.1910.tb05571.x
Thakur, M., Schättin, E. W., and Mcshea, W. J. (2018). Globally common, locally rare: revisiting disregarded genetic diversity for conservation planning of widespread species. Biodivers. Conserv. 27, 3031–3035. doi: 10.1007/s10531-018-1579-x
Thiers, B. M. (2019). Index Herbariorum: A Global Directory of Public Herbaria and Associated Staff. Bronx, NY: NYBG.
Thomas, J. R., Gibson, D. J., and Middleton, B. A. (2005). Water dispersal of vegetative bulbils of the invasive exotic Dioscorea oppositifolia L. in southern Illinois1. J. Torrey Bot. Soc. 132, 187–196. doi: 10.3159/1095-5674(2005)132[187:wdovbo]2.0.co;2
Vargas, P., Arjona, Y., Nogales, M., and Heleno, R. H. (2015). Long-distance dispersal to oceanic islands: success of plants with multiple diaspore specializations. AoB Plants 7:lv073. doi: 10.1093/aobpla/plv073
Vargas, P., Nogales, M., Jaramillo, P., Olesen, J. M., Traveset, A., and Heleno, R. (2014). Plant colonization across the Galápagos Islands: success of the sea dispersal syndrome. Bot. J. Linn. 174, 349–358. doi: 10.1111/boj.12142
Walters, T., and Osborne, R. (2004). Cycad Classification: Concepts and Recommendations. Wallingford: CABI.
Wang, F. (1996). Cycads in China. Guangzhou: Guangdong Science and Technology Press.
Wang, J., and Caballero, A. (1999). Developments in predicting the effective size of subdivided populations. Heredity 82, 212–226. doi: 10.1038/sj.hdy.6884670
Wang, P., Li, Q., and Li, C.-F. (2014). Geology of the China Seas. Dev. Mar. Geol. 6, 469–553.
Wang, X. H., Wu, W., and Jian, S. G. (2019). Transcriptome analysis of two radiated Cycas species and the subsequent species delimitation of the Cycas taiwaniana complex. Appl. Plant Sci. 7:e11292. doi: 10.1002/aps3.11292
Waters, J. M., Fraser, C. I., and Hewitt, G. M. (2013). Founder takes all: density-dependent processes structure biodiversity. Trends Ecol. Evol. 28, 78–85. doi: 10.1016/j.tree.2012.08.024
Wegmann, D., Leuenberger, C., Neuenschwander, S., and Excoffier, L. (2010). ABCtoolbox: a versatile toolkit for approximate Bayesian computations. BMC Bioinformatics 11:116. doi: 10.1186/1471-2105-11-116
Whitelock, L. M. (2002). The Cycads. Portland: Timber Press.
Wiley, E. O. (1978). The evolutionary species concept reconsidered. Syst. Zool. 27, 17–26. doi: 10.2307/2412809
Wolf, J. B., and Ellegren, H. (2017). Making sense of genomic islands of differentiation in light of speciation. Nat. Rev. Genet. 18, 87–100. doi: 10.1038/nrg.2016.133
Wu, C. I. (2001). The genic view of the process of speciation. J. Evol. Biol. 14, 851–865. doi: 10.1098/rstb.2008.0078
Wu, L.-W., Yen, S.-H., Lees, D. C., and Hsu, Y.-F. (2010). Elucidating genetic signatures of native and introduced populations of the Cycad Blue, Chilades pandava to Taiwan: a threat both to Sago Palm and to native Cycas populations worldwide. Biol. Invasions 12, 2649–2669. doi: 10.1007/s10530-009-9672-4
Xiao, L.-Q., and Moller, M. (2015). Nuclear ribosomal ITS functional paralogs resolve the phylogenetic relationships of a late-Miocene radiation cycad Cycas (Cycadaceae). PLoS One 10:e0117971. doi: 10.1371/journal.pone.0117971
Yamamoto, Y. (1928). Supplementa Iconum Plantarum Formosanarum. Taihoku: The Department of Foresty, Goverment Research Institute.
Yang, J.-D., Lin, H.-C., and Liu, H.-Y. (2017). The Red List of Vascular Plants of Taiwan. Taiwan: Endemic Species Research Institute, Forestry Bureau, Council of Agriculture, Executive Yuan and Taiwan Society of Plant Systematics.
Yang, S.-F., Komaki, S., Brown, R. M., and Lin, S.-M. (2018). Riding the Kuroshio Current: stepping stone dispersal of the Okinawa tree lizard across the East Asian Island Arc. J. Biogeogr. 45, 37–50. doi: 10.1111/jbi.13111
Zachos, F. E. (2018). (New) Species concepts, species delimitation and the inherent limitations of taxonomy. J. Genet. 97, 811–815. doi: 10.1007/s12041-018-0965-1
Zhang, J., Kapli, P., Pavlidis, P., and Stamatakis, A. (2013). A general species delimitation method with applications to phylogenetic placements. Bioinformatics 29, 2869–2876. doi: 10.1093/bioinformatics/btt499
Zheng, X., Li, A., Kao, S., Gong, X., Frank, M., Kuhn, G., et al. (2016). Synchronicity of Kuroshio Current and climate system variability since the Last Glacial Maximum. Earth Planet. Sci. Lett. 452, 247–257. doi: 10.1016/j.epsl.2016.07.028
Zheng, Y., Liu, J., Feng, X., and Gong, X. (2017). The distribution, diversity, and conservation status of Cycas in China. Ecol. Evol. 7, 3212–3224. doi: 10.1002/ece3.2910
Zonneveld, B. J. M. (2011). Genome sizes for all genera of the Cycadales compared with earlier cladistic analysis. Plant Biol. 14, 252–256.
Keywords: continental island, Cycas, Kuroshio, long distance dispersal, speciation, species concept
Citation: Chang J-T, Chao C-T, Nakamura K, Liu H-L, Luo M-X and Liao P-C (2022) Divergence With Gene Flow and Contrasting Population Size Blur the Species Boundary in Cycas Sect. Asiorientales, as Inferred From Morphology and RAD-Seq Data. Front. Plant Sci. 13:824158. doi: 10.3389/fpls.2022.824158
Received: 29 November 2021; Accepted: 15 March 2022;
Published: 09 May 2022.
Edited by:Angelica Cibrian-Jaramillo, Instituto Politécnico Nacional de México (CINVESTAV), Mexico
Reviewed by:Loreta Brandão de Freitas, Federal University of Rio Grande do Sul, Brazil
Diego Garfias-Gallegos, National Polytechnic Institute of Mexico (CINVESTAV), Mexico
Christine D. Bacon, University of Gothenburg, Sweden
Copyright © 2022 Chang, Chao, Nakamura, Liu, Luo and Liao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Pei-Chun Liao, email@example.com