Data on the distribution of R. dauricum in seven forestry bureaus in the Greater Khingan Mountains was retrieved ( Supplementary Table S1 ). Plants from the following 13 locations were selected for phenotyping and molecular analysis ( Supplementary Figure S1 , Supplementary Table S2 ): Daling (DL), Hongwei (HW), Jiagedaqi (JGDQ), Tahe (TH), Huzhong (HZ), Taoshan (TS), Jinshantun (JST), Hongxing (HX), Meixi (MX), Youhao (YH), Mordauga (MEDG), Haolibao (HLB), and Delbuer (DEBE). They covered the Greater Khingan (DL, HW, JGDJ, TH, and HZ), the Lesser Khingan (TS, JST, HX, MX, and YH), and Inner Mongolia regions (MEDG, HLB, and DEBE) that were representative of the natural habitats of R. dauricum. For each population, 30 sampling blocks (5 m x 5 m for each) with an inter-block space of 30 m were set. Data were collected on individual plants of 10-12 years. Sampling was performed on May, 2022.

A digital caliper or metric tape was used to measure plant height, ground diameter, and various leaf/flower parameters. Three randomly selected branches were used to measure the ground diameter. Leaves from the current-year branches were measured according to the DUS guidelines developed for R. dauricum. Three flowers from each plant were selected to measure the flower diameter. Color determination was performed according to the standards of the Royal Horticultural Society. Specifically, scores 1-3 were assigned to plants with the following flower groups ( Figure 1 ): the red-purple (level I), purple (level II), and purple-violet (level III) groups, respectively. In addition, a score of 4 was assigned to white variants (level IV).

We used nine traits for our survey: flower color, flower amount, flower diameter, plant height, ground diameter, branch number, leaf length, leaf width, and leaf aspect ratio. Flower traits were selected because they are the key features determining the ornamental value of R. dauricum. Other traits were selected due to their importance to the overall plant shape and the commercial potential of R. dauricum.

Dormant branches were sampled in February, 2022 and then stored at room temperature to break dormancy. Flower buds were removed, and the branches were hydro-cultured. Newly emerged leaves (3 to 4 cm long) were used for DNA extraction. Thirty samples were collected for each of the 13 populations. DNA extraction was conducted using the CTAB method (Xiao et al., 2007) with 0.5 g of leaves per sample.

A total of 68 SSR markers were obtained from species closely related to R. dauricum (Yang et al., 2018; Zhang et al., 2016; Zhou et al., 2019). These primers were first tested in two plants randomly selected from the 13 populations. Out of these 68 primers (the primer information is shown in Supplementary Table S3 ), 13 showed bright and repeatable bands with polymorphism and were used for further analysis of the 390 samples (30 samples per population; 13 populations). Information on the available primer sequence, length, repeat motifs, and annealing temperature were summarized in Supplementary Table S4 . For PCR amplification, reactions were performed in 20 μL consisting of 5X PCR buffer, forward primer, reverse primer, DNA template, PCR Enhancer, probe, and H₂O. The reactions were carried out with the following steps: pre-denaturation at 95°C for 15 min, followed by 35 cycles of amplification (95°C: 30 s, 56°C: 30 s, and 72°C: 30 s); and final extension at 72°C for 3 min. After amplification, electrophoresis detection was performed using the ABI3130xl gene sequencer.

Morphological data was processed with Excel, Spass, and NTSYS (Rohlf, 2000). The degree of dispersion was measured using the coefficient of variation (CV). A two-factor nested design variance analysis was performed to determine the inter- and intra-population variations (Li et al., 2002). The unweighted paired arithmetic mean (UPGMA) method was used for clustering analysis based on the phenotypical data. The neighbor-joining (NJ) method with pair-wise Euclidean distance was used to construct a dendrogram. Associations between phenotypical data and ecological factors including longitude, latitude, altitude, temperature, and precipitation were determined by calculating the Pearson's Correlation coefficient.

The data matrix of the SSR experiment was used to calculate various parameters to infer genetic diversity using Popgene32 and Powermarkersoftware (Liu and Muse, 2005; Yeh et al., 1999). These parameters included number of alleles (Na), number of effective alleles (Ne), expected heterozygosity (He), observed heterozygosity (Ho), Shannon's information index (I), percentage of polymorphic sites (PPL), Nei's gene diversity index (H) (Masatoshi, 1973), gene flow (Nm), gene frequency, site polymorphism information content (PIC) value, genetic differentiation index (Fst). The Chi-square test was used to determine statistical significance. TASSEL was used to determine the correlation between genotype data and phenotypic traits.

We first surveyed the distribution of R. dauricum in seven forestry bureaus in the Greater Khingan Mountains ( Supplementary Table S1 ). A total area of 860,994 hectares of R. dauricum was recorded across the seven forestry bureaus. Wild populations of R. dauricum are distributed within 122°18'28''E to125°18'37''E and 50°12'21''N to 53°24'17''N, which spans the middle temperate zone and cold temperate zone. R. dauricum can grow in multiple habitats including mining sites, volcanic rocks, sunny slopes, forests (Betula platyphylla and Larix gmelinii), and riverbanks ( Figure 2 ). The altitude range is 200 m to 1300 m, although 400 m to 800 m is the preferred altitude. R. dauricum thrives in full sun but can tolerate shade. It was usually found in patches with very few individual plants scattered.

