Phenotype correlation analysis and excellent germplasm screening of herb Bletilla Rchb. f. based on comprehensive evaluation from thirty-three geographic populations

Background The Bletilla genus of Orchidaceae includes plants with great economic value, among which B. striata is the main traditional medicinal plant, and its pseudobulb, known as BaiJi, was first recorded in Shennong’s Classic of Materia Medica. However, there has been little systemic evaluation of the germplasm quality of Bletilla plants in China. In order to comprehensive evaluate the Bletilla resources in China and screen out the candidate phenotypic traits determining yield and/or quality of Bletilla, the variation of phenotypic indicators (pseudobulb, leaf, stem, inflorescence, flower) and active ingredients contents (polysaccharide, total phenolics and militarine) in different populations of B. striata and B. ochracea were investigated through 4 years’ common-garden experiment. Results There were abundant phenotypic variations and significant differences among different populations in the morphological phenotypes, pseudobulb weight and main active ingredient contents. AHBZ, HBLT and HBSN populations showed good prospects for industrial development, presenting higher quality in terms of yield and main active ingredient content. Pseudobulb yield, polysaccharide and total phenol content are positively correlated with phenotypic traits. Militarine content is negatively correlated with almost all indexes. Plant height, leaf width and stem diameter may be important indicators of potential excellent germplasms. Conclusions Bletilla is not strictly geoauthentic medicinal plants. B. ochracea could be accepted as an alternative resource to B. striata. The best harvest period of Bletilla is the third year after cultivation. Plant height, leaf width and stem diameter may be important indicators of potential excellent germplasms. These results provide important information required for the efficient screening and utilization of Bletilla germplasm resources. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-022-03540-w.

Common garden experiments retaining the same cultural environment can be used to analyze phenotypic variation among different populations to eliminate phenotypic variation induced by environmental factors and have become an important means of resource evaluation [20][21][22]. Through the collection of wild resources for common garden cultivation, quality evaluation was performed to determine the high-quality resources in the tested varieties and has been widely used in many medicinal plants, such as Epimedium sagittatum and Lycium chinense [22,23]. Gray relation analysis (GRA) is a multifactor statistical analysis method that includes dimensionless processing, correlation coefficient calculation, correlation degree and relative correlation degree calculation and is widely used in the comprehensive evaluation of the quality of medicinal materials [24]. The correlation is used in the context of a linear relationship between 2 continuous variables and the Pearson correlation coefficient is typically used for jointly normally distributed data [25]. Correlation analysis can be used to measure the relationship among yield, quality and agronomic traits to determine which trait affects yield and/or quality [26,27].
Systematic collection of wild and cultivated resources, investigation of phenotypic traits and active ingredient content, and selection of excellent germplasm are an important part of the effective utilization of medicinal plants [3]. Small-scale collection and resource evaluation of Bletilla have been performed. For example, Chinese researchers evaluated the germplasm resources of B. formosana in Yunnan Province, determined the total phenol content of B. striata in northwest of Hubei Province, and determined the polysaccharide content of B. striata from different areas in Anhui and Guizhou Provinces, while these studies have focused on phenotypic traits or multiple active components in a small region, or single component analysis in multiple provinces, lacking a large-scale collection and comprehensive evaluation of Bletilla resources in China [28]. Therefore, we selected 33 Bletilla populations, which covered their major distribution areas in China, to systematically evaluate the main phenotypic traits, pseudobulb yields and main active ingredients among different species and populations. We aimed to screen higher-quality germplasm resources of Bletilla and the key indicators related to economic characteristics, which will provide important information for the efficient breeding of B. striata. Our study will enhance our ongoing efforts in the sustainable development of traditional medicinal plants.

