Genetic architecture of microhabitat adaptation traits in a pair of sympatric Primulina species

Abstract

Exploring the genetic basis of adaptive divergence at fine spatial scales can broaden our basic understanding of evolution and how organisms may adapt to changing environments in the future. Cave-associated microhabitats provide a unique opportunity to gain insight into microgeographic adaptation. We studied the genetic architecture of microhabitat-related divergence in flower phenology and leaf traits between two sister species of Primulina, P. depressa and P. danxiaensis, which live in sympatry but occupy contrasting microhabitats. We identified 40 significant quantitative trait loci (QTLs) associated with the interspecific differences in these microhabitat adaptation traits. Flowering time was controlled by one major-effect and six minor-effect QTLs, while leaf traits were influenced by 9–12 QTLs of small to moderate effect. The genetic architecture of the flowering time and the specific leaf area was genetically independent of other traits. Our results suggest that microhabitat adaptation in sympatric populations of Primulina differs according to different traits, with leaf traits diverging with the accumulation of many small changes and flowering phenology being driven by major effect variance.

Keywords:

• Primulina
• microhabitat adaptation
• QTL mapping
• species divergence
• leaf traits
• flowering phenology

Introduction

Adaptive divergence at fine spatial scales (i.e. microgeographic adaptation) has broad ecological and evolutionary implications [1]. In the face of local environmental conditions, plants cannot escape to a favorable habitat but cope with environmental factors directly. This is particularly relevant for flowering plants confined to terrestrial habitat islands, such as inselbergs, granite outcrops and caves, where species occur in small, restricted landmasses and are generally isolated from a surrounding community [2]. This phenomenon is apparent in southern China, where both karst and Danxia landforms (a red terrigenous landform named after Danxia Mountain in Guangdong, China) have generated abundant cave and cave-associated habitats [3, 4]. These microhabitats support a remarkably high level of species endemism in flowering plants, as observed for Begonia, Pilea, and Primulina. Although the species diversity of individual caves is usually low, the plant composition shows considerable variation from one cave to another and, hence, low α diversity but high β diversity across the fragmented landscape [5]. This has led to the recognition of southern China as one of the world’s centers of plant diversity [6]. Both ‘South China Karst’ and ‘China Danxia’ are special landforms and are rich in biodiversity, and they have been listed as World Natural Heritage sites by UNESCO and are in great need of protection. The ability of plants to adapt to cave and cave-associated habitats may have important consequences for the generation, maintenance, and distribution of species diversity.

As highlighted previously [7], cave and cave-associated habitats often provide abiotic and biotic selection agents for a wide range of plant traits. These include pollinator limitation, shallow soil with low-water holding capacity, and reduced availability of essential plant nutrients. Flowering phenology is a key ecological trait that has been shown to be under both biotic and abiotic selection [8, 9]. It is a critical life history trait that ensures the seed production required for the survival of species [10]. In plants that depend on animal pollinators for reproduction, the flowering times have evolved alongside the emergence period of pollinators [11, 12]. In terms of abiotic factors, water deficiency results in earlier flowering to escape drought [13, 14]. Temperature is also an important selection factor for flowering time, with warmer climates resulting in earlier flowering [15].

Leaf traits are important adaptive traits that experience numerous selective pressures. In addition to their major role in photosynthesis, leaves are involved in nutrient storage, defense and stress responses. Small leaves are associated with harsh conditions such as nutrient deficiency, aridity, high light, or a combination of these factors [16]. Even though leaf shape varies less predictably across environments or biomes than size [17], it has often been reported to be strongly influenced by leaf water relations (reviewed by Nicotra et al. [16]). Moreover, the specific leaf area (SLA) is one of the most predictive traits for plant physiological strategies because of its close relationship with leaf life span, growth rate, photosynthetic rate, nutrient availability and water availability [18]. Thus, the diversity of both flower and leaf traits contributes to microhabitat adaptation. Lambrecht and Dawson (2007) demonstrated that leaf morphology and physiology influence flower traits (pollination system) through moisture availability and water usage in five plant species [19].

Although many studies have explored the ecological and evolutionary significance of flower and leaf traits in various natural systems, their genetic architecture has received less attention. Quantitative trait locus (QTL) mapping is effective in detecting genetic architecture, and it has been used widely for genetic analyses in crop species [20–22]. Crops are not suitable for the study of the genetic basis of adaptation to natural environments because their traits are often shaped from artificial selection. The genetic architecture of flower and leaf traits in wild populations provides us a candidate method to detect the evolutionary processes of plant traits in nature [23, 24].

