The G. lucidum from Fujian was provided by Fujian Xian-Zhi-Lou Technology Co., Ltd. The strain of G. lucidum was JB-A119-20191010. It was incubated for 7 days at 28°C on potato dextrose agar (PDA) (medium: potato, 200 g; dextrose, 20 g; agar, 20 g; and deionized water, 1 L). The inoculated substrate was incubated at 28°C for mycelial colonization. The air temperature in the cultivation room was maintained at 28°C, and the relative humidity was set at 70%–80%. The grown mycelium was homogenized after 15–20 days of cultivation. The culture medium was autoclaved at 121°C for 2 h. The fruiting body samples were cultivated by the solid method using wood (the wood of Castanopsis sclerophylla used for cultivating G. lucidum was taken from local forests). Before planting, the cultivation site had been aired, cleaned, and slightly lime-disinfected and had proper ventilation (CO₂ 0.1%–0.3%). The experiment was conducted at the three growth stages of bud development, growth, and maturing. Bud development stage: After inoculation, the mycelia of G. lucidum began to germinate. Primordium culture conditions were maintained at 28°C–30°C, 90%–95% relative humidity, 0.03%–0.1% CO₂, and a light intensity of 100 lux. The initial sampling was conducted approximately 35–40 days after inoculation, which is when the bud-developing stage and differentiation had begun. Growth stage: The second sampling was conducted at the growth stage, approximately 55–60 days after inoculation, which is when the primordial grows long and stripes form. The growth stage began after elongation. The color of the caps deepens as they grow, demonstrating a deep center and white edges. The culture conditions of fruiting bodies were maintained at 25°C–28°C, relative humidity at 80%–95%, CO₂ of 0.03%–0.1%, and a light intensity of 150–200 lux. Maturing stage: The last sampling was performed at the maturing stage approximately 70–80 days after the inoculation. The culture conditions of fruiting bodies were maintained at 25°C–28°C, relative humidity at 80%–95%, CO₂ at 0.03%–0.05%, and a light intensity of 150–200 lux. G. lucidum at different stages, when the white edge of the cap disappears and hardens, is shown in Figure 1 .

G. Lucidum collected from Zhejiang, Anhui, and Jilin was grown on wood logs and purchased from the local “Daodi” market (maturity stage). The origin of G. lucidum includes latitude, longitude, and specific parameters of the planting process as shown in Table 1 . The pictures of G. lucidum from different origins and the specific sampling locations are presented in Figures 2 and 3 , respectively. The Pseudostellaria heterophylla samples used for heterogeneous validation were collected from Anhui (n = 3) and Fujian (n = 3). All the samples were dried at 45°C in a blast drying oven (9240A, Shanghai Jing Hong Laboratory Instrument Co., Ltd.). The dried samples were then ground into a fine powder with a ball mill, sieved through a 100-mesh sieve, and then subjected to stable isotope ratio analysis.

The stable isotope ratios (C and N) and the carbon and nitrogen contents were measured with stable isotope mass spectrometry (Element analyzer-isotope ratio mass spectrometers (EA-IRMS), Thermo). The EA-IRMS conditions were set as follows: the reactor temperature was 960°C, and gases produced by the oxidation reactor were changed into He in the reducing furnace. He was then used as a carrier gas at a flow rate of 90 mL·min⁻¹. The gases were subsequently separated and measured by EA-IRMS. Reference materials, including IAEA-600 and USGS40, were used for calibration. The sample powders were accurately weighed at 1 mg and then packed in a tin cup for the analysis of EA-IRMS.

The stable isotope ratios (O and H) were measured with elemental analysis isotope ratio mass spectrometry (high-temperature pyrolysis element analyzer and gas isotope ratio mass spectrometer (HT-IRMS), Thermo). The HT-IRMS conditions were set as follows: the reactor temperature was 1,380°C, He was used as carrier gas at a flow rate of 100 mL·min⁻¹, and the column temperature was 600°C. Reference materials, including EMA-P1, EMA-P2, and IAEA-601, were used for calibration. The sample powders were accurately weighed at 3 mg and then packed in a tin cup for the analysis of HT-IRMS.

Stability checks of the used reference gases were continuously performed by determining International Atomic Energy Agency standards with defined ¹³C/¹²C, ¹⁵N/¹⁴N, ¹⁸O/¹⁶O, and D/¹H ratios. Isotope ratios were expressed in per mil (‰) deviation relative to the V-PDB, V-SMOW, and N₂ international standards. Values were calculated using the following formula:

where R_sample is the ratio of the heaviest to the lightest stable isotope in the sample, and R_standard is the ratio of the heaviest to the lightest stable isotope in the international reference standards, namely, V-PDB, air, and V-SMOW for δ ¹³C, δ ¹⁵N, and δ ¹⁸O, and δD values, respectively.