Morphological measurement and resource evaluation of pseudobulb weight
The morphological phenotypes of the aboveground parts of 33 populations (including 28 B. striata populations and 5 B. ochracea populations) were measured for three years ( Fig. 1; Table S2). Statistical analysis showed that there were abundant phenotypic variations and significant differences among different populations, especially the plant height, leaf blade length, stem diameter and inflorescence height with high F values and extremely low P values ( Table 1). The plant height, leaf blade length, inflorescence height and length of B. ochracea were greater than those of B. striata. In addition to leaf blade length in B. striata, the plant height, leaf blade length and width of B. striata and B. ochracea increased with increasing planting years, while the fruit number decreased in B. ochracea and B. striata (Fig. 2a).
There were significant differences in pseudobulb weight among different species and different populations, and the coefficient of variation (CV) of pseudobulb weight was higher than the CVs of other phenotypic traits (Fig. 2, Table 1). In detail, the weight of the oneyear-old pseudobulb of different populations varied from 5.02 to 86.08 (g), and the CV (%) was 69.55. Equal group classification method reveals that population categorized in class 5 represents the maximum level of pseudobulbs weight, and GXHZ, HBXG and HBLT populations were categorized in class 5, while the last five populations were HBXS, HBSY, HBWF (B. ochracea), SXXX and GXLY, which were categorized in class 1. For the two-year-old pseudobulb, the weight in different populations varied from 12.58 to 176.35 (g), and the CV (%) was 66.25. GXHZ, HBXG, HBLT, AHBZ and GZZY were ranked as the top five populations, which were categorized in class 5, while the last five populations were SCLS, HBSY, HBWF (B. ochracea), SXXX and GXLY. For the three-year-old pseudobulb, the weight in different populations varied from 23.13 to 329.42 (g), and the CV (%) was 72.34. GXHZ, AHBZ, HBXG and HBLT were categorized in class 5, while the last five populations were SYZS, HBSY, HBWF (B. ochracea), SXXX, and GXLY. For the four-year-old pseudobulb, the weight in the different populations varied from 31.63 to 442.77 (g), and the CV (%) was 71.06. HBLT, AHBZ, HBWH and HBSN were categorized in class 5, while the HBWF (B. ochracea), HBSY, SXXX, HNCL and GXLY populations showed lower weights. For pseudobulb weight, HBLT, AHBZ, HBXG, GXHZ and HBSN populations can be prioritized (Table S3). In general, except in GXHZ, HNCL and SCWY, the weight of the pseudobulbs increased with the growing period, while the average annual growth rates of the pseudobulb were 2.42, 1.69 and 1.31, respectively, showing a decreasing trend. There was a significant difference between the two species, and the average one-, two-, three-and four-year-old pseudobulb weights of B. ochracea were 19

Determination of active ingredients content
The polysaccharide, total phenol and militarine contents of pseudobulbs were determined among different species and populations. The results showed that significant differences were observed among different populations, with a minimum F value of 284.97 and P values much lower than 0.05 (Table 1). Overall, there were no differences in the content of the main active ingredients among different species and different planting years (Fig. 2c).
Most populations had an excellent phenotype in terms of polysaccharide content, which varied from 8.52 to 59.77 (%). There are 8, 12, 15 and 13 populations accounted for more than 50% of the one-, two-, three-and four-year-old pseudobulbs, and 17, 12, 13 and 9 populations ranged from 40 to 50%, respectively (Table 3). Equal group classification shows that  Table S1. This map was drawn by ArcGIS software based on longitude and latitude in population information (Table S1) (Table S5). Obvious variation was observed in the content of total phenol, which varied from 0.90 to 5.00 (%). There were 14, 11, 16 and 13 populations that were higher than 4% of the one-, two-, three-and four-year-old pseudobulbs, and 12, 15, 12 and 12 populations ranged from 3 to 4%, respectively ( Table 3). The average total phenol content of one-, two-, three-and four-year-old pseudobulbs of B. ochracea were 3.96, 3.57, 3.89 and 3.35 (%), which were lower than the 3.98, 4.14, 4.33 and 3.43 of B. striata, respectively. Therefore, the order of the contents of total phenol in the pseudobulb of B. striata was three-year-old > two-year-old > one-year-old > four-yearold, while that in B. ochracea was one-year-old > threeyear-old > two-year-old > four-year-old ( Fig. 2c). In general, the top five populations in class 5 with the highest average total phenol contents in pseudobulbs of four years were JSNJ, HBYL, HBXG, SCWY and HBLT, while YNBS and HBJM populations were categorized in class 1 (Table S6).
The average milliarine contents of the one-, two-, threeand four-year-old pseudobulbs of B. ochracea were 1.43, 1.62, 1.85 and 1.64 (%), which were generally higher than those of B. striata (1.44, 1.39, 1.45 and 1.60 (%), respectively). That is, the order of the militarine contents of the pseudobulbs of B. striata was four-year-old > threeyear-old > one-year-old > two-year-old, while that in B. ochracea was three-year-old > four-year-old > two-yearold > one-year-old (Fig. 2c). Militarine contents varied from 0.51 to 4.12 (%), 4, 4, 6 and 5 populations accounted for more than 2% of the one-, two-, three-and fouryear-old pseudobulbs, and 7, 5, 9 and 12 populations Table 2 The fresh weight of pseudobulbs of 28 populations of B. striata and 5 populations of B. ochracea No. 1-28 are B. striata, No. 29-33 are B. ochracea. Red represents the top five populations of pseudobulb weight in different years, while green represents the last five populations. Values represent the mean ± SD of ten individuals, and different lowercase letters with superscripts indicate significant differences between populations with the parameter of P < 0.05 accounted for 1.5% to 2%, respectively (Table 3) (Table S7).
GRA of the pseudobulb weight, pseudobulb growth ratio and main active ingredient content showed that there are 2 populations greater than 0.6 and 15 populations greater than 0.5. Meanwhile, HNSZ, AHBZ, HBLT, HBSN and JSNJ in class 5 were ranked as the top five, indicating that these populations have good characteristics in comprehensive yield and quality, while the last five populations of class 1 were SXXX, SYZS, GXLY, HBJM and YNBS (Table 4, Table S8). Meanwhile, HCA showed that HBWH, HBSN, AHBZ and HBLT belonged to group I, and GZZY, HBXG and GXHZ belonged to group II when the parameters were chosen to divide the population into 5 categories (Fig. S1).