Primulina depressa and P. danxiaensis are sister species that are narrowly endemic to northern Guangdong, south China. These two species have a sympatric distribution (<100 m) in Danxia Mountain, a National Park and UNESCO World Heritage site in Guangdong, south China. However, these two species differ in their microhabitats, with P. danxiaensis tending to occupy the twilight zone of caves and P. depressa tending to grow on outcrops (Figure S1 of Feng et al. [25]). Despite their contrasting morphological differences, phylogenetic analysis of individuals from multiple sites indicates that the P. danxiaensis population sympatric with P. depressa shows greater similarity to P. depressa than to conspecific populations [26]. Field investigations have revealed that the two have partially overlapping flowering times for about one month, and their pollinators are partially identical [27]. In a greenhouse, fertile hybrids can be generated by using P. depressa as the pollen donor. We hypothesized that phenotypic differentiation in these species in sympatry is maintained by local adaptation to local microhabitats. Identification of the number, size and distribution of genomic regions underlying phenotypic differentiation is a key step in the analysis of the genetic basis of microgeographic adaptation.

Figure 1. Fully expanded leaves of P. depressa (1), P. danxiaensis (2), the F₁ hybrid (3) and some individuals of the F₂ population (4–8).

In our previous studies of this species pair, we detected the genetic architecture of some flower and leaf traits and hybrid male fertility with QTL methods based on an F₂-interspecific cross between P. depressa and P. danxiaensis [25, 27]. For this report, we conducted QTL mapping for leaf size and shape traits (maximum leaf area and leaf length-width ratio), a leaf physiology trait (SLA), and a phenology trait (flowering time). We characterized the number, mode of action, digenic interactions, and phenotypic effects of loci. We then detected whether QTLs underlying these microhabitat adaptation traits were coincident with other flower and leaf or post-zygotic isolation trait QTLs. Our primary objective was to dissect the genetic architecture of these traits that play important roles in microhabitat adaptation in this system.

Materials and methods

Plant materials and sequencing

As previously described [27], we collected parental lines (P. depressa and P. danxiaensis) at the sympatric site of Danxia Mountain National Park. The voucher specimens were deposited at the South China Botanical Garden Herbarium (IBSC; the serial numbers for P. depressa and P. danxiaensis are DXS02 and DXS04, respectively). The parents exhibited differences in leaf morphology, leaf physiology and flowering time traits [25] (Figure 1; Table 1). The F₂ mapping population of 201 individuals was derived by self-fertilizing one F₁ plant that was descended from a cross between a paternal DXS02 individual and a maternal DXS04 individual. We propagated these plants by leaf cuttage, and planted them in a greenhouse in the South China Botanical Garden under fluorescent light to provide a 14-h day length, with a temperature of ∼26 °C. Plants were sub-irrigated as needed and fertilized weekly. F₁ hybrids, male and female parents, and 201 F₂ individuals (the mapping population) were genotyped by restriction site-associated DNA sequencing (RAD-seq). The linkage map was constructed with 195 individuals and it was of high-density with 2,484 markers spaced across 18 linkage groups. For further details, see Feng et al. [25, 27].

Table 1. Means and standard deviations for the four studied traits in P. depressa and P. danxiaensis and F₁ hybrids, and phenotypic values range of F₂ population.

Trait	P. depressa population (n = 10)	P. danxiaenesis population (n = 9)	Paternal individual	Maternal individual	F1 (n = 3)	F2 phenotypic range	t-test (two-sided)
Maximum leaf area (LA; cm2)	55.76 ± 10.76	23.56 ± 7.30	55.76 ± 10.76	23.56 ± 7.30	43.37 ± 7.45	6.01-51.72	t = 3.1824 p < 0.05
Leaf length-width ratio (LWR)	2.47 ± 0.032	0.83 ± 0.025	2.47 ± 0.032	0.83 ± 0.025	0.96 ± 0.10	0.75-1.44	t= −14.884 p < 0.001
Specific leaf area (SLA; cm^2/g)	241.78 ± 35.99	414.06 ± 37.84	241.78 ± 35.99	414.06 ± 37.84	284.31 ± 4.99	223.07-532.79	t= −2.9769 p < 0.05
Flowering time (FT; days)	257.30 ± 13.64	312.60 ± 12.70	257.30 ± 13.64	312.60 ± 12.70	170	229-326	t= −12.667 p < 0.05

