Genetic Divergence and Relationship Among Opisthopappus Species Identified by Development of EST-SSR Markers

Opisthopappus Shih is an endemic and endangered genus restricted to the Taihang Mountains that has important ornamental and economic value. According to the Flora Reipublicae Popularis Sinicae (FRPS, Chinese version), this genus contains two species (Opisthopappus longilobus and Opisthopappus taihangensis), whereas in the Flora of China (English version) only one species O. taihangensis is present. The interspecific phylogenetic relationship remains unclear and undefined, which might primarily be due to the lack of specific molecular markers for phylogenetic analysis. For this study, 2644 expressed sequence tag-simple sequence repeats (EST-SSRs) from 33,974 unigenes using a de novo transcript assembly of Opisthopappus were identified with a distribution frequency of 7.78% total unigenes. Thereinto, mononucleotides (1200, 45.39%) were the dominant repeat motif, followed by trinucleotides (992, 37.52%), and dinucleotides (410, 15.51%). The most dominant trinucleotide repeat motif was ACC/GGT (207, 20.87%). Based on the identified EST-SSRs, 245 among 1444 designed EST-SSR primers were selected for the development of potential molecular markers. Among these markers, 63 pairs of primers (25.71%) generated clear and reproducible bands with expected sizes. Eventually, 11 primer pairs successfully amplified all individuals from the studied populations. Through the EST-SSR markers, a high level of genetic diversity was detected between Opisthopappus populations. A significant genetic differentiation between the O. longilobus and O. taihangensis populations was found. All studied populations were divided into two clusters by UPGMA, NJ, STRUCTURE, and PCoA. These results fully supported the view of the FRPS, namely, that O. longilobus and O. taihangensis should be regarded as two distinct species. Our study demonstrated that transcriptome sequences, as a valuable tool for the quick and cost-effective development of molecular markers, was helpful toward obtaining comprehensive EST-SSR markers that could contribute to an in-depth assessment of the genetic and phylogenetic relationships between Opisthopappus species.