Correlation analysis of phenotypic traits, pseudobulb weight and active ingredient content
To investigate which phenotypic traits may determine the yield of pseudobulb and active ingredient content, and whether there is a correlation between the yield of pseudobulb and active ingredient content, Pearson correlation analysis was performed for two-, three-and four-year-old plants among 13 indexes, including 9 morphological indexes of aboveground parts, pseudobulb weight indexes and 3 active ingredient indexes. The results showed that except for the negative correlation between the militarine content and other indexes, the other 12 indexes were positively correlated (Fig. 3).
The correlation coefficient of plants in 2018 ranged from -0.401 to 0.861. The correlation coefficients between plant height and leaf blade length, and between inflorescence height and inflorescence length were greater than 0.8 (0.806 and 0.861, respectively). There were 4 other correlation coefficients greater than 0.7, 3 correlation coefficients greater than 0.6, 5 correlation coefficients greater than 0.5, and 7 correlation coefficients greater than 0.4, and all of them showed extremely significant differences with P values of less than 0.01. Overall, there were significant correlations between phenotypic traits except leaf number, flower number and fruit number. There were significant correlations between pseudobulb weight and plant height, leaf blade width, stem diameter, and total phenol content. There were significant correlations between polysaccharide content and leaf blade width, stem diameter, and total phenol content. There were significant correlations between the total phenol content and leaf blade width, stem diameter, pseudobulb weight, and polysaccharide content. There was a significant negative correlation between militarine content and leaf blade length (Fig. 3a).
The correlation coefficient of plants in 2019 ranged from -0.282 to 0.920. The correlation coefficient of inflorescence height and inflorescence length was 0.920. The Table 3 Active ingredient contents of 28 populations of B. striata and 5 populations of B. ochracea No. 1-28 are B. striata, No. 29-33 are B. ochracea. Red represents the top five populations in terms of active ingredient contents in different years, while green represents the last five populations. Values represent the mean ± SD of three biological replicates, and different lowercases letters with superscript indicate significant differences between populations with the parameter of P < 0.05 correlation coefficients of plant height with inflorescence height and inflorescence length were greater than 0.8. There were 6 other correlation coefficients greater than 0.7, 7 correlation coefficients greater than 0.6, and 7 correlation coefficients greater than 0.5, and all of them had extremely significant differences with a P value less than 0.01. Overall, there were significant correlations between phenotypic traits except leaf number, and fruit number. There were significant correlations between pseudobulb weight and all phenotypic traits except for fruit number; however, the three indexes with the strongest correlations were the same as those in 2018, which were plant height, leaf width and stem diameter. The correlation between polysaccharide content and phenotypic traits was lower than that in 2018, and the same trends were observed for total phenol and militerine. However, there was still a significant correlation between the polysaccharide and total phenol content (Fig. 3b).
The correlation coefficient of plants in 2020 ranged from -0.331 to 0.931, and the highest value was obtained for inflorescence height and inflorescence length. There were 2 other correlation coefficients greater than 0.8, 3 correlation coefficients greater than 0.7, 6 correlation coefficients greater than 0.6, and 6 correlation coefficients greater than 0.5, and all of them had extremely significant differences with a P value less than 0.01. Overall, there were significant correlations between phenotypic traits except leaf number and fruit number. There were significant correlations between pseudobulb weight and plant height, leaf blade width, leaf number, and stem diameter. There were significant correlations between polysaccharide content and leaf blade width, stem diameter, total phenol content. There were significant correlations between the total phenol content and leaf blade width, and polysaccharide content. There was also a nearly negative correlation between the militarine content and other indexes (Fig. 3c).