Trait measurements

Four traits were measured, including two leaf morphological traits (leaf area, and leaf length-width ratio), one leaf physiology trait (SLA) and one phenology trait (flowering time). Plants for leaf trait measuring were consistent with our previous work [25], which include 10 P. depressa individuals, nine P. danxiaensis individuals, three F₁ hybrids, and 201 F₂ individuals. The genotypes of each individual were unique.

We scanned three mature and fully expanded leaves collected from each individual (Figure 1). We used the ImageJ v.1.49 [28] software to quantify leaf area, leaf length and leaf width. The largest mature leaf of each individual was sampled to measure the maximum leaf area (LA; cm²), length (LL; cm), and width (LW; cm). A derived trait, leaf length-width ratio (LWR) = LL LW⁻¹, was calculated. Following this measurement, the leaves were then placed in an air oven to dry for a minimum of 72 h at 80 °C, following which the final dry mass (MD; g) was recorded. The SLA (in units of cm²g⁻¹) of each leaf was calculated as: SLA = LA MD⁻¹. The flowering time was scored as the number of days between seedling transplantation and the day the first corolla was fully expanded. Traits data for analysis are available at the FIGSHARE repository: https://doi.org/10.6084/m9.figshare.7270715.v3.

We assessed whether parental trait values differed significantly using two-sided tests in the R statistical environment [29]. We calculated Spearman’s correlations between all pairwise combinations of traits in the F₂ population using R [29]. In the phenotypic correlation analysis, we assessed the four traits mentioned above together with 13 previously reported traits [25, 27]. Statistical significance was evaluated using the sequential Bonferroni correction to adjust significance levels for multiple comparisons [30].

QTL mapping

A linkage map containing 2,484 markers with an average marker distance of 0.96 cM was developed in our previous work [25]. We mapped the QTLs on the linkage map by using MapQTL v.6.0 software [31], as described by Feng et al. [27]. In brief, we determined the logarithm of the odds (LOD) score threshold by permutation testing. Multiple QTL Model (MQM) analysis was repeated to characterize loci with LOD scores higher than the threshold. The QTLs were finally estimated as 1.5-LOD support intervals for each MQM cofactor. The QTL maps were drawn using MapChart v.2.2 software [32].

Locus analysis

To evaluate the significance of the correspondence between QTLs in different traits, a previously described statistical test based on the hypergeometric probability distribution [33] was used to calculate the probability of obtaining the observed number of matching QTL by chance alone [34]. The equation is as presented in Wu and Davis (1993) [35]. In this analysis, we evaluated the QTLs underlying traits in this study as well as other traits that we studied previously [25, 27].

A QTL was considered as minor or major depending on the phenotypic variance it expressed. A QTL expressing a phenotypic variance of more than 15% was considered as a major QTL [36]. We further analyzed the pairwise interactions between QTLs for each trait with markers closer to QTL peaks using R/qtl [37].

Results

Analysis of phenotypic variance and correlations

We found significant phenotypic differentiation between all four measured traits (Table 1). The leaf size of the paternal lines (P. depressa) was larger than that of the maternal lines (P. danxiaensis) (Figures 1 and 2). The paternal individual leaves were long and narrow, and thus had much larger LWR values than that of the maternal individuals (Figures 1 and 2). The SLA of P. depressa was lower than that of P. danxiaensis (Figure 2), which may indicate that P. depressa has higher water-use efficiency. The F₁ plants showed intermediate phenotype between the two parental phenotypes, except for flowering time (Figure 2).

Figure 2. Histograms of the four analyzed traits in the F₂ population. Phenotypic means for P. depressa, P. danxiaensis and F₁ hybrids are indicated with blue, red and orange arrows, respectively.

The four traits were weakly correlated with each other (Figure 3). However, leaf area was positively significantly correlated with maximum leaf length and width (Supplemental Figure S1). Strong significant correlations were also observed among leaf physiological traits, with a negative correlation between SLA and leaf pigment concentrations (Supplemental Figure S1). These strong pairwise correlations may indicate that the trait divergence between P. depressa and P. danxiaensis is controlled by pleiotropic and/or tightly linked loci.

Figure 3. Correlations between and within the four traits (LA: maximum leaf area; LWR: leaf length-width ratio; SLA: specific leaf area; FT: flowering time) in the F₂ population. Scatterplots of all traits are given below the diagonal and Spearman’s correlation coefficients and associated significance are given above the diagonal.