A great degree of variation in flowers was observed in the natural populations of R. dauricum. First, the corolla was usually shallowly lobed, but deeply or completely lobed corolla was also observed. Five-lobed corollas were the most common ones, while six- or seven-lobed were rare. The flower diameter was 2.5 to 3.5 cm with some reaching up to 4 cm. The number of flowers on a single branch ranged from 1-4. Based on the geographic locations, R. dauricum can bloom between late April and early May (Early), from May 1 to May 10 (Middle), and after middle May (Late). In addition, a secondary bloom may occur after the end of September (Secondary).

The flower color varied from white, light purple, pink-purple, purple-red, and purple ( Figure 1 ). Interestingly, we found a white variety, which was distinct from R. dauricum in terms of corolla color, anther color, filament color, branch color, spot color, and leaf color. This variety had been named R. dauricum var. Album (Wilson, 1923). In addition, we found yellow-green spots on some of the white-flowered R. dauricum, which had been used to develop 'Ao Xue' via tissue culture and propagation (Wang et al., 2022). Another variety we found showed spherical inflorescences with 8-14 flowers on the same branch, which were significantly more than that in a typical R. dauricum. These plants were used to develop a new variety 'Yanricai' (Wang et al., 2023), which had been registered to the Royal Horticultural Society. The numerous flowers make it a variety of high ornamental value.

To determine the phenotypic trait diversity of the 13 populations, we measured nine traits and summarized their mean values ( Table 1 ) and coefficient of variations (CVs, Table 2 ). Based on the inter-population variation, the nine traits were ranked in the order of flower color > flower number > ground diameter > leaf width > branch number > leaf length > leaf aspect ratio > flower diameter > plant height. Of particular interest to us were flower color and number, which showed CVs of 0.2651 and 0.1974, respectively, indicating a large degree of heterogeneity among the populations. It is noteworthy that variations existed for traits showing a small degree of CVs. For instance, the plant height ranged from 158.64 cm (HZ) to 185 cm (MEDG) despite a low inter-population CV (0.0688). We also found that the overall inter-population CV was population-dependent, with HZ and HX showing the largest degree of variations (CVs of 0.1838 and 0.1799 respectively).

To further determine the inter- and intra-population variations, we performed the nested design variance analysis and found that all the nine traits showed significant variations at both the inter- and intra-population levels (α = 0.01; Table 2 ). At the inter-population level, the F statistics ranked the nine traits in the order of plant height > flower number > branch number > leaf length > leaf width > flower diameter > flower color > ground diameter > leaf aspect ratio. At the intra-population level, the order of F values from large to small was plant height > flower color > flower number > number of branches > leaf length > flower diameter > ground diameter > leaf width > leaf aspect ratio.

Correlations between the nine traits and various ecological factors of the sampling site were summarized in Table 3 . First, several traits including plant height, ground diameter, branch number, and flower number were negatively correlated to longitude, while flower diameter showed a positive correlation with longitude. Second, plant height, ground diameter, branch number, flower number, and leaf length showed significant positive correlations with the latitude, confirming that R. dauricum prefers habits of high latitudes. Third, plant height, branch number, and flower number were extremely significantly positively correlated with the altitude, whilst leaf length and flower diameter decreased as the altitude increased. Fourth, the annual average temperature showed a positive impact on several traits (flower diameter and leaf length) and a negative one on others (plant height, ground diameter, branch number, flower color, and flower number). Lastly, the annual precipitation showed a positive correlation with flower diameter and negative correlations with plant height, ground diameter, and branch number.

The greater the genetic distance, the smaller the genetic consistency, which proves that the farther the genetic relationship between the two; on the contrary, when the genetic distance is smaller, the greater the genetic consistency, the higher the genetic relationship. It can be seen from the Table 4 that the genetic distance between HX and DL was 0.0417, the genetic identity was 0.9592, and the genetic relationship between them was the closest. The genetic distance and genetic identity between JGDQ and MEDG were 0.4463 and 0.64, and the genetic relationship between them was the farthest.

The UPGMA dendrogram constructed by similarity coefficient data obviously showed that three clusters were formed when the similarity coefficient was 1.30 ( Figure 3 ). The resulting dendrogram revealed three major groups. The first group consisted of MEDG, HLB, HZ, DL, and HX. The second group comprised MX, YH, TS, and JST. The third group included JGDQ, DEBE, HW, and TH. The group 1 accessions were mainly R. dauricum collected from the western Greater Khingan Mountains, HX accessions were also assigned into this groups. All four Lesser Khingan Mountains accessions were classed to group 2. The R. dauricum genotypes of group 3 were primarily collected from the east of the Greater Khingan Mountains. While most populations clustered according to geographical distance, there were exceptions. This indicated a discontinuity in the phenotypic trait variation of the 13 populations.