Discussion
Stable yield and quality are the most important value indicators of medicinal materials, and good germplasm resources are the basis for ensuring the quality of medicinal materials. In our research, we continuously evaluated the main value species, including B. striata and B. ochracea, for 4 years. Based on our yearly results, we screened excellent populations with higher yields in terms of four-year-old pseudobulb weight ( Table 2, Table S3-4) and higher contents of polysaccharide, total phenol and militarine compounds, respectively ( Table 3, Table S5-7), indicating that these populations were good germplasm resources for the single-component utilization of Bletilla. Furthermore, HNSZ, AHBZ, HBLT, HBSN and JSNJ populations in class 5 showed good characteristics in terms of pseudobulb weight, the pseudobulb growth ratio and the main active ingredients contents according to GRA, indicating that these populations showed higher quality with good values of all the indicators ( Table 4, Table  S8). Due to the difference in calculation method, PCA showed that HBWH, HBSN, AHBZ and HBLT in group I, and GZZY, HBXG and GXHZ in group II were high placing. Combining the two analyses, AHBZ, HBLT and HBSN populations can be prioritized.
Bletilla is not strictly geoauthentic medicinal plant. The different geographical populations from different  weight of pseudobulbs and the contents of three active ingredients among different populations, whether near Shennongjia Forest District or Enshi Tujia and Miao Autonomous Prefecture, and there were abundant variations among populations. For example, the weights of the pseudobulbs of the HBXS and XSSC populations were significantly lower than those of HBSN, and the polysaccharide and total phenol contents in HBXS were higher than those in HBSN. The content of militerine was highest in HBLC, while that in HBXE and HBHF was very low, and the HBHF population had high contents of polysaccharide and total phenol. The determination of total polyphenol in B. striata from 16 areas of northwestern of Hubei Province showed that the content of total polyphenol varied in different areas. The analysis of the polysaccharide content of 18 samples of B. formosana from 8 populations in Baoshan City of Yunnan Province showed that there were significant differences among different samples. The polysaccharide contents of B. striata from different areas of Guizhou and Anhui Province were also significantly different [18]. All these studies showed that it is necessary to understand the resource distribution, habitat and situation of Bletilla and then to evaluate germplasm resources, which can provide a theoretical basis for the protection and utilization of Bletilla. Traditionally, although other species of Bletilla are also used as local medicinal plants, B. striata is the only species officially approved for use [3]. Considering the high medicinal value, high demand and gradually diminishing wild resources of B. striata, to protect traditional resources and develop new medicinal plant resources, it is necessary to conduct in-depth research on B. striata and its related species, to excavate alternative resources of B. striata and reduce the waste of effective ingredients and resources. However, the distribution of B. sinensis is narrow and limited to the south of Yunnan Province in China. The pseudobulb size of B. formosana is significantly smaller than the pseudobulb size of B. ochracea and B. striata, and the polysaccharide content of B. formosana was slightly lower than the polysaccharide content of B. ochracea and B. striata [29]. The pseudobulbs of B. ochracea are larger, and its flowering period is later than the flowering period of B. striata. Currently, Bletilla is widely planted as a bonsai plant, so in addition to ornamental plants, can underground pseudobulbs of B. ochracea also be harvested as medicinal plants?
We chose 5 populations of B. ochracea and 28 populations of B. striata to demonstrate the phenotypic quality of B. striata with higher average pseudobulbs weights, while there was no difference between the two species in terms of polysaccharide, total phenol and militarine contents. These results suggested that B. ochracea could be accepted as an alternative resource to B. striata.
Efficient screening of excellent genetic resources is important for new germplasm creation and breeding, and efficiently discovering the phenotypic indicators significantly related to yield and quality is especially critical for perennial herbs [30]. Usually, B. striata exhibits the best comprehensive traits and is harvested in the third year after cultivation [31,32], which requires at least 3 years to estimate the genetic material and preliminarily screen out a quality germplasm. In this paper, Pearson correlation analysis of two-, three-and four-year-old plants showed that there were significant correlations between pseudobulb weight and other phenotypic traits, especially plant height, leaf blade width and stem diameter, which implies that plants with excellent growth can produce better pseudobulb yields (Fig. 3). Interestingly, there was a strong significant correlation between the polysaccharide and total phenol contents, and polysaccharide and total phenol showed a similar trend to pseudobulb yield, which was significantly related to leaf blade width and stem diameter. These results may provide important and more direct indicators for quickly screening and breeding potential good germplasms. The uncommon positive correlation between yield and the main quality factors may provide hopeful prospects for excellent germplasm breeding of Bletilla herbs; of course, this positive relation and the genetic foundation of this relation need further demonstration and in-depth investigation [26]. In addition, the militarine content was almost negatively correlated with other indexes, which indicated that the accumulation of militarine may occur more often under the condition of limited growth. We isolated one-, two-, three-and four-year-old pseudobulbs of four-yearold plants from four populations. The militarine content in pseudobulbs shows an increasing trend year by year, implying that sufficient growth years are required to accumulate more militarine content, and further research is needed to reveal the synthesis and metabolic mechanisms of militarine in Bletilla plants (unpublished). Next, integrating phenomics, genomics and metabonomics approaches for high-throughput mining of the key phenotypes and revealing the genetic mechanism is our aim in recent years.
The growth capacity of yield characteristics can also be selected to evaluate different genetic resources. Populations of HBSN, HBWH, and AHBZ showed excellent yield phenotypes with high yield and growth capacity, and the pseudobulb weight of the HBSN population could even reach 7.29 and 10.13 times in the third and fourth years after cultivation, respectively. Other populations such as HNSZ, SXZP (B. ochracea), HBXS, and HBJM, also showed a good growth capacity. Our results indicated that the mean annual growth ratios of pseudobulbs in 3 years were 2.42, 1.69 and 1.31, which gradually decreased with the cultivation years. Even in some populations, such as GXHZ, HNCL and SCWY, the pseudobulb weight in the fourth or even third year began to decrease. Yearly content variation showed that the average contents of both polysaccharides and total phenol in the one-, two-, three-and four-year-old pseudobulbs of Bletilla varied in the following order: three years old > two years old > one year old > four years old, which confirmed the conclusion regarding the best harvest period of Bletilla [31,32]. However, the average militarine content of one-, two-, three-and four-yearold pseudobulbs of Bletilla varied in the following order: four-year-old > three-year-old > one-year-old > two-yearold, which again indicated differences in polysaccharide and total phenol, and the synthesis and metabolic mechanisms of militarine in Bletilla plants need to be revealed.