QTL mapping analyses

A total of 40 QTLs were identified for the four phenotypic traits examined in this study (Figure 4), with the number of QTLs detected per trait ranging from seven for flowering time to 12 for either maximum leaf area or leaf length-width ratio (Supplemental Table S1). The phenotypic variance explained (PVE) for the identified QTLs ranged from 3.1% to 15.8%, with the majority (95%) of the QTLs having small effects (PVE < 10%; Supplemental Table S1). These QTLs accounted for 57.4–66.0% of the total proportion of the trait variation. The QTL for maximum leaf area explained the highest proportion of F₂ variation, whereas the QTL for SLA had the lowest explanatory power (Supplemental Table S1).

Figure 4. Genome location of significant QTLs for four traits. Confidence intervals of the QTLs are illustrated on the right side of the linkage groups.

Each QTL for three leaf traits (maximum leaf area, leaf length-width ratio, and SLA) explained their corresponding trait variance of less than 15% (Figure 5; Supplemental Table S1). This suggests that these leaf traits were controlled by some minor-effect loci. For flowering time, FT02 on LG02 is a major-effect locus with a LOD value of 9.44 and a variance explained of 15.8% (Figure 5; Supplemental Table S1). The P. depressa allele at FT02 was associated with earlier flowering by 17.89 d (Supplemental Table S1). These results suggest that the genetic architecture of flowering time includes a major-effect QTL and many minor-effect QTLs.

Figure 5. The phenotypic variance explained value of QTLs for four traits. LA: maximum leaf area; LWR: leaf length-width ratio; SLA: specific leaf area; FT: flowering time.

In these 40 QTLs, LA04 for the maximum leaf area and LWR09 for leaf length-width ratio on LG10 were overlapped (Figure 4). We calculated the significance of the correspondence between the overlapping QTLs and found that this QTL pair did not show significantly more overlap than expected by chance (p > 0.05) (Supplemental Table S2). This is consistent with the low phenotypic correlations among these traits. This indicated that these traits were controlled by a different genetic basis. Compared with previous QTLs for 12 flower and leaf traits and one reproductive isolation trait (hybrid pollen fertility) for the same population [25, 27], we found that both the maximum leaf area and leaf length-width ratio harbored some QTLs (LA01, LA05, LA07, LA12, and LWR03, LWR04, LWR10, and LWR11) that overlapped with those for other leaf size traits (maximum leaf length and maximum leaf width) (Supplemental Figure S2). Additionally, they exhibited significantly more overlap than expected by chance (Supplemental Table S2). These results were consistent with the strong correlations of the phenotypic values. However, the overlap between SLA09 for SLA and Cc11 for leaf carotenoid concentration on LG18 was more likely to be by chance. In addition, we did not find any QTL for flowering time that overlapped with other traits (Figure 4).

Table 2. Summary of significant epistatic interaction between QTLs.

Interacted loci ID	LOD	Var (%)	F value	p value	Significance level
LA02&LA10	2.677006	2.930819	2.700273	0.031979	*
LA03&LA04	2.552528	2.798692	2.570812	0.039315	*
LA06&LA07	2.983151	3.254096	3.020333	0.019114	*
LA08&LA10	2.408477	2.645293	2.421479	0.049819	*
LWR03&LWR05	3.208731	3.770538	3.257681	0.013011	*
LWR05&LWR06	2.585602	3.060950	2.605173	0.037222	*
SLA02&SLA08	2.440450	3.308638	2.545490	0.040930	*
SLA03&SLA09	2.757467	3.724322	2.887280	0.023689	*
SLA04&SLA08	2.502818	3.390668	2.612525	0.036788	*
FT01&FT06	3.025253	5.027619	3.217425	0.014197	*
FT01&FT07	2.649116	4.425249	2.802542	0.027583	*

Epistasis analysis among the identified QTLs was conducted with R/qtl. Only QTLs that had significant epistatic interactions are shown. Var (%), F value and p value represent percentage of variance explained, F statistics and F distribution, respectively.

*p < 0.05.

Loci interaction

