Introduction
Salvia miltiorrhiza Bunge of the family Lamiaceae, is an outcrossing insect-pollinating perennial herb native to China. The red root of S. miltiorrhiza is a traditional Chinese medicine and play an important role in the treatment and recovery of patients with hypertension and coronary heart disease. Moreover, it has recently begun to be used for the clinical treatment of organ damages caused by COVID-19 (Zhang et al., Reference Zhang, Li, Li, Wong, Wang and Zhang2021). Therefore, the demand for S. miltiorrhiza is increasing remarkably. As a medicinal plant, it has been cultivated in China, but it is generally believed that the content of major secondary metabolites from wild herb is higher than those from cultivated herb (Gao et al., Reference Gao, Jia, Gao, Wang and Xiao2005; Zhang et al., Reference Zhang, Yu, Yang, Qi, Liu, Deng, Cai, Li, Sun and Liang2018). The natural resources of S. miltiorrhiza in central China have rapidly declined in recent years because of the over-mining of its roots.
The Qinling Mountains in central China is one of the world's major centres for plant diversity. The well-preserved glacial relics on Taibai Mountain (3767 m) shows that glacier occurred in the region during the Quaternary (Rost, Reference Rost2000; Zhang et al., Reference Zhang, Liu, Chen, Liu, Harbor, Cui, Liu, Liu and Zhao2016). As suggested by previous studies, the Qinling Mountains may have acted as glacial refugia for many species including walnut and peonies during glacial period due to their heterogeneous topography (Feng et al., Reference Feng, Zhou, Zulfiqar, Luo, Hu, Feng, Malvolti, Woeste and Zhao2018; Xu et al., Reference Xu, Cheng, Peng, Sun, Hu, Li, Xian, Jia, Abbott and Mao2019). The Dabie Mts are composed of a chain of ancient and isolated massifs. To date, there was no evidence of glacier in the middle-low mountainous areas of eastern China, but temperatures were estimated to have been much lower during glacial period than at present (Wang and Sun, Reference Wang and Sun1994). The Dabie Mts are often postulated as refugia for warm-temperate forest species (Liu et al., Reference Liu, Hao, Liu, Wei, Cun and Wang2014; Yang et al., Reference Yang, Dick, Yao and Huang2016), and several new species have been identified recently only in the area (Sun et al., Reference Sun, Orozco-terWengel, Chen, Sun, Sun, Wang, Shi and Zhang2021). Today, natural populations of S. miltiorrhiza can be found in restricted areas such as China's Qinling Mts in Henan and Shannxi, the Mengshan Mts in Shandong, the Dabie Mts in Anhui and other places. Most of the previous studies reported the population genetic diversity of cultivated S. miltiorrhiza (Song et al., Reference Song, Li, Wang and Wang2010; Xu et al., Reference Xu, Liu, Huang, Wang, Zhang, Liu, Liao, Yuan and Zhou2013; Zhang et al., Reference Zhang, Li and Wang2013; Peng et al., Reference Peng, Ru, Wang, Wang, Li, Yu and Liang2014) and wild populations located in Anhui and Hubei using 2b-RAD genome-wide genotyping method (Zhou et al., Reference Zhou, Zhang, Huang, Xiao, Wu, Qi and Wei2021). However, no effort has been made to assess the genetic variability of wild S. miltiorrhiza using microsatellite markers.
To better understand and make use of this medicinal plant resources, we performed a sampling of S. miltiorrhiza naturally distributed in the mountainous regions in central eastern China. The objectives of this study were to evaluate the genetic variation of wild S. miltiorrhiza using microsatellite markers, explain the genetic structure based on historical climate and develop the conservation strategies.
Materials and methods
Sampling
A total of 215 specimens of S. miltiorrhiza were collected from 11 natural locations in central eastern China, wild populations of S. miltiorrhiza usually grows on the slopes and understories in broad-leaved deciduous forests. Of these locations, HXX, HLS, NNZ, SYX, SZS and SZA were located in the Qinling Mts of Henan and Shaanxi Provinces, SMS was located in the Mengshan Mts of Shandong Province, and AHS, AHX, ALQ and ALQ2 were located in the Dabie Mts and its surrounding areas of Anhui Province (Table 1, online Supplementary Fig. S1). Species morphological characters such as leaf and root were investigated carefully. Samples of fresh leaves were dried using silica gel and maintained at room temperature.
Table 1. Sampling locations of wild Salvia miltiorrhiza populations in central eastern China
DNA extraction and microsatellite genotyping
Genomic DNA was extracted from leaves using Plant Genomic DNA Extraction Kit (Tiangen, Beijing) according to the manufacturer's protocol and stored at −20°C. Twelve microsatellite markers (online Supplementary Table S1) were selected for this analysis: SoUZ012, SoUZ014, SoUZ020 (Radosavljević et al., Reference Radosavljević, Jakse, Javornik, Satovic and Liber2011), SoUZ024 (Radosavljević et al., Reference Radosavljević, Satovic, Jakse, Javornik, Greguraš, Jug-Dujaković and Liber2012), ssr31, ssr34, ssr37, ssr40, ssr43, ssr65, ssr66, and ssr68 (Xu et al., Reference Xu, Liu, Huang, Wang, Zhang, Liu, Liao, Yuan and Zhou2013), and then labelled with the fluorescent dyes FAM, HEX, or TAMARA on the 5′ end of each forward primer. PCR were performed in 20 μl of reaction containing 30–50 ng genomic DNA, 2 μl 10 × PCR buffer, 0.2 mM dNTPs, 0.2 μM of each primer, 0.5U Takara Taq DNA polymerase (Takara Biotechnology (Dalian) Co., Ltd.).
PCR amplifications were conducted under the following conditions: 94°C (5 min), and then 30 cycles at 94°C (30 s)/annealing temperature (online Supplementary Table S1) (30 s)/72°C (45 s), and a final extension step at 72°C (8 min). PCR products were separated and visualized using agarose gel stained with Gelred, and then the ABI 3730XL DNA sequence analyser and GeneMaker were used to determine the allele sizes.
Genetic diversity
The null alleles and genotyping stutter were tested using Microchecker version 2.2.3 (Van Oosterhout et al., Reference Van Oosterhout, Hutchinson, Wills and Shipley2004). FreeNA was used to estimate null allele frequencies for each locus with 500 bootstrap replicates (Chapuis and Estoup, Reference Chapuis and Estoup2007). The number of alleles (Na) and effective number of alleles (Ne) were calculated in Popgene version 1.31 (Yeh and Boyle, Reference Yeh and Boyle1999). Popgene was also used to calculate Shannon's information index (I), observed heterozygosity (H O) and expected heterozygosity (H E) of populations. Nei's gene diversity (H), the allelic richness (A R) and Wright's inbreeding coefficient (F IS) were calculated by FSTAT version 2.9.3 (Goudet, Reference Goudet1995).
Genetic differentiation and structure
The fixation index (F ST) was used as measures of population differentiation. Pairwise F ST values between populations and the effective number of migrants (Nm) were estimated in Arlequin with 1000 permutations (Excoffier and Lischer, Reference Excoffier and Lischer2010) and the correlation between F ST and geographical distance was examined by Mantel test in the ‘ecodist’ R package (Goslee and Urban, Reference Goslee and Urban2007).
Genetic variations within and among populations were calculated by analysis of molecular variance (AMOVA) in Arlequin. Population structure was analysed using STRUCTURE version 2.3.4 (Evanno et al., Reference Evanno, Regnaut and Goudet2005) based on an admixture model with K varying from 2 to 11, ten independent runs were performed for each value of K with a burn-in of 10,000, followed by 100,000 Markov chain Monte Carlo (MCMC) iterations, and then bar plot of assignments was generated using DISTRUCT version 1.1 (Rosenberg, Reference Rosenberg2004). The delta K was applied to the microsatellite data in order to infer the optimal K value. Population structure was also examined by the Discriminant Analysis of Principal Components (DAPC) using the ‘adegenet' package implemented in R (Jombart, Reference Jombart2008).
Species distribution modelling
To better understand the genetic structure of the wild populations of S. miltiorrhiza, the potential present and Last Glacial Maximum (LGM) distribution was modelled by MAXENT ver. 3.3.3k (Phillips et al., Reference Phillips, Anderson and Schapire2006). The present and LGM bioclimatic variables were obtained from Worldclim database at a resolution of 30 arc-seconds and 2.5 arc-minutes, respectively. The suitability scores were classified into four intervals from low (0.0–0.25), moderate (0.26–0.50), high (0.51–0.60) to very high (>0.60) (Liu et al., Reference Liu, Berry, Dawson and Pearson2005), and then, the models were visualized in ArcGIS ver. 9.3. (Esri, Redlands, CA, USA).
Results
Genetic diversity
Genetic diversity of 215 individuals from 11 locations was examined. Microchecker revealed that one locus (ssr65) showed evidence for a null allele with a frequency ⩽ 0.2 by expectation and maximization (EM) method. Thus, all loci were used for genetic analysis. The number of alleles per locus (Na) ranged from 2 (SoUZ012) to 13 (SoUZ014, ssr31 and ssr65) among all sites with an average of 8.083, and the effective number of alleles per locus (Ne) were 1.02 (SoUZ020) to 8.40 (SoUZ014). The main indexes describing genetic diversity are shown in Table 2. In the Qinling Mts, the highest level of allelic richness was recorded in HLS (4.889). In the Dabie Mts, the allelic richness was highest for population AHS (3.999), while it was the lowest for ALQ (3.034). Nei's gene diversity (H) ranged from 0.368 in ALQ to 0.588 in HLS with an average of 0.493. The average Shannon's information index was 0.947. In HLS, SZA, SYX, AHS and ALQ populations, the observed heterozygosity values were slightly higher than the expected heterozygosity, suggesting an excess in heterozygosity within populations and a possibility of outcrossing.
Table 2. Genetic diversity indexes for wild populations of Salvia miltiorrhiza by 12 microsatellite loci
N, number of individuals; Na, average number of observed alleles; H, Nei's gene diversity; A R, allelic richness; I, Shannon's Information index; H O, observed heterozygosity; H E, expected heterozygosity; F IS, inbreeding coefficient.
Genetic differentiation and structure
The AMOVA analysis revealed that as much as 81.99% of the total genetic diversity was attributable to differences among individuals within populations, which is characteristic for outcrossing and long-lived plants, and only small but significant percentage of the total genetic variation was attributed to genetic differences between the two clusters (F CT = 0.073; P < 0.01) (Table 3).
Table 3. Analysis of molecular variance (AMOVA) of wild Salvia miltiorrhiza in central eastern China based on microsatellites
Pairwise estimates of F ST among the wild populations of S. miltiorrhiza varied from 0.036 (between ALQ and ALQ2) to 0.312 (between ALQ and SYX) with an average of 0.145, which demonstrated that genetic differentiation existed among the populations (Fig. 1a). In contrast, gene flow (Nm) estimated among the populations varied from 0.552 (between ALQ and SYX) to 6.588 (between ALQ and ALQ2), the overall gene flow value obtained in this study was 1.796.
Figure 1. Pairwise F ST (a), mantel test between genetic distance and geographical distance (in kilometers) (b) and genetic similarity of wild Salvia miltiorrhiza populations in central eastern China inferred by discriminant analysis of principal components (DAPC) (c). The scatter plot was created to assist with visualizing the first two principal components, points represent individuals, eigenvalues represent the amount of genetic variation captured by the discriminant factors plotted as the x- and y-axes.
A mantel test of F ST against geographical distance (in kilometers) (Fig. 1b) showed a significant positive linear relationship between two variables among populations of S. miltiorrhiza (mantel r = 0.584, P < 0.01) and the F ST significantly increased with increasing geographic distance between populations, suggesting that part of the variance could be explained by isolation by distance.
Results of the STRUCTURE analysis (online Supplementary Fig. S1) showed that for S. miltiorrhiza, the highest delta K was obtained when k = 2, which is generally in accordance with the geographical position of the collecting sites. Among these populations, AHS, ALQ and ALQ2 from the Dabie Mts formed an independent cluster, whereas AHX and populations from the Qinling and Mengshan Mts were grouped in another. Moreover, the results of the DAPC based on the first two principal components provided further evidence for the validity of the two lineage groups (Fig. 1c).
Modelling of species distributions
The area under receiver operating characteristic curve values over 0.9 suggested the good performance of the models. The model for the present was congruent with the natural distribution S. miltiorrhiza in China (Fig. 2a); the generated model for the potential species distribution during the LGM implied the possibility of existence of multiple refugia in central eastern China (Fig. 2b), which might have led to the genetic difference among the wild populations of S. miltiorrhiza.
Figure 2. Species distribution modelling of wild Salvia miltiorrhiza populations at present (a) and during the Last Glacial Maximum (b). Colours indicate habitat suitability, from low (blue) to high (red), as calculated by the model.
Discussion
Spatial genetic structure
The microsatellite marker analysis of 215 samples of wild medicinal plant S. miltiorrhiza in central eastern China revealed a significant genetic structure and subsequent analysis showed that genetic differentiation between populations was most likely because of historical climate change, wild S. miltiorrhiza populations have been characterized by division into two subclusters corresponding to separate refugia. Distribution modelling showed the higher suitability for S. miltiorrhiza in these two regions for LGM. This was also supported by the significant correlation between geographic distance and genetic variation, and this type of pattern is often associated with postglacial colonization (Rešetnik et al., Reference Rešetnik, Baričevič, Batîr Rusu, Carović-Stanko, Chatzopoulou, Dajić-Stevanović, Gonceariuc, Grdiša, Greguraš, Ibraliu, Jug-Dujaković, Krasniqi, Liber, Murtić, Pećanac, Radosavljević, Stefkov, Stešević, Šoštarić and Šatović2016).
Additionally, subsequent gene flow between wild S. miltiorrhiza populations might have been influenced by the dispersal of seed and pollen. The mericarps of Salvia are small and smooth enough to be transported by wind, water, animals or gravity, but unlikely to achieve long-distance dispersal (Zona, Reference Zona2017), making gene flow primarily dependent on pollen dispersal between populations. It is well known that the genus Salvia is characterized by the famous staminal lever mechanism of the flower (Clasenbockhoff et al., Reference Clasenbockhoff, Speck, Tweraser, Wester, Thimm and Reith2004; Hu et al., Reference Hu, Takano, Drew, Liu, Soltis, Soltis, Peng and Xiang2018); pollen mediated gene flow in the presence of pollinators, therefore, can be expected, and the efficiency of pollinator-mediated gene flow depends on the flying ability of pollinators such as bee, bird and fly (Ohashi, Reference Ohashi2002; Benitez-Vieyra et al., Reference Benitez-Vieyra, Fornoni, Perez-Alquicira, Boege and Dominguez2014, Reference Benitez-Vieyra, Perez-Alquicira, Sazatornil, Dominguez, Boege, Perez-Ishiwara and Fornoni2019; Celep et al., Reference Celep, Atalay, Dikmen, Dogan and Classen-Bockhoff2014, Reference Celep, Atalay, Dikmen, Doǧan, Sytsma and Claßen-Bockhoff2020; Fragoso-Martínez et al., Reference Fragoso-Martínez, Martínez-Gordillo, Salazar, Sazatornil, Jenks, García Peña, Barrera-Aveleida, Benitez-Vieyra, Magallón, Cornejo-Tenorio and Granados Mendoza2018; Cairampoma et al., Reference Cairampoma, Tello and Claßen-Bockhoff2020).
In China, bees such as Xylocopa appendiculata and Bombus trifasciatus are excellent pollinators of S. miltiorrhiza, as these species are long-lived and have a wide geographical range (Gerling et al., Reference Gerling, Velthuis and Hefetz1989; Huang et al., Reference Huang, Liu and Huang2015; Liu et al., Reference Liu, Wang, He, Ning and Guo2015). Previous studies recorded a foraging distance of 1200 m for Xylocopa violacea (Gathmann and Tscharntke, Reference Gathmann and Tscharntke2002), and the distances travelled by Xylocopa flavorufa ranged from a minimum of 50 m to a maximum of 6040 m, with a median of 720 m (Pasquet et al., Reference Pasquet, Peltier, Hufford, Oudin, Saulnier, Paul, Knudsen, Herren and Gepts2008), which may contribute to higher gene flow and lower differentiation between neighbouring populations. However, in this research, genetic differentiation between the samples collected at different mountains was still detected.
The result was generally in accordance with earlier population genetic studies of Salvia. Based on enzyme electrophoresis method, researchers reported that genetic differentiation (G ST) in 14 populations of Salvia pratensis was 0.156 (Van Treuren et al., Reference Van Treuren, Bijlsma, Van Delden and Ouborg1991). In plastid DNA phylogeographic study of eight Balkan populations, Stojanović et al. (Reference Stojanović, Aleksić, Jančić and Jančić2015) detected the existence of two lineages in Salvia officinalis (F ST ranging from 0 to 0.934). Similar sub-clusters have been reported by Rešetnik et al. (Reference Rešetnik, Baričevič, Batîr Rusu, Carović-Stanko, Chatzopoulou, Dajić-Stevanović, Gonceariuc, Grdiša, Greguraš, Ibraliu, Jug-Dujaković, Krasniqi, Liber, Murtić, Pećanac, Radosavljević, Stefkov, Stešević, Šoštarić and Šatović2016) in S. officinalis populations distributed across the Balkan Peninsula using microsatellite markers, the global F ST value was 0.153, as well as Salvia rosemarinus (Zigene et al., Reference Zigene, Asfaw and Bitima2021). Amplified fragment length polymorphism analysis by Jug-Dujaković et al. (Reference Jug-Dujaković, Ninčević, Liber, Grdiša and Šatović2020) showed that S. officinalis from 25 localities distributed along eastern Adriatic was divided into two distinct ancestral clusters. Moreover, analysis of inter-simple sequence repeats (ISSR) also showed evidence for the existence of genetic differentiation among Salvia multicaulis populations (Talebi et al., Reference Talebi, Askary, Khalili, Matsyura, Ghorbanpour and Kariman2021).
Genetic diversity and conservation implications
The genetic diversity parameters analysed for the wild S. miltiorrhiza populations reflected the presence of higher variability (mean A R and H was 3.891 and 0.493, respectively) than those reported for cultivated S. miltiorrhiza (mean H: 0.159, Song et al., Reference Song, Li, Wang and Wang2010; mean H: 0.156, Xu et al., Reference Xu, Liu, Huang, Wang, Zhang, Liu, Liao, Yuan and Zhou2013; mean H: 0.279, Peng et al., Reference Peng, Ru, Wang, Wang, Li, Yu and Liang2014), cultivated S. officinalis (mean H: 0.195, Sarrou et al., Reference Sarrou, Ganopoulos, Xanthopoulou, Masuero, Martens, Madesis, Mavromatis and Chatzopoulou2017) and Salvia lachnostachys (mean H: 0.250, Erbano et al., Reference Erbano, Schuhli and Santos2015), these findings are typical for cultivated species as the cultivation always results in loss of genetic diversity because of artificial selection and founder effects, but lower than indigenous S. officinalis across the Balkan Peninsula (mean A R was 7.920, Rešetnik et al., Reference Rešetnik, Baričevič, Batîr Rusu, Carović-Stanko, Chatzopoulou, Dajić-Stevanović, Gonceariuc, Grdiša, Greguraš, Ibraliu, Jug-Dujaković, Krasniqi, Liber, Murtić, Pećanac, Radosavljević, Stefkov, Stešević, Šoštarić and Šatović2016). This is probably related to the over-mining of its roots.
To date, over collection have led to the decline and regional extinction of natural populations of this medicinal herb. Although the wild S. miltiorrhiza has not be listed as an endangered species, it is essential to conserve the most diverse wild populations for future breeding programmes. In situ conservation focused on these two units are recommended, and the HLS population in the Qinling Mts and AHS and ALQ2 in the Dabie Mts should be given priority for conservation due to their higher genetic diversity.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S1479262124000340
Acknowledgements
We are thankful to our 129 group members for their assistance with performing experiments. This work was financially supported by the Science and Technology Key Projects of Henan province (232102110061, 192102310176), Sanmenxia Programs of Soft Science (2022003031) and the National Key R&D Program of China during the 13th Five-Year Plan Period (2018YFD0600204-04).
Competing interests
None.