Conclusion
There were abundant phenotypic variations and significant differences among different populations in the morphological phenotypes, pseudobulb weight and main active ingredient contents. Bletilla is not strictly geoauthentic medicinal plants. Excellent germplasm resources may be distributed in distant geographical locations. There may be a large difference between different geographical populations even in the same region. B. ochracea could be accepted as an alternative resource to B. striata. AHBZ, HBLT and HBSN populations showed good prospects for industrial development, presenting higher quality in terms of yield and main active ingredient content. Pseudobulb yield, polysaccharide and total phenol content are positively correlated with phenotypic traits. Militarine content is negatively correlated with almost all indexes. Plant height, leaf width and stem diameter may be important indicators of potential excellent germplasms. The best harvest period of Bletilla is the third year after cultivation. These results provide important information required for the efficient screening and utilization of Bletilla germplasm resources.

Experimental material
Thirty-three natural populations of Bletilla (including 28 B. striata populations and 5 B. ochracea populations) were randomly collected from its main natural range across ten provinces in China (Table S1; Fig. 1). In 2017, 30 individuals of new pseudobulb of the year with similar size were collected from each population without affecting the continuous growth of wild resources and transported to a germplasm resources nursery in the Wuhan Botanical Garden (WBG) of the Chinese Academy of Sciences (CAS) for common-garden cultivation with the same soil formula of 20% sand, 30% humus and 50% garden soil.
WBG is an institution for species conservation and scientific research. The ex-situ conservation of plant resources is the main job of WBG. Introducing species from other area only requires a work certification issued by WBG and does not require government or local administrative permission. The representative specimens of plants from each population were identified by XiaoDong Li