The results were presented as mean values ± SD (standard deviation). One-way analysis of variance was utilized for the statistical analysis (ANOVA). The FLDA was examined with SPSS 26.0 (SPSS Inc., Chicago, IL, USA). Chemometric analyses, including principal component analysis (PCA) and orthogonal partial least squares-discriminant analysis (OPLS-DA), were applied using the SIMCA-P Version 14.1 (Umetrics AB, Umea, Sweden).

Six variables, including C, H, O, and N stable isotopes and the C and N contents in G. lucidum samples from different provinces (ZJ, AH, JL, and FJ), are compiled in Table 2 . G. lucidum isotope ratios showed δ ¹³C values ranging from -29.61 to -12.01, δD values ranging from -31.9 to 4.1, δ ¹⁸O values ranging from 21.3 to 24.7, and δ ¹⁵N values ranging from -4.36 to 3.59. The distribution of four stable isotopes is shown in Figure 4 to demonstrate the differences among origins.

Using the data of isotope ratios (δ ¹³C, δD, δ ¹⁸O, and δ ¹⁵N), an unsupervised multivariate data analysis method, that is, PCA, was performed to minimize dimension and visualize the separation effects of geographical origins. The 2D score plot of the PCA divided the samples into four independent clusters ( Figure 5A ). The results of the classification were consistent with those of the original data. The results are based on the 95% confidence interval ( Figure 5B ).

OPLS-DA is a modification of the PLS-DA method, where the variation in the data (X matrix) is separated into two parts: one that is linearly related to the response variable (Y matrix) and one that is orthogonal or unrelated. In SIMCA version 14.1 by Umetrics, the OPLS-DA method is implemented to allow multigroup data analysis by performing several binary (one-against-all) comparisons. In the context of a four-group study, four separate OPLS-DA models would be constructed, each focusing on the discrimination of one group against the remaining three. Therefore, these models can then be collectively interpreted to understand the differences among the four groups. OPLS-DA models in SIMCA 14.1 also allow users to visualize and interpret the predictive and orthogonal components separately, which aids in effectively understanding the underlying structure of the dataset. Therefore, OPLS-DA is a supervised discriminant analysis model that can assist in understanding inter-class separation and identifying potential categorical variables, increasing its suitability for multivariate discriminant analysis. The OPLS-DA model was developed utilizing stable isotope ratios as variables to visualize and investigate their intrinsic correlation and explore their contribution to the origin discrimination of G. lucidum. The OPLS-DA 2D score plot showed that samples were divided into four categories ( Figure 6A ). The results are based on the 95% confidence interval ( Figure 6B ). The model can explain 91.6% of the variation (R ²X = 0.916, R ²Y = 0.726), demonstrating a high predictive capability (Q ² = 0.67) and showing a good fit and predictive capability of the training set. Theoretically, in OPLS-DA, VIP values larger than 1 indicate that variables are essential. The VIP ( Figure 6C ) suggested that the δ ¹³C value was the most important contributor to the discrimination (VIP > 1). Cross-validation with 200 permutation tests showed that the OPLS-DA model has good robustness (R² = 0.0295, Q² = −0.276, Figure 6D ).

G. lucidum samples (ZJ = 2, AH = 2, JL = 2, and FJ = 3) were used for external source validation. Figure 7 shows the predicted and actual values of the PLS validation according to the isotope ratios (RMSEE = 0.62, RMSECV = 0.66). The results revealed that the model could identify nine external samples in their classified group with 100% accuracy.

Six Pseudostellaria heterophylla heterogeneous samples were incorporated at random, including three samples from Fujian and others from Anhui, to further confirm the reliability of the model. δ ¹³C, δD, δ ¹⁸O, and δ ¹⁵N were acquired by using the same isotope ratio mass spectrometry for verification of other heterogeneous samples ( Figure 8 , RMSEE = 0.61, RMSECV = 0.60). Figure 8 shows the classification of G. lucidum and Pseudostellaria heterophylla samples with a cutoff value of 4.5. The result indicated that G. lucidum has its own unique stable isotopes. Moreover, the results revealed that the model could identify other heterogeneous samples (Pseudostellariae Radix) in other groups with 100% accuracy. By contrast, these samples (G. lucidum and Pseudostellaria heterophylla) were again analyzed by PCA ( Figure S1 ). The analysis results showed that G. lucidum was clearly isolated from samples of Pseudostellaria heterophylla. Therefore, G. lucidum is efficiently and precisely identified using the established model.