INTRODUCTION
Opisthopappus Shih is a perennial herbaceous genus of the Asteraceae family, which is endemic to and endangered within China, that grows primarily in the cliff cracks and rock gaps of the Taihang Mountains (Shih, 1979;Ding and Wang, 1997;Wu et al., 1999;Hu et al., 2008). This genus exhibits excellent drought tolerance and leanness resistance that make it an important wild Chrysanthemum resource (Chai et al., 2018). Being a diploid, Opisthopappus is a 2A asymmetrical karyotype that has an intimate relationship with the Artemiisinae subtribe (Zhao et al., 2009). AFLP analysis has revealed that Opisthopappus is clustered with a majority of Ajania taxa, located on the basal branch of Dendranthema and Ajania. This indicates that Opisthopappus might have been a more primitive group during the eastward spread and evolutionary process of Dendranthema ancestor groups (Zhao et al., 2010).
According to the Flora Reipublicae Popularis Sinicae (FRPS) 1 (Shih and Fu, 1983), Opisthopappus contains two species, Opisthopappus longilobus Shih and Opisthopappus taihangensis (Ling) Shih. The morphological differentiation between these two species presents mainly on the leaves and bracts. For O. longilobus, which has hairless leaves, there is one pinnatisect, except the basal stem leaves, and a pair of bracteal leaves beneath the involucres. For O. taihangensis, there is appressed puberulent on both surfaces of the leaves, two pinnatisect stem leaves, and no bracteal leaves. Nevertheless, according to the Flora of China 2 , O. taihangensis was the only one species in the genus of Opisthopappus and O. longilobus was merged into it. Its basal leaves are one-or two-pinnatisect, both surfaces sparsely pubescent or glabrous, and the stem leaves are similar to the basal leaves, uppermost leaves pinnatifid.
Based on the observation on the pollen morphology of Opisthopappus, the results showed the pollen of O. longilobus at spine length and density, ostiole, density, colpus depth is more highly evolved than that of O. taihangensis (Gao et al., 2011). The micromorphological characteristics of leaves revealed that these two species differed in their characteristics and density of glandular hairs, size and density of stomata, anticlinal wall patterns, wax characteristics, as well as the size and density of salt vesicles (Jia and Wang, 2015). A comparison between the biological characteristics of O. longilobus and O. taihangensis in their natural habitats and artificial populations indicated stable inheritable characters and significant differences in their leaf morphologies and other features (Hu et al., 2008). And above all, the evolutionary relationships and taxonomic status of the species within Opisthopappus remain debated and inconclusive. Further investigations and analysis are required to definitively determine whether O. longilobus and O. taihangensis should be regarded as two distinct species or a single species.
Our previous research endeavored to explore the relationships between Opisthopappus species through various molecular 1 Chinese version: http://www.iplant.cn/frps2019/frps?id=%e5%a4%aa%e8%a1% 8c%e8%8f%8a%e5%b1%9e 2 English version: http://www.efloras.org/florataxon.aspx?flora_id=2&taxon_id= 200163 markers, such as Inter-Simple Sequence Repeat (ISSR), Sequence Related Amplified Polymorphism (SRAP), Internal Transcribed Spacer (ITS), and chloroplast gene sequences (Wang, 2013;Wang andYan, 2013, 2014;Wang et al., 2015). However, cluster analysis revealed that neither O. longilobus nor O. taihangensis formed a monophyletic clade (Wang, 2013;Wang and Yan, 2013), particularly, some populations of O. longilobus always gathered with O. taihangensis populations. These cluster analysis results indicated that the two species would seemingly be merged into one species, or the limited gene data offered from the above molecular markers might not be sufficient for phylogeny analysis. Zhang et al. (2019) also pointed that the phylogenetic relationships and genetic components between species were unclear due to the limited detection capacity of specific or unique nuclear gene sequences. Therefore, the development of transferable and orthologous molecular markers between phylogenetically related species of Opisthopappus is essential and critical.
Recently, RNA sequencing (RNA-Seq) has been employed in the development of a high-throughput sequencing technology that has the ability to profile the entire transcriptome of any organism, particularly for non-model plants that lack sequenced genomic data (Gross et al., 2013;Wu et al., 2014;Du et al., 2015;Zhang et al., 2019). Meanwhile, transcriptome sequencing holds great potential as a platform for the generation of molecular markers . Further, de novo RNA-Seq assembly can provide a cost-effective approach to identify molecular markers, such as microsatellite markers (SSRs) and singlenucleotide polymorphisms (SNPs) (Davey et al., 2011;Karan et al., 2012;Chen et al., 2013;Zhang et al., 2019). What is more, both of SSRs and SNPs markers are codominant and the reflection of specific segment or sequence information, which are recommended by UPOV-BMT molecular testing guidelines as optimal molecular markers for authenticity identification of varieties and database construction (UPOV-BMT, International Union for the Protection of New Varieties of Plants-Molecular Biology Technology Working Group 2007 3 ). Based on the predicted SNPs, we have successfully developed eight SNP markers for population genetic analysis in Opisthopappus (Chai et al., 2018). However, with the help of the new technology and effective method, more information and further evidences are needed for more comprehensively and accurately solving the problem on the taxonomy of Opisthopappus.
Simple sequence repeats derived from ESTs (EST-SSRs) identified from transcribed RNA sequences typically contain mono-, di-, tri-, tetra-, penta-, or hexa-nucleotide units, and possess abundant information, codominance, and loci specificity (Song et al., 2012;Zhou et al., 2016). With codominant features between distantly related species, EST-SSRs cross-transferability can shed new light on the evolution of plant genomes, changes in gene location, and genome organization (Bazzo et al., 2018;Zhang et al., 2019). Moreover, this transferability provides a costeffective source of markers for related species, which is important for taxa with low microsatellite frequencies, or for those whose microsatellites are difficult to isolate (Bazzo et al., 2018;Zhang et al., 2019). Previous transcriptome sequences in Opisthopappus, which generated and identified 33,974 unigenes and more than 3410 differentially expressed transcripts (DETs), provided crucial genetic resources and sequence data toward the development of EST-SSR markers (Chai et al., 2018). Therefore, in this study, we developed EST-SSR markers from the Opisthopappus transcriptome and characterized their frequency, distribution, and function, which were further employed to elucidate the genetic differentiation and phylogenetic relationships between O. longilobus and O. taihangensis.