Phenotypic measurements
The weight of the pseudobulb was measured and marked as one-year-old in winter 2017. The morphology of the aboveground parts and the weight of the pseudobulbs of 10 individuals were measured every year for the next three years (2018-2020).
The recorded morphological parameters of the aboveground parts included plant height (cm), leaf blade length (cm), leaf blade width (cm), total leaf number, stem diameter (mm), inflorescence height (cm), inflorescence length (cm), total flower number and total fruit number. Leaf phenotypes and stem diameter were measured in August or September. Plant height was the vertical distance between the top of the highest leaf and the ground. Leaf length and leaf width were measured using the largest leaf. The phenotypes of the inflorescences were measured in April or May. Inflorescence height was the vertical distance from the base to the tip of the inflorescence stem.
In early November 2018, the withered aboveground parts of 10 individuals were removed, and the pseudobulbs (two years old) were excavated and cleaned to measure the weight. The corresponding fraction of pseudobulbs was removed for the subsequent determination of polysaccharide, total phenol and militarine content. The same method was used to measure the weight of pseudobulbs in 2019 (three years old) and in 2020 (four years old) using the other 10 individuals.

Polysaccharide content analysis
After removing the fibrous roots, the sampled pseudobulbs were sliced, steamed in a water bath, and dried in the shade according to the Chinese Pharmacopeia. Dried pseudobulb slices were crushed for polysaccharide content, total phenol content and militarine content analyses. All the samples were collected in early November 2017, 2018, 2019 and 2020.
The phenol-sulphuric acid method was used to analyze polysaccharide content as described previously with minor modifications [23]. Approximately 2 g of dried pseudobulb powder was treated with 50 mL H 2 O for 120 min at 100 °C twice using the conventional heating reflux extraction method. The 100 mL extract was concentrated to 10 mL by rotary evaporation, 25 mL of 95% ethanol solution was added to half of the extract solution (5 mL) for precipitation, and they had to rest overnight at room temperature. After centrifugation (3000 rpm/ min) for 15 min, the precipitate was dissolved in 25 mL H 2 O to obtain the crude polysaccharide solution and quantified by the phenol-sulfuric acid method. 10 mg of standard glucose was dissolved in 100 mL H 2 O at a concentration of 100 µg/mL and diluted to six different concentrations of 0, 20, 40, 60, 80 and 100 µg/mL. 1 mL glucose solution, 1 mL 5% phenol and 6 mL concentrated sulfuric acid were added into a plug test tube in turn. After shaking and mixing, the samples were left at 25 °C for 40 min, and the absorbance was determined at 484 nm (A 484 ). The regression equation was A 484 = 0.0088 C (µg/mL) + 0.0122, r = 0.9998 (n = 6). Then, 1 mL polysaccharide solution was measured by the same method as glucose solution. The absorbance was compared with the regression equation, multiplied by the dilution multiple (50) and divided by the weight of dried pseudobulb powder. The results were expressed as polysaccharide content (%).

Total phenol content analysis
The sampling procedure was the same as above. Total phenol content was determined by Folin-Ciocalteu's phenol reagent [33]. Gallic acid was used as the standard to draw the regression equation, and 100 mg gallic acid was dissolved in 100 mL 60% ethanol at a concentration of 1000 µg/mL and diluted to five different concentrations of 10, 20, 30, 40 and 50 µg/mL. Then, 1 mL of gallic acid solution and 5 mL 10% (v/v) Folin-Ciocalteu's phenol reagent were added into a test tube to react for 10 min, and 5 mL 2% (m/v) Na 2 CO 3 was added. After shaking, the sample was allowed to rest at 25 °C for 1 h, and the absorbance of each sample was determined at 760 nm (A 760 ). The regression equation was A 760 = 0.0132 C (µg/ mL) − 0.0004, r = 0.9999 (n = 6). Approximately 2 g dried pseudobulb powder was treated with 80 mL 60% ethanol for 120 min at 90 °C by using the conventional heating reflux extraction method. Cooling down to RT, the solution was made up to 80 ml with 60% ethanol and then filtered using filter paper. The residue was extracted with 60% ethanol again, and the filtrates were mixed as the test solution, which was measured by the same method as the standard sample, to obtain the A 760 , compared with the regression equation, multiplied by the dilution multiple (160) and divided by the weight of dried pseudobulb powder. The results were expressed as total phenol content (%).