Cross-species microsatellite markers are powerful for studying the genetic diversity of species whose genome information is unavailable. We randomly selected 68 pairs of published SSR primers from the same genus and found that the majority resulted in positive PCR bands (41 out of 68). This supports the conserveness in the flanking sequences of SSR among closely related species (Guo, 2007; Nie, 2013; Wang et al., 2010a, 2016, 2010b). Based on the band clarity, we selected 13 SSR markers to further evaluate the genetic diversity of 390 individuals from the 13 populations ( Table 5 ). Four markers including Rd-SSR-43, Rd-SSR-60, Rd-SSR-9, and Rd-SSR-65 resulted in a single band, indicating monomorphism of the locus. Nine pairs of SSR primers (Rd-SSR-53, Rd-SSR-34, Rd-SSR-51, Rd-SSR-38, Rd-SSR-59, Rd-SSR-57, Rd-SSR-16, Rd-SSR-66, and Rd-SSR-19) showed multiple bands and thus polymorphism among the 13 populations.

The nine pairs of polymorphic primers covered 25 alleles, with 1-7 alleles per marker. The number of effective alleles ranged from 1 to 4.0563 (mean of 1.8662). The Shannon information index I ranged from 0 to 1.6410 (mean of 0.4752). Polymorphism information index (PIC) ranged from 0.0712 to 0.7605, identifying Rd-SSR-57 and Rd-SSR-66 as the preferred markers for genetic diversity studies in R. dauricum. The observable heterozygosity (Ho) ranged from 0 to 0.5385 (mean = 0.2514) and the expected heterozygosity (He) ranged from 0 to 0.7862 (mean = 0.3722), suggesting the presence of inbreeding and genetic variation among the populations. This was further confirmed by the average heterozygosity (Ave_Het), which ranged from 0 to 0.2692. In addition, the genetic diversity index Nei ranged from 0 to 0.7535 (mean = 0.3575) and indicated that the genetic diversity varied greatly among the alleles.

We also determined the genetic diversity among the 13 R. dauricum populations ( Table 6 ). The average number of alleles is about 15 with HZ and TH populations showing 18 alleles. There was a general agreement between the observed heterozygosity and the expected heterozygosity, indicating random breeding of these populations. JGDQ and DL showed the lowest observed heterozygosity, and HZ and TH showed the largest. Accordingly, HZ and TH showed the highest Nei's gene diversity index values (0.1923 for both).

Pearson correlation analysis of genetic and geographic distances was performed using Origin. The results are shown in Figure 4A (r=0.1112) and Figure 4B (r=0.0499), which shows that there is a certain correlation between genetic distance and geographic distance between populations, but the correlation is not significant. This indicates that genetic variation has little correlation with geographic distance. This may result in lower gene exchange between populations, leading to genetic differentiation between populations.

Next, we used Nei's F statistic to analyze the SSR data and found that the Fst values ranged between 0.35 and 1 ( Table 7 ), with an average inter-population Fst of 0.6556. Thus, the within-population variance (34.44%) was lower than the between-population variance (65.56%). Clustering analysis revealed three groups in the phylogenetic tree ( Figure 5 , Supplementary Table S5 ) with a clear separation among HLB, MEDG, and the remaining 11 populations. The relatively large genetic distance between HLB/MEDG and other populations can also be reflected in the principal coordinate plot ( Figure 6 ). Overall, we observed an increase in genetic distance among the populations as their geographical distance increased, consistent with the geographical distance model (Jiang et al., 2016). Similar genetic backgrounds were found in relatively close populations, for example, between HW and TH, HLB and MEDG, JST and TS, and MX and YH. However, deviation from this model was also common, such as HZ clustered closely with MX although they were geographically distant. Many factors including human activities and insect pollination may be behind this, which requires further investigation.

We next used GLM to determine the association between SSRs and phenotypic traits and found that 3 loci with statistical significance (p < 0.05, Table 8 ): Rd-SSR-57 for leaf aspect ratio, ground diameter, branch number, and flower number; Rd-SSR-66 for flower color and plant height; and Rd- SSR-34 for plant height and flower diameter ( Supplementary Figure S2 ). The explanation rate of phenotypic variation by individual SSRs ranged from 0.3179 to 0.7110 ( Table 8 ). Among them, Rd-SSR-66 showed an explanation rate of 0.7110 for flower color, while both Rd-SSR-66 and Rd-SSR-57 were significantly correlated with flower color and leaf aspect ratio (Q < 0.1). Further analysis revealed a significant variance in flower based on Rd-SSR-66 genotypes, with genotypes 244:244, 230:230 and 227:230 showing darker flower color in TH, HW, DEBE, YH and 227:227 and 227:236 corresponding to lighter flower color values in MEDG and HLB ( Figure 7A ). For Rd-SSR-57, a significant difference was found in leaf length-to-width ratio between TS/JST and the remaining populations (P =0.0008, Figure 7B ).