Development of EST-SSR Markers From Opisthopappus Transcriptome
The cDNA libraries were constructed using Opisthopappus tissue samples, and sequenced based on next-generation sequencing (NGS) (Chai et al., 2018). The obtained raw reads were preprocessed to filter adaptor sequences, low-quality sequences, and reads with qualities of less than Q30 using the FASTX toolkit (Gordon and Hannon, 2010). Following the removal of raw reads, high-quality clean reads were assembled de novo into contigs and unigenes using a short read assembling program called Trinity RNA-Seq Assembler (Haas et al., 2013). The assembled unigene sequences were directly identified by using MicroSatellite (MISA) identification tools (Thiel et al., 2003), for which the parameters were set at mono-, di-, tri-, tetra-, penta-, and hexa-nucleotide motifs as a minimum repeat number of 12, 6, 5, 5, 4 and 4, respectively. Based on the MISA results, the EST-SSR primers were designed using Primer3 software (Untergasser et al., 2012).

Plant Materials for Genetic Analysis
Nineteen populations of Opisthopappus were sampled, including 12 populations of O. taihangensis and 7 populations of O. longilobus, which were collected during August and September of 2017 that essentially covered their distribution ranges ( Table 1 and Figure 1). Each individual from the same population was collected from different locations at least 5 m apart. Fresh young leaves were selected in each of the sample sites, where 10-15 individuals from each population were collected, for a total of 214 samples from 19 populations. The collected samples were quickly dried with silica gel and stored at 20 • C for future use. A global positioning system (GPS) was employed to demarcate each sample site and record the longitude, latitude, and altitude of each population.

DNA Extraction and PCR Amplification
The genomic DNA was extracted from leaves though modified CTAB methods (Doyle et al., 1987). The quantity and quality of the DNA samples were determined using a NanoDrop ND1000 spectrophotometer (Thermo Scientific, Waltham, MA, United States) and 0.8% agarose gel electrophoresis, which were then diluted to 30 ng/µl and stored at −20 • C until further use.
Two hundred twenty-five EST-SSR primer pairs were initially designed based on the potential SSR loci identified in transcriptome. Subsequently, 63 primer pairs were selected to amplify 19 random DNA samples (one for each population) for the pre-experiments (primer screening). These abided by the following principles: the primer length was controlled between 18 and 25 bp, the GC content was from 40 to 60%, and the annealing temperature was from 50 to 60 • C. Hairpin structures, dimers, hat structures, and mismatches were avoided as much as possible.
The primers were synthesized by Shanghai Sangon Biological Engineering Technology (Shanghai, China). Among these 63 primer pairs, 11 primers were ultimately selected based on high transferability, polymorphism, and repeatability to generate clear and reproducible bands with expected sizes, which were used to further amplify all individuals of the Opisthopappus populations ( Table 2).

Data Analysis
The presence or absence of each SSR fragment was coded as "1" and "0, " respectively, and only clear and reproducible bands were considered. For the neutrality test of microsatellite alleles, we used the Bayesian-based outlier analysis of the multinomial-Dirichlet model (BayeScan analysis) to assess the F ST distribution among all alleles by BAYESCAN v 2.1 (Beaumont and Nichols, 1996;Hsiung et al., 2017).
Genetic diversity parameters, including the observed number of alleles (Na), effective number of alleles (Ne), Nei's genetic diversity (H), Shannon's information index (I), and the percentage of polymorphic bands (% PPB) were calculated using the POPGENE v 1.31 program (Yeh et al., 1999).
A principal coordinate analysis (PCoA) was constructed using GENALEX v 6.5 (Peakall and Smouse, 2012) and visualized with the R package "ggplot2" (Wickham, 2015). The unweighted pair-group method with arithmetic means (UPGMA) tree and neighbor-joining (NJ) trees were constructed using MEGA v 7 (Kumar et al., 2016), through which the probabilities for each node was assessed by PUAP with 1000 replicates bootstrap (Swofford, 2000). The analysis of molecular variance (AMOVA) was employed to investigate the total genetic variation between and within species, and among and within populations using GENALEX v 6.5 (Peakall and Smouse, 2012).
Bayesian clustering analysis (BCA) was utilized to examine the similarity and divergence of genetic components between populations, and was performed using STRUCTURE 2.3.4 (Hubisz et al., 2009). The analysis was run for K-values ranging from 1 to 19 possible clusters with 15 independent runs each and a burn-in of 10,00,000 iterations and 20,00,000 subsequent MCMC steps to evaluate consistency. The combination of an admixture and a correlated-allele frequencies model was used. The posterior probability number of grouping was evaluated  using K by STRUCTURE HARVESTER 0.6.94 (Earl, 2012). These runs were aligned and summarized using CLUMPP 1.1.2 (Jakobsson and Rosenberg, 2007), and the visualization of the results was plotted using DISTRUCT 1.1 (Rosenberg, 2004). Pearson correlation analysis was used to test the correlations between the genetic parameters and geographic factors (longitude, latitude, altitude). Correlations between pairwise genetic distance and adjusted pairwise geographic distances were calculated with GENALEX v 6.5 (Peakall and Smouse, 2012) using the mantel test with 9999 permutations. The generalize linear regression model (GLRM) was conducted between the two distances (Nelder and Wedderburn, 1972).