Militarine content analysis
The sampling procedure was the same as above. HPLC was used to measure the militarine content according to the Chinese Pharmacopeia 2020 edition. 11.4 g millitarine was dissolved in 6 ml 52.9% ethanol to obtain a 1.9 mg/mL standard solution, which was diluted to six different concentrations of 1900, 950, 475, 237.5, 118.75 and 59.375 µg/mL. The standard solutions were analyzed by reversed-phase HPLC using a mobile phase of acetonitrile and water with 0.1% phosphoric acid (elution ratio: 22:78, v/v; wavelength: 223 nm; flow rate: 1 mL/min; injection volume: 10 µL) and the column oven temperature was set at 30 °C. The regression equation of the standard curve was peak area = 22,287,199 C (µg/mL) + 137,858, r = 1 (n = 6). Approximately 0.2 g dried pseudobulb powder was extracted with 25 ml 52.9% ethanol at 30 °C for 30 min by an ultrasound system (power = 300 W, frequency = 37 Hz). The solution was brought up to 25 ml with 52.9% ethanol and then filtered through filter paper. The filtrate was analyzed by the same method as the standard sample, compared with the regression equation, multiplied by the dilution multiple (25,000) and divided by the weight of the dried pseudobulb powder. Ultimately, the results were expressed as militarine content (%).

Prioritization and elite identification
The value of each parameter, including pseudobulbs weight, pseudobulb growth ratio, polysaccharide content, total phenols content, militerane content and GRA (Tables 2, 3, and 4), in each population was categorized into five equal classes, and population in class 5 with higher values were considered as a prioritized group of the population [34]. At the same time, the sum of the one-, two-, three-and four-year-old pseudobulb weight, and the average of polysaccharide, total phenols and militerane content in one-, two-, three-and four-yearold pseudobulb were used for prioritization and elite identification.

Gray relation analysis
The means of the values in Tables 2 and 3 were used as the raw data for GRA. For pseudobulb weight, in addition to considering the weight of the pseudobulb measured directly in every year, pseudobulb growth ratios of one, two and three years were also used as indexes. A total of 19 indexes constituted the evaluation matrix, {X ik } (i = 1, 2, 3, … m; k = 1, 2, 3, … n, m = 33, n = 19) and calculated the average value (‾X K ) of each index. The original data were normalized by the formula Y ik = X ik /‾X K . The maximum and minimum values of each index were selected as the optimal {Y sk } and the least {Y tk } reference sequences respectively. The correlation coefficients relative to the optimal and the least reference sequences were calculated with the formulas (∆min + ρ∆max)/(Y sk -Y ik + ρ∆max) and (∆min + ρ∆max)/(Y ik -Y tk + ρ∆max), respectively, where ρ is the discrimination coefficient, and its value is generally approximately 0.5, ∆min = min | Y sk -Y ik |, ∆max = max | Y sk -Y ik |. The uniform average of the correlation coefficient in each index was used to calculate the correlation degree relative to the optimal and the least reference sequences, which is expressed as r i(s) and r i(t) , respectively. The relative correlation degree (r i ) of the evaluation matrix {X ik } relative to the optimal and least reference sequences was defined as r i = r i(s) / (r i(s) + r i(t) ), which wes used to rank (Table 4).

Hierarchical clustering analysis (HCA)
The data used for GRA (pseudobulb weight, pseudobulb growth ratio, polysaccharide content, total phenol content and militarine content) was also used for HCA by IBM SPSS Statistics 25 software.

Correlation analysis
The means of the values in Table S2, Tables 2 and 3 were used to calculate the Pearson correlation coefficient with IBM SPSS Statistics 25 software, which was used to construct a heat map in Microsoft Office Excel 2019. * or ** represent significant differences with the parameters of P < 0.05 and P < 0.01, respectively.

Data analysis
Two-tailed ANOVA was performed to test significant differences in measured variables among different populations and different growth years. ANOVA and Pearson correlation analysis were carried out with IBM SPSS Statistics 25.
Abbreviation GRA : Gray relation analysis.