FLDA was employed to establish a predictive model to facilitate the clustering of group members based on observed characteristics. The discriminant function based on linear combinations of predictor variables was utilized to discriminate further and classify the unknown members. The FLDA model of the geographical origin of G. lucidum was established, and the classification results are shown in Table 3 . The FLDA function is as follows: Y(ZJ) = −18.131 δ ¹³C + 3.131 δ ¹⁵N + 70.205 δ ¹⁸O − 0.316 δD − 932.103, Y(AH) = −47.559 δ ¹³C − 3.887 δ ¹⁵N + 62.028 δ ¹⁸O − 6.324 δD − 1371.246, Y(JL) = −46.931 δ ¹³C − 2.976 δ ¹⁵N + 57.780 δ ¹⁸O − 7.040 δD − 1275.817, and Y(FJ) = −43.875 δ ¹³C − 1.462 δ ¹⁵N + 67.994 δ ¹⁸O − 4.505 δD − 1394.239. The FLDA linear discriminant with stable isotopes as discriminant variables has been evaluated with back-replacement accuracy of 100% and cross-validation reaching 95.7%.

Scatter plots were used to evaluate the distribution of the geographical origins to distinguish the discrepancies explicitly. As shown in Figure 9 , geographical origins could be successfully divided into four groups. The discriminant function is based on linear combinations of predictor variables to discriminate further and classify the unknown members. The model for geographical origin discrimination was used to verify the origins of nine samples outside the model. The nine G. lucidum samples were accurately recognized, demonstrating the precision and consistency of the model.

Thus far, studies that have exhaustively examined the relationship between stable isotopes and geographical origins at various phases are currently unavailable. Four stable isotopes and C and N contents are shown in Table 4 and Figure 10 , respectively. G. lucidum stable isotope showed δ ¹³C values ranging from −29.61 to −24.30, δD values ranging from −31.8 to 4.1, δ ¹⁸O values ranging from 22.8 to 25.0, and δ ¹⁵N values ranging from −4.22 to 1.53. The average δD values of samples at the maturing stage were highest (−3.0‰), followed by the growth and bud development stages.

The correlation between carbon and nitrogen can reflect the nutrient acquisition of G. lucidum during the growth process, which plays an important role in its growth (Hsieh and Yang, 2004). Correlations of δ ¹³C and δ ¹⁵N values with C and N contents in G. lucidum are shown in Table 5 . A significantly negative relationship is observed between the C/N ratio and the N content (r = −0.894, p < 0.001). The relationship between the carbon-to-nitrogen ratio and nitrogen content is one of mutual inhibition.

The model based on stable isotope ratios fails to distinguish samples at different stages and overlaps in some areas ( Figure 11A ). Compared with the model of stable isotopes (four variables, including δ ¹³C, δ ¹⁵N, δD, and δ ¹⁸O), the PCA model, after increasing C% and N%, showed that G. lucidum samples at each stage were highly centralized ( Figure 11C ). Different stages were clustered on the left of PC1. The results are based on the 95% confidence interval ( Figures 11B, D ). The OPLS-DA score plot showed that the samples were divided into three categories, with significant differences ( Figure 12A ). The model could explain 99.0% of the variation (R²X = 0.991, R²Y = 0.74) and had high prediction capability (Q² = 0.638). The results are based on the 95% confidence interval ( Figure 12B ). The VIP of six variables showed that δ ¹⁸O, N%, and δ ¹⁵N are the most important contributors to the discrimination (VIP > 1, Figure 12C ). Cross-validation with 200 permutation tests showed that the OPLS-DA model has good robustness (R² = 0.0915 and Q² = −0.508, Figure 12D ).

Growth stage samples were used for the model verification. Figure 13 shows the predicted and actual values of the PLS model according to the four isotope ratios (δ ¹³C, δD, δ ¹⁸O, and δ ¹⁵N) and C% and N% (RMSEE = 0.27, RMSECV = 0.34). All external sample validations of G. lucidum were classified correctly.

The eigenvalues of the sample correlation matrix were calculated, and the eigenvalue demonstrates the importance of each new factor in explaining data variability. The contribution of six variables, namely, δ ¹³C, δD, δ ¹⁸O, δ ¹⁵N, C%, and N%, was analyzed by PCA at the PCA model, and the results are presented in Table 6 . Three major factors generally accounted for 87.578% (factors 1, 2, and 3 explained 47.078%, 22.266%, and 18.234% of the data variability, respectively) of the six variations inherent in the input model data. Among them, N% was strongly associated with factor 1 (0.878), and C% was highly correlated with factor 2 (−0.854), indicating that the contribution rates of C% and N% were highly correlated with the first two factors, thereby boosting the model’s accuracy, especially N% (see Table 7 ). The contribution of six variables was also analyzed using the VIP factor in the OPLS-DA model, and the results are presented in Figure 12C . δD, N%, and δ ¹⁵N VIP values >1 indicate that these markers contribute to differentiation at different stages. This finding also showed the importance of N% in the model.