Characterization and Distribution of EST-SSR Markers
Through the removal of adaptors, approximately 34.9 million clean reads were obtained. The average Q30 was >95%, which reflected the high quality of the sequencing data, whereas the percentage of unknown nucleotides was only ∼0.01%. The obtained unigenes were 33,974 having paired-end reads with average lengths of 801 bp. Among these unigenes, the lengths of 14,810 unigenes (43.59%) ranged from 200 to 500 bp, and 1893 (5.57%) unigenes were longer than 2 kb. Of the 33,974 unigenes, all were successfully annotated in the NR and GO databases. About 990 of the unigenes (2.91%) were simultaneously annotated in all of the databases (NR, NCBI nonredundant protein sequences; Go, Gene Ontology; KEEG, Kyoto Encyclopedia of Genes and Genome; eggNOG, evolutionary genealogy of genes: non-supervised Orthologous Groups; Swiss-Prot) (Supplementary Table S1).

Genetic Analysis of Opisthopappus Species
There were total 420 alleles generated by 11 EST-SSR primers across 214 samples. Eight alleles were removed with the neutrality tests by BAYESCAN 2 and the remaining 412 alleles used for further genetic population assessment and divergence tests (Supplementary Figure S1).
For the genus Opisthopappus, the PPB was 100%. Among all populations, the observed allele number (Na) ranged from 0.35 to 1.539, with an average of 1.13805. The expected allele number (Ne) ranged from 1.045 to 1.377, with an average of 1.2505. Nei's genetic diversity (H) ranged from 0.0315 to 0.2326, with an average of 0.128185. Shannon diversity index (I) was from 0.056 to 0.35, with an average of 0.24435.
In general, O. longilobus had a relatively higher level of genetic diversity (Na = 1.937, Ne = 1.356, H = 0.137, I = 0.327, Analysis of molecular variance ( Table 4) was employed to evaluate variance components between and within different populations and species. The results revealed that 67% of the total variation occurred within Opisthopappus populations, whereas 33% of the total variation existed between the genus populations. For each species, variations also significantly occurred within populations of both O. longilobus (62%) and O. taihangensis (74%). This was consistent with the hierarchical AMOVA results, which revealed that most of the genetic variation resided within populations (65%) rather than among groups (7%).
Clustering analysis well distinguished all populations into two major clusters in both UPGAM and NJ tree (Figure 2): O. taihangensis cluster and O. longilobus cluster. The results we observed from the PCoA were in agreement with the clustering analysis (Figure 3). The first two principal coordinates detected 64.14% of the total variation. The first coordinate axis explained 43.04% of the total genetic components, whereas the second coordinate axis explained 21.10% of the genetic components. A moderate but clear separation between the populations of different species was observed (Figure 3).
The genetic structures of 19 populations of Opisthopappus were analyzed using STRUCTURE software. Based on maximum-likelihood and K ( K)-values, the optimal number of groups was two (K = 2). As shown in Figure 4, the O. longilobus species could be separated from the O. taihangensis species, although some individuals were mixed between populations. When K = 3 (the second-best K-value), 7 populations of O. longilobus still gathered into one group, while 12 populations of O. taihangensis were divided into two subgroups. Specifically, the WWS, YTS, WML, BJY, LSZT, and XYG populations were gathered together, whereas the other 6 populations FXF, SNS, GS, QLX, JYS, and DSC were clustered into another subgroup (Figure 4).
The Mantel test was employed to investigate the correlations between genetic distances and geographic distances (Figure 5). Positive correlations between the genetic and geographic distances were clearly found in the GLRM (r = 0.4565, P = 0.001; y = 0.0011x + 0.2117, R 2 = 0.2082) (Figure 5), which indicated that the genetic differentiation of the populations was linearly correlated with geographic distances.

The Development of EST-SSRs Markers Based on Opisthopappus Transcriptome Database
Due to the lack of sufficient genetic information and effective molecular marker systems, the phylogenetic relationship between Opisthopappus (Asteraceae) species remains unclear and unresolved to a certain extent. SSR markers are considered to be one of the most important marker systems for plant genetic studies, with high polymorphisms, high abundance, co-dominance, and genome-wide distribution (Thakur and Randhawa, 2018;Zhang et al., 2018). Compared with the genomic SSRs, EST-SSRs are effective and useful markers for non-model plants, which are easily transferred between related species owing to these regions being more evolutionarily conserved than non-coding sequences (Varshney et al., 2005;Ellis and Burke, 2007;Tahan et al., 2009;Bazzo et al., 2018;Park et al., 2019).
The application of EST-SSR markers for genetic diversity, phylogenetic analyses, the construction of genetic maps, and the marker-assisted selection of important traits has recently been reported for several plant species (Dutta et al., 2011;Jo et al., 2015;Huang et al., 2016;Vatanparast et al., 2016;Taheri et al., 2019). At present, the transcriptome sequencing technology (RNA-Seq) is a rapid, reliable, and cost-effective tool for the characterization of gene content, and the identification of polymorphic markers in non-model plants (Taheri et al., 2019;Zhang et al., 2019). An increasing number of successful examples have supported the  strategy of using transcriptome data to predict SSR molecular markers (Zalapa et al., 2012;Han et al., 2018;Taheri et al., 2019). For our work, among 33,974 assembled unigene sequences, 2644 potential EST-SSRs were identified, which represented ∼7.78% of the transcriptomic sequences (Supplementary Table  S2), and the distribution density was one SSR per 10.30 kb. The SSRs frequency and distribution density were higher than some species, such as Argyranthemum broussonetii (2.3%, 27 kb) (White et al., 2016) and Zingiber officinale (2.7%, 25.2 kb) (Awasthi et al., 2017), and lower than Arachis hypogaea (17.7%, 3.3 kb) , Curcuma longa (14.6%, 5.3 kb; 14.9%, 5.2 kb; 20.5%, 4.8 kb) (Annadurai et al., 2013), and Curcuma alismatifolia (12.5%, 6.6 kb) (Taheri et al., 2019). Differences in the frequency of SSRs in ESTs could be partially attributed to the different genetic basis of various plant species, the size of the unigene assembly dataset, SSR search criteria, sequence redundancy, as well as the mining tools utilized (Kumpatla and Mukhopadhyay, 2005;Chen et al., 2015;Wei et al., 2016;Liu et al., 2018;Taheri et al., 2019).
Based on the number of bases, five different repeat motifs from mono-to penta-nucleotide were identified, in which the most abundant SSR motifs were mononucleotide repeats (45.39%), followed by trinucleotide repeats (37.52%), and dinucleotide repeats (15.51%). In contrast, hexa (1.48%) and penta (0.11%) were rare (Supplementary Table S2). This result was similar to previous findings of di-and tri-nucleotide motifs, which were reported as the more frequent SSR motif types within the transcriptome sequences of many other plants (Zheng et al.,   The K plot shows that K = 2 got the highest K-value and K = 3 got the second-best K-value, meaning that the most probable grouping number could be two or three. (B) Estimated genetic structure for K = 2 obtained for Opisthopappus. (C) Estimated genetic structure for K = 3 for Opisthopappus.
2013; Awasthi et al., 2017;Park et al., 2019;Taheri et al., 2019). In many organisms, the extensive distribution of trinucleotide repeats in coding sequences is an indicator of the effects of selection and evolution. Because open-reading frames do not disturb insertions and deletions within translated regions, this signifies that these SSRs do not change the coding frame of the gene regions (Feng et al., 2012;Bazzo et al., 2018). When there were variations in SSR lengths of other repetitions, they would lead to frame shifts and induce negative mutations, which might restrict the development of other motif types (Barrero et al., 2011;Park et al., 2019;Taheri et al., 2019).
Abundant trinucleotide repeats generally encoded for threonine, isoleucine, lysine, asparagine, and serine, respectively, and except for isoleucine, these were primarily small hydrophilic amino acids. It is possible that the codon repeats corresponding to hydrophilic amino acids more easily tolerated selection pressures, which might eliminate codon repeats that encoded hydrophobic and basic amino acids (Katti et al., 2001;Han et al., 2018;Taheri et al., 2019). Consequently, the high occurrence level of these motifs was significant, as amino acids generated by them were observed in proteins to a high degree (Bazzo et al., 2018;Taheri et al., 2019;Zhang et al., 2019). The Opisthopappus species typically grows within the cliff cracks and rock gaps of the Taihang Mountains, have good drought tolerance, and are leanness-resistant (Hu et al., 2008;Chai et al., 2018). Abundant repetitions of the acids from the trinucleotide motif might explain these characteristics.
The 63 designed EST-SSR primers were obtained and verified, which offered an informative and applicable approach for the evaluation of genetic relationships between Opisthopappus species. Of the 63 primer pairs, 11 (17.5%) of the pair primers produced clear bands. This rate was relatively lower than that of previously reported species (Feng et al., 2017;Bazzo et al., 2018;Taheri et al., 2019;Zhang et al., 2019). Since there were no specificity or assembly errors following careful checking, the other 52 pair primers failed to generate the expected amplification size or non-amplification PCR products, which may have been due to the presence of introns and indels.
Meanwhile, the EST-SSR markers developed from Opisthopappus transcriptome data exhibited good transferability across different Opisthopappus species, indicating that these EST-SSR markers were useful tools for further genetic diversity analysis for this genus. The transferable nature of the EST-SSR markers between related species extended their usefulness in plant genetic studies (Varshney et al., 2005;Mathithumilan et al., 2013;Jia et al., 2016;Wang et al., 2018;Park et al., 2019).

Genetic Differentiation and Phylogenetic Relationships Between Opisthopappus Species
Opisthopappus, which is an endemic and endangered genus that is geographically distributed in the Taihang Mountains, includes two species, according to the FRPS, while these two species are merged into a single species in the Flora of China. Therefore, the exact phylogenetic status of the Opisthopappus species remains under debate. Although the systematic classification and phylogenic research on Opisthopappus has been reported previously (Gao et al., 2011;Wang and Yan, 2013;Jia and Wang, 2015), O. longilobus populations are not well separated from O. taihangensis populations.
A number of reports have been published regarding the interspecific relationships between Opisthopappus species based on molecular markers (Wang, 2013;Wang and Yan, 2013); however, such research has been limited due to the lack of specific genomic information on Opisthopappus. Compared with other EST-SSR markers of plants (Kudapa et al., 2012;Onda et al., 2015;Wei et al., 2016;Li et al., 2018), the developed Opisthopappus transcriptomic SSR markers showed a high level of genetic diversity (Table 2), which indicated that the identified microsatellite loci possessed high allele numbers and heterozygosity, which provided efficient genetic markers for Opisthopappus populations.
Genetic diversity and fitness play a critical role in influencing the genetic structures of populations. High genetic diversity within populations is a successful strategy for adapting to various habitats and environmental conditions (Booy et al., 2000), which also increases the capacity of individuals plants, or populations thereof, to reach and colonize new habitats (Crawford and Whitney, 2010;Drummond and Vellend, 2012). Therefore, as one of the three forms of biodiversity, genetic diversity is recognized by the World Conservation Union (IUCN) as deserving conservation along with species and ecosystems (Reed and Frankham, 2003).
As an endangered perennial herb that is native to the Taihang Mountains, Opisthopappus genus populations are facing the prospect of completely disappearing in China. However, our results revealed that the genetic diversity of the populations evaluated with EST-SSR markers was more robust (Table 1). Moreover, O. longilobus populations exhibited a relatively higher level of genetic diversity than that of O. taihangensis populations ( Table 1). O. longilobus primarily distributed along the northern regions of the Taihang Mountains, which has relatively low temperatures and high precipitation (Hu et al., 2008;Gong, 2010). In contrast, O. taihangensis is mainly distributed in the Southern Taihang Mountains (Hu et al., 2008;Gong, 2010). In comparison, the environmental conditions in which O. longilobus grown should be relatively harsher than that for O. taihangensis. Thus, the relatively inhospitable surroundings of O. longilobus might have cumulatively imbued it with high levels of genetic diversity.
In addition, based on the newly developed EST-SSR markers, it was found that there was a significant differentiation between Opisthopappus species (Figures 2-4). This differentiation might have resulted from long-term adaptions within a heterogeneous environment, as there are many gully valleys and gulches in the Taihang Mountains. Meanwhile, the complex topographies and diverse localized climates have also led to the varied surroundings of different species or populations.
The mantel test and GLRM ( Figure 5) revealed a significant correlation between genetic and geographical distances; thus, geological impacts may have played a major role in driving and governing the genetic divergence of Opisthopappus. As new mountains formed, a long established population might have been broken up into spatially isolated subpopulations, which gradually accumulated genetic differences via adaptations to the local environment caused by this geographic segmentation (Meng et al., 2015). Originally, with changing monsoon patterns in East Asia during the Middle Miocene, the genus Opisthopappus have begun to differentiate into two major lineages: O. longilobus and O. taihangensis (Ye et al., 2019) and might have been homogeneously distributed across Northern China. Subsequently, with the uplift of Taihang Mountains and the glacial-interglacial cycles during Quaternary period (Gong, 2010;Zhang and Li, 2014), the population was spatially and geologically decollated into several different subpopulations situated in specific and distinct ecological niches, which initiated genetic differentiation and variation.
Furthermore, the Northern Taihang Mountains belong to a temperate continental semi-humid monsoon climate, while the Southern Taihang Mountains have a warm temperate continental monsoon climate (Wu et al., 1999;Gong, 2010). In the southern regions, the annual mean temperature is 12.7 • C and the annual mean precipitation is 606.4 mm, whereas in the northern domains, the annual mean temperature is ∼10 • C and the annual mean precipitation is 700 mm. Climatic and environmental diversity would result in geographical and ecological isolation between the species/populations, and promote the differentiation and speciation of O. longilobus and O. taihangensis. Therefore, with this accumulation of genetic differences, these isolated subpopulations began to evolve into genetically distinct lineages as the result of their adaptations to localized environmental conditions, to form the speciation pattern of Opisthopappus.
Clustering analysis (Figure 2) showed that all Opisthopappus populations were grouped into two major clusters, where O. longilobus populations gathered into one cluster, and O. taihangensis congregated into another (Figure 2). PCoA results also showed an obvious difference in the two groups (Figure 3), which was consistent with the clustering analysis (Figure 2). Moreover, BCA suggested that the best grouping number (K) was 2, and the second-best K was 3 based on the K (Figure 4). All studied O. longilobus populations were gathered into one group in both K = 2 and K = 3. Nevertheless, the populations of O. taihangensis were inferred into two subgroups when K = 3. Therefore, we inferred that the O. taihangensis may retain a divergence potentiality to differentiation into distinct groups. The two inferences (K = 2 or K = 3) were similar and consistently indicated that the populations should be divided into two major groups, namely O. taihangensis and O. longilobus.

CONCLUSION
For the current study, we provided a number of EST-SSR markers to elucidate the genetic relationships of Opisthopappus species, based on a de novo RNA-Seq assembly. To the best of our knowledge, this is the first attempt to develop SSR markers for this endemic and endangered genus. In the light of the above and previous results, the strategy for obtaining EST-SSRs and SNPs from the transcriptome proved to be efficient. Using the designed EST-SSRs and SNPs markers, a significant genetic differentiation was revealed between O. longilobus and O. taihangensis (Chai et al., 2018).
All of the studied Opisthopappus populations were separated into two distinct groups, which corresponded to O. longilobus and O. taihangensis species. These results indicated that O. longilobus and O. taihangensis should be regarded as two independent species, which was well supported by the FRPS. The developed and identified EST-SSR markers were shown to be efficient and significantly contributed to the genetic and evolutionary study of Opisthopappus.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the NCBI and the GenBank accession numbers were SUB4002238 and PRJNA471103.

AUTHOR CONTRIBUTIONS
MC, HY, and ZW performed the experiments. HY, ZW, JW, and YG collected and analyed the transcriptome data. WH, EZ, WR, and HZ analyed the population genetics of Opishopappus species. MC, HY, and YW wrote the manuscript. MC, HY, YZ, GS, and YW edited the manuscript. All authors contributed to manuscript revision, read and approved the submitted version.