There was evidence of epistasis for all four traits analyzed (Table 2). A significant pairwise epistatic interaction (p < 0.05) was found for 45% (18 of 40) of QTLs. Genetic epistasis could partly explain the phenotypic variation for these measured traits (Table 2 and Figure 6). The fraction of the phenotypic variation explained by epistasis for each trait ranged from 24.3% to 66.4% (Figure 6). There were seven loci for the maximum leaf area and three for the leaf length-width ratio involved in digenic interaction (Supplemental Figure S3). The effect of LA03 in the background of LA04 with the homozygous male parent genotype (AA) was evident from the fact that the maximum leaf area of all genotypes was the same. For LA08, the maximum leaf area of the heterozygote (AB) was high under the background of homozygous LA10 (BB) (Supplemental Figure S3). The phenotypic values were similar between LA10 heterozygotes (AB) and homozygotes (AA), which is consistent with the high degree of dominance of this locus. For SLA, five out of nine loci were involved in three pairs of digenic interactions: SLA02/SLA08, SLA03/SLA09 and SLA04/SLA08 (Supplemental Figure S3; Table 2). Genetic conflicts were obvious in all epistatic loci pairs; in other words, paired loci that inherit alleles from different parents will result in low phenotypic values. For flowering time, digenic interactions were identified between two QTL pairs involving three loci: FT01/FT06 and FT01/FT07 (Supplemental Figure S4; Table 2).

Figure 6. Proportion of phenotypic variation explained by all of the single QTLs and epistatic interactions. LA: maximum leaf area; LWR: leaf length-width ratio; SLA: specific leaf area; FT: flowering time.

Discussion

The genetic architecture of microgeographic adaptation

We detected 40 QTL for the four microhabitat adaptation traits, with each trait influenced by at least seven loci. These results indicate a highly polygenic genetic architecture involved in microhabitat adaptation in P. depressa and P. danxiaensis in sympatry. For leaf morphological and physiological traits, each QTL explained small to moderate effects (3.1%–9.5%) of the phenotypic variance. In our previous research, we found that differences in leaf size and pigment contents were also dominated by substantial QTLs with small to moderate effects in these two sympatric species [25].

We found four pairs of digenic interactions among leaf size QTLs, which indicated that the genetic architecture of leaf size was more complex than the other traits. Furthermore, this was consistent with the leaf size traits (leaf length and leaf width) we studied previously by our group [25]. The polygenic genetic basis of leaf shape (leaf length-width ratio) is in agreement with that observed in other plant species. For example, Zhang et al. (2015) identified seven QTLs, each of which explained up to 10% of the variance for the leaf length-width ratio in rice [38]. In poplar, nine leaf-shape traits were controlled by QTLs with small to moderate effects [39]. The small-to-moderate effect size of the leaf shape QTLs in these studies is also consistent with our research. However, the leaflet width-length ratio in tomato is controlled by two large-effect QTLs, each explaining more than 19% of the variance [40]. This suggests that the genetic architectures of leaf shape in plants are lineage-specific. SLA is an important indicator for water usage efficiency in plants [38, 41]. Although the ecological significance of SLA has received attention in many crops, there has been less research into the genetic architecture of SLA in natural systems. Our results indicated that the differences in SLA between P. depressa and P. danxiaensis were dominated by some QTLs with small effects, consistent with the observation in soybean [42] but in contrast to that in maize [43] and rice [44].

For flowering time, our analysis identified one major-effect QTL, which explained 15.8% of the trait variation. Six minor-effect individual QTLs jointly explained 4.9–10.7% of the flowering time variance. These findings are consistent with Ferris et al. (2017), who reported that few QTLs of large effect underlie interspecific divergence in the phenotypic traits of sympatric species of the Mimulus guttatus complex [45]. Similarly, Dittmar et al. (2016) identified two candidate genes associated with local adaptation, underlying a major flowering time QTL in Arabidopsis thaliana [46]. In addition, we tried to search the SNP locus underling the 1.5-LOD peak region of this major-effect flowering time QTL (FT02 on LG02) on the P. huaijiensis reference genome [7]. We found a compelling candidate gene copy, FRIGIDA-LIKE PROTEIN (gene number of Phu03062 located on Chr06 in the reference genome), at 0.9 Mb downstream of this SNP locus. FRIGIDA regulates the floral repressor FLC in Arabidopsis [47, 48]. This gene may be an important adaptive gene for the divergence of this pair of sister species (P. depressa and P. danxiaensis). However, this hypothesis should be further explored experimentally, as the reference genome we used was not of the parental species.

Reproductive asynchrony (e.g. different flowering times) is a prevalent isolation barrier between closely-related sympatric species [45, 49, 50]. In this study, P. depressa and P. danxiaensis had partially overlapping flowering times for about one month [25], and the mean flowering time of P. depressa was much earlier than that of P. danxiaensis (Figure 2). We identified a large-effect QTL underlying flowering time divergence in the present study. Theoretically, the local adaptation of closely related populations in sympatry should be due to few loci of large effect because of the homogenizing effects of gene flow [51]. In other words, large-effect loci may maintain the purity of species. These results indicated that flowering time may play an important role in the microhabitat adaptation and maintenance of species identity in the sympatric pair of P. depressa and P. danxiaensis.

We demonstrated that leaf traits and flowering phenology had contrasting patterns of genetic architecture and were relatively independently genetically controlled in this species pair. This genetic independence will result in independent selection on those traits. This is consistent with that found in Arabidopsis [52], but differs from the results in maize [21] and turnip [53]. Different patterns of genetic architecture also indicated that microhabitat adaptation in the sympatric Primulina was unequable for different traits, with leaf traits diverging via the accumulation of many small changes, and flowering phenology divergence being driven by major-effect variance.

Microhabitat adaptation and reproductive barriers in sympatric Primulina species

Sympatric sister species, particularly those connected by some degree of gene flow, provide an ideal system to investigate microhabitat adaptation and reproductive barriers. The overlapping flowering period provides an opportunity for the gene flow between P. depressa and P. danxiaensis in Danxia Mountain. This was supported by the phylogenetic results that show that the P. danxiaensis population sympatric with P. depressa shows greater genetic similarity to P. depressa than it is to conspecific populations [26].

The divergence of species with gene flow requires more barriers to accelerate or maintain speciation. In our previous work, we found that flower size significantly diverged between these two species and might be due to natural selection [25]. Hybrid sterility contributes to the post-zygotic barriers [27]. In this study, we found that flowering phenology may also play an important role in reproductive barriers. Both flowering phenology and hybrid fertility were controlled by similar genetic architecture patterns, with one major effect locus plus multiple small-to-medium effect loci. Two QTLs for hybrid fertility overlapped with QTLs for flower size; however, no overlapping QTL was identified between flowering phenology and other traits. This indicated that the genetic architecture of flowering phenology was independent of the other traits. Using the current data, it was not possible to distinguish if the divergence of flowering phenology was caused by independent selection or genetic drift.

Traits contributing to reproductive isolation include pre-zygotic barriers such as ecogeographic isolation [54–56], pollinator-mediated isolation [57–59], and post-zygotic barriers such as hybrid inviability and sterility [60–62]. In the sympatric sister species of Primulina, we previously identified the genetic architecture of pollinator-mediated isolation (flower size) [25] and post-zygotic isolation (hybrid sterility) [27]. In this study, both flowering phenology and leaf traits were associated with ecogeographic isolation. Our work improved the understanding of the genetic architecture of reproductive isolation and microhabitat adaptation in sympatric Primulina species.

Conclusions

Primulina species are cave-associated organisms that are an excellent model for research into microgeographic adaptation. In this study, we investigated the genetic architecture of flowering phenology and leaf morphology and physiology between a pair of sympatric sister species of Primulina. We found that flowering time was controlled by one major- and six minor-effect QTLs, while leaf traits were influenced by nine to 12 QTLs of small to moderate effect. The genetic architecture of flowering time and SLA was genetically independent of other traits. The different genetic architecture patterns of microhabitat adaptation traits also indicated that leaf traits diverged with the accumulation of many small changes, while the divergence of flowering phenology was driven by major-effect variance.

Authors’ contributions

CF designed this project, conducted the experiments and wrote the manuscript. CF, SZ and JZ analyzed the data. SZ and JZ revised the manuscript.

Acknowledgments

We are sincerely grateful to Prof. Ming Kang and Dr. Lihua Yang (South China Botanical Garden, CAS) for plant collection and advise for project designing.

Data availability statement

All the data and plant material are available from the first author (CF) upon reasonable request. Traits data for analysis are available at the FIGSHARE repository: Feng C, Feng C, Yang L, Kang M, Rausher M. qtlST script and phenotypic data. figshare. Dataset; 2019. https://doi.org/10.6084/m9.figshare.7270715.v3

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China (31900278), the Biological Resources Programme, Chinese Academy of Sciences (KFJ-BRP-007-013) and the Special Project for Lushan Plants (2021ZWZX18).