Accurate Traceability of Stable C, H, O, N Isotope Ratios and Multi-Element Analysis Combined with Chemometrics for Chrysanthemi Flos ‘Hangbaiju’ from Different Origins

Abstract

Chrysanthemi Flos ‘Hangbaiju’ (HBJ) is a common Chinese medicinal material with the same origin as the medicinal and edible cognate plant in China, whose quality is seriously affected by the place of origin. In this study, four stable isotope ratios (δ¹⁵N, δ²H, δ¹³C, and δ¹⁸O) and 44 elements were detected and analyzed in 191 HBJ flower samples from six locations in China to trace the origin of HBJ. An ANOVA analysis of δ¹⁵N, δ²H, δ¹³C, and δ¹⁸O values, as well as milti-elements, showed that there were significant differences among the six places of origin. Partial least squares discriminant analysis (PLSDA) and one-class partial least squares discriminant analysis (OPLS-DA) models were established to trace the origin of HBJ from these six locations. The results showed that the classification effect of the PLSDA model is poor; however, the established OPLS-DA model can distinguish between products of national geographic origin (Tongxiang City, Zhejiang Province, China) and samples from other origins, among which Ni, Mo, δ¹³C, Cu, and Ce elements (VIP > 1) contribute the most to this classification. Therefore, this study provides a new method for tracing the origins of HBJ, which is of great significance for the protection of origin labeling of products.

1. Introduction

Chrysanthemi Flos (often called “Juhua” in China) is the dried capitulum of Chrysanthemum morifolium Ramat and is one of the most famous traditional Chinese medicines (TCM), and widely distributed in many countries in Asia [1,2]. As a common medicine and food homologous TCM, “Juhua” is widely used [3]. “Juhua” has the effect of dispersing wind and clearing heat, calming the liver and clearing the eyes, and clearing heat and detoxifying, and is mainly used for wind and heat colds, headache and dizziness, eye swelling and pain, eye faintness, and other symptoms [4]. In addition, it can be eaten as “Juhua” tea, “Juhua” wine, “Juhua” cake, “Juhua” soup, etc., all of which are very popular [5]. Chrysanthemi Flos ‘Hangbaiju’ (HBJ), the dried capitulum of Chrysanthemum morifolium Ramat ’Hangabaiju,’ is one of the most famous medicinal and edible varieties of “Juhua.” It is mainly produced in Tongxiang City, Zhejiang Province, China, and is a national geographical indication (GI) product in China [6]. HBJ has a long history of cultivation in China for more than 1500 years [7]. Its flowers are widely used to make healthy beverages and herbal teas, and are very popular among consumers for their various beneficial chemicals and unique aroma [8]. Therefore, in order to meet the market demand, many farmers have introduced HBJ from Tongxiang City for planting and cultivation. However, due to the influence of natural factors (soil factors, precipitation, altitude, etc.) or human factors (fertilization, planting methods, processing methods, etc.) in the production area, the quality of the introduced HBJ entering the market is not as good as that of the GI products produced in Tongxiang City. In order to protect the interests of consumers and ensure the quality of GI products, it is very necessary to trace the origins of HBJ.

In recent years, stable isotope ratio and multi-element analyses have been widely used in food and drug origin traceability, authenticity discrimination, variety identification, etc. [9,10,11]. Zhao et al. used four stable isotope ratios of δ¹³C, δ¹⁵N, δ²H, and δ¹⁸O, and seven element concentrations, to identify organic pork and conventionally-grown pork from four origins, and the results showed that the method could distinguish between organic pork and conventionally-grown pork from different origins [12]; McLeod et al. successfully used stable isotope ratios and trace elements to discriminate the authenticity of goat milk powder [13]. Wang et al. studied and analyzed the stable isotope ratios (δ¹⁵N, δ²H, δ¹³C, and δ¹⁸O) and 32 elements of Cassiae Semen tea, and found that stable isotope ratios and multi-element combined chemometric methods can identify the geographical origin of Cassiae Semen tea [14]. In terms of drug analysis, Du et al. used the stable isotope ratio and the content of various elements to discriminate the origin of three different varieties of Rhizoma Coptidis, conducted authenticity identification of true and false Rhizoma Coptidis, and found that the combination of stable isotope and multi-element analyses can accurately perform origin traceability and authenticity identification of Rhizoma Coptidis [15]. Fu et al. used δ¹⁵N, δ²H, δ¹³C, and δ¹⁸O, ⁸⁷Sr/⁸⁶Sr values and multiple elements to successfully trace the origins and discriminate the authenticity of Huangjing herbs [16]. These studies show that stable isotope ratios combined with multi-element analysis technology is a powerful method for origin traceability, authenticity discrimination, and variety identification in food and medicinal products [9,17].

Chemometrics, as a strong and effective analytical method, is often used in combination with other analytical techniques to perform differentiation and identification [18,19,20,21,22]. Therefore, it is also widely used in the tracing of the origins and identifying the authenticity of food and drugs.

In recent years, with the increase in people’s demand for HBJ, the research on it has gradually increased. Among studies on HBJ, there is a relatively large number of research reports on the chemical constituents of HBJ [23,24,25]. Stable isotope ratios often reflect useful information about the climate and anthropogenic environment a plant was cultivated in, while fingerprint information about plant growth environment such as the soil (including anthropogenic and geological factors) is often reflected by trace elements [16]. However, at present, there are only a few element-based studies on HBJ; these studies involve fewer elements and a smaller number of origins [26], and there is no single study related to the stable isotope ratios of HBJ. Therefore, in order to more systematically study the climatic factors, soil environments, and human influences on the growth of HBJ, more stable isotopic ratios, elements, and production areas should be considered.

The purpose of this paper is to analyze the stable isotope ratio and multi-element content of HBJ flowers from different origins, and to establish a simple, fast, and reliable method for tracing the origin of HBJ samples in combination with the chemometrics method. In order to achieve this research purpose, we studied the stable isotope ratios of C, H, O, and N, and 44 elements in different HBJ flowers, and combined these chemometrics analysis methods to identify six different origins of HBJ. The ultimate purpose of this study is to detect the authenticity and reliability of the information on HBJ products that flow into the market, to safeguard the legitimate rights and interests of consumers, and to protect the geographical indication of products.

2. Materials and Methods

2.1. Sampling and Preparation

“Juhua” samples were collected from local farmers or planting bases from October to November 2020, and samples from each production area were sampled according to the five-point sampling method. The details are shown in Figure 1. The collected fresh samples were air-dried in the shade to remove surface moisture, and then placed in a forced-air oven at 60 °C for 24–48 h until the samples were dry. They were then crushed, passed through a 50-mesh sieve, labeled according to the number of samples, and stored in a desiccator for later use.

2.2. Chemicals and Instruments

Nitric acid (analytical grade, Merck, Rahway, NJ, USA) and hydrogen peroxide solution (analytical grade, Tianjin Damao Chemical Reagent Factory, Tianjin, China) were the chemicals used in this study. The mineral element standard materials used in the analysis of 44 elements include GNM-M321850-2013 standard material, GSB 04-1724-2004 (Cs, 1000 μg/mL), GSB 04-1789-2004 (100 μg/mL), GSB 04-1750-2004 (Sc, 1000 μg/mL), and GSB 04-1751-2004 (1000 μg/mL) (National Nonferrous Metals and Electronic Materials Analysis and Testing Center, Beijing, China), all used to prepare standard curves, as well as GBW07697 (Rb, 1000 μg/mL, China Institute of Metrology, Beijing, China) to prepare working curves.

A total of 44 kinds of elements were analyzed by inductively-coupled plasma atomic emission spectrometry (ICP-OES, ICAP 7400, Thermo Fisher Scientific, Waltham, MA, USA), inductively-coupled plasma mass spectrometry (ICP-MS, Thermo Fisher Scientific, Waltham, MA, USA), and atomic fluorescence spectrometry (AFS, AFS-8500, Beijing Haiguang Instrument Co., Ltd., Beijing, China). The carbon and nitrogen stable isotope ratios of the samples were analyzed by an elemental analyzer (EA, HT 2000, Thermo Fisher Scientific, Waltham, MA, USA) coupled to a stable isotope ratio mass spectrometer (Delta V Advantage, Thermo Fisher Scientific, Waltham, MA, USA); hydrogen and oxygen isotopes were analyzed using stable isotope ratio mass spectrometry and an AUW120D electronic balance (Shimadzu Corporation, National, Kyoto, Japan).

2.3. Stable Isotope Ratio Analysis

The carbon and nitrogen stable isotope ratios of the samples were analyzed using an ICP-MS and a stable isotope ratio mass spectrometer (Delta V Advantage, Thermo, Waltham, MA, USA). Samples of about 2 mg were weighed out with a one-millionth electronic balance and wrapped in a tin cup for testing. The elemental analyzer conditions were as follows: the furnace temperature was 960 °C, the oxygen supply time was 3 s, the flow rate was 250 mL/min, and the helium flow rate was 50 mL/min. In the elemental analyzer, the carbon and nitrogen elements of the sample were oxidized to CO₂ and N₂, respectively. The column temperature was 70 °C. The isotope standard used IAEA-600 (δ¹³CV-PDB = −27.771‰, δ¹⁵NAir = +1‰) as the quality control for carbon stable isotope analysis. The test accuracies of carbon and nitrogen stable isotope ratios were 0.1‰ and 0.2‰, respectively.

Hydrogen and oxygen isotopes were analyzed by stable isotope ratio mass spectrometry. Samples that were about 0.3 mg were weighed and wrapped in a silver cup and placed in a sample pan, and then the whole sample pan was placed in a desiccator to dry for 24 h to reduce the influence of moisture in the air. The conditions of the pyrolysis analysis method were as follows: the temperature of the combustion furnace was 135 °C, the temperature of the chromatographic column was 90 °C, the flow rate of helium gas was 100 mL/min, and the flow rate of the injector purge was 100 mL/min. The hydrogen and oxygen elements of the sample were pyrolyzed by the pyrolysis tube to generate H₂ and CO in the carbon environment, and then sent to the mass spectrometer through a He flow to analyze the stable isotope ratio of hydrogen and oxygen. The stable isotope standards USGS 43 (δ²H_V-SMOW = −50.3‰, δ¹⁸O_V-SMOW = +14.11‰), USGS 54 (δ²H_V-SMOW = −150.4‰, δ¹⁸O_V-SMOW = +17.79‰), USGS 56 (δ₂H_V-SMOW = +17.79‰), and SMOW = −44‰ (δ¹⁸O_V-SMOW = +27.23‰) were used as quality control for data correction of hydrogen’s and oxygen’s stable isotopes (Table S1). The long-term laboratory analysis accuracies of hydrogen and oxygen stable isotopes were 1.2‰ and 0.3‰, respectively.

The stable isotope calculation formula is as follows:

δ(‰) = [(R_sample/_Rstandard) − 1]

Here, R_sample represents the heavy/light isotope ratio in the measured sample, and R_standard represents the standard heavy/light isotope ratio. δ(C) is relative to the internationally recognized Vienna Pee Dee Belemnite standard (V-PDB), δ(C) is relative to atmospheric N₂ (AIR), and Vienna Standard Mean Ocean Water was used for H and O isotope analysis [27,28].

2.4. Multi-Element Mineral Analysis

Solid samples weighing 0.25 g were loaded into the microwave digestion inner tank. Then, 6 mL of nitric acid was added (or 6 mL of nitrous nitric acid and 1 mL of MOS-grade hydrogen peroxide), and the tank was covered with the lid and placed into the microwave digestion apparatus to digest according to certain procedures. After cooling, the sample was removed from the tank, placed on an acid plate, steamed at 150 °C to a small volume (about 0.5 mL), diluted to 10 mL with ultrapure water, well-shaken, and set aside to clarify.

(1)

Nineteen multi-elements (including V, Cr, Co, Ni, Cu, Zn, Sr, Cd, Ba, Tl, Mo, Li, Be, Ga, Ge, Rb, Nb, Cs, and Th) and 16 rare earth elements (including La, Ce, Nd, Y, Sc, Pr, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, and Lu) were determined by ICP-MS;

(2)

ICP-OES detected a total of seven elements of Al, Ca, Fe, K, Mg, Mn, and Na;

(3)

Four milliliters of the solution was dispensed into a colorimetric tube, one milliliter of hydrochloric acid was added, the solution was well-shaken, and Se was finally measured using AFS. Following that, another 4 mL of solution was dispensed into a colorimetric tube, along with hydrochloric acid and thiourea-ascorbic acid solution up to 10 mL; this solution was well-shaken, and As was measured using AFS. The methodological verification is shown in Table S2. For all instrument testing parameters, see Support Information Tables S3–S6.

2.5. Statistical Analysis

All data processing and calculations were performed in WPS Office (11.1.0.12302) (Beijing Kingsoft Office Software, Inc., Beijing, China), OriginPro 9.0 (OriginLab, Northampton, MA, USA), SIMCA 14.1 (MKS Data Analytics Solutions, UMEA, Sweden), and IBM SPSS Statistics 21 software (IBM, Armonk, New York, NY, USA). One-way analysis of variance was used to test the differences in the C, N, H, and O stable isotope ratios and between the 44 elements in the flowers of HBJ at the 5% significance level (p < 0.05). To identify differences between sample groups from different geographic origins for a single variable, a Dunnett-T3 multiple comparison test was performefd on the data. Four stable isotope ratios and 44 multi-elements were analyzed using the supervised chemometric partial least squares discriminant analysis (PLSDA) method for the purpose of tracing origins. One-class partial least squares discriminant analysis (OPLS-DA) was used to distinguish the geographical indication products, namely Tongxiang HBJ, from other producing areas.

3. Results and Discussion

3.1. Stable Isotope Ratio Analysis

The δ¹⁵N, δ²H, δ¹³C, and δ¹⁸O of the HBJ samples from Yancheng (n = 29), Rudong (n = 30), Tongxiang (n = 29), Jinhua (n = 25), Jiangchang (n = 38), and Huangtan (n = 40) were detected and analyzed, as shown in

Table 1. Mean values of N, C, H, and O stable isotopes of HBJ flower from different origins.

Elements	Yancheng (n = 29)	Rudong (n = 30)	Tongxiang (n = 29)	Jinhua (n = 25)	Jiangchang (n = 38)	Huangtan (n = 40)
δ15N	2.71 c	10.93 a	2.87 c	1.11 d	5.06 b	4.98 b
δ13C	−28.26 c	−27.26 b	−26.56 a	−28.40 c	−28.97 d	−28.92 d
δ2H	−78.11 c	−72.26 b	−78.06 c	−85.85 d	−62.94 a	−77.93 c
δ18O	23.25 c	22.26 d	21.02 e	20.83 e	23.97 b	28.40 a

Note: Different letters (a–e) in the same column indicate statistically significant differences (p < 0.05).

and Figure 2. Regardless of the origin, the δ¹⁵N values of all samples ranged from −0.39‰ to 15.87‰. Among the six origins, the Rudong sample has the highest δ¹⁵N value of 10.93‰, followed by Jiangchang (5.06‰), Huangtan (4.98‰), Tongxiang (2.87‰), Yancheng (2.71‰), and Jinhua (1.11‰). As shown in Figure 2A, the δ¹⁵N of Rudong was significantly higher than that of other origins, and there was a significant difference from the other origins (p < 0.05), which proved that the δ¹⁵N value could identify HBJ from different origins. According to previous research [16,29], it is speculated that such differences may be caused by soil fertility, which may also be closely related to artificial fertilization during planting.

Different biochemical pathways of CO₂ fixation in different plants lead to different δ¹³C values. The difference in the δ¹³C values of plants may also be affected by factors such as climate, temperature, sunshine intensity, precipitation, and relative air humidity [30]. Comparing the mean values of δ¹³C in each production area (Table 1 and Figure 2B), it is found that the average value of the six production areas has little difference, ranging from −29.0‰ to −26.6‰, and the order from large to small is Tongxiang (−26.6‰) > Rudong (−27.3‰) > Yancheng (28.3‰) > Jinhua (−28.4‰) > Huang Tan (−28.9‰) > Jiangchang (−29.0‰). It has been confirmed by previous studies that the δ¹³C value of C₃ plants is about −26–−27‰, and the δ¹³C value of C₄ plants is about −12‰ [31]. In this study, Chrysanthemum morifolium Ramat ’Hangabaiju’ is a C₄ plant; however, the δ¹³C in this study was between −26.6‰ and −29.0‰, which is much lower than the previous report. This may be caused by the relatively high CO₂ in these cities [32]. To further understand the influence of climate on δ¹³C values, climate factor information (from the Worldclim database at http://www.worldclim.org/ (accessed on 15 August 2022)) was introduced, as shown in

Table 2. Average information table of climate factors from different origins.

Origin	Yancheng	Rudong	Tongxiang	Jinhua	Jiangchang	Huangtan
Humidity/%	43.5	22.8	38.9	59.9	37.5	37.8
Mean annual temperature/°C	14.2	15.1	16.2	17.1	16.6	16.5
Isothermal (BIO2/BIO7) (×100)	26.1	24.0	24.4	27.7	25.5	25.5
Seasonality of temperature (SD × 100)	913.4	878.3	878.4	815.8	891.1	892.1
Annual precipitation/mm	998.0	1028.8	1143.3	1479.8	1089.0	1086.0
Precipitation seasonality (coefficient of variation)	79.8	55.9	48.0	50.2	52.5	53.4

Note: The data in the table represent the average climatic factors of each origin.

. It was found that the climate information in the table alone had little effect on δ¹³C, which may be caused by the relatively small differences in the latitude, longitude, and altitude of these six producing areas. It is also possible that the difference in the δ¹³C of HBJ may be caused by a combination of multiple climatic factors, rather than only one climatic factor. According to the one-way analysis of variance (p < 0.05) (Table 1), the δ¹³C can be expected to be an effective method for different-origin identification.

According to previous research findings, δ²H and δ¹⁸O in plants reflect the proportion of water absorbed by plants, the evaporation and diffusion effects during transportation, and the biosynthesis in plants. The composition is related to latitude, altitude, distance from evaporation source, temperature, and precipitation [33]. Among the samples from these six origins, the order of the mean value of δ²H is: Jiangchang (−62.9‰) > Rudong (−72.3‰) > Huangtan (−77.9‰) > Yancheng = Tongxiang (−78.1‰) > Jinhua (−85.8‰); the order of the mean value of δ¹⁸O is: Huangtan (28.4‰) > Jiangchang (24.0‰) > Yancheng (23.3‰) > Rudong (22.3‰) > Tongxiang (21.0‰) > Jinhua (20.8‰). Previous studies have confirmed that δ¹³C and δ¹⁵N are positively correlated with longitude and altitude, respectively, and δ¹⁸O and δ²H values are negatively correlated with latitude and altitude, respectively [34,35]. δ²H and δ¹⁸O decreased from low-latitude, low-altitude coastal areas to inland, high-latitude, and increasing altitude [36]. Therefore, it can be seen from Figure 1 and Table 2 that the latitudes of Huangtan and Jiangchang are lower, and the values of δ²H and δ¹⁸O are larger, whereas Yancheng, Rudong, and Tongxiang have higher latitudes, and smaller δ²H and δ¹⁸O ratios. Although the latitude of the Jinhua origin is low, the δ¹⁸O and δ²H ratios are also low, which may be because the Jinhua altitude (161.0 m) is much higher than other areas [35]; therefore, even the annual precipitation (1479.8 mm) is higher and δ¹⁸O and δ²H are lower. The results of one-way variance of δ²H and δ¹⁸O (p < 0.05) showed that the δ²H and δ¹⁸O values of the six production areas were significantly different. This indicated that the values of δ²H and δ¹⁸O could be used as potential isotopic markers to distinguish the different cultivation origins of HBJ.

3.2. Multi-Element Mineral Analysis

Fourty-four multi-elements were determined in 191 samples from six origins, and the average concentrations are listed in

Table 3. Concentration (μg/g) of 44 mineral elements in HBJ samples from six producing areas in China.

Elements	Yancheng (n = 29)	Rudong (n = 30)	Tongxiang (n = 29)	Jinhua (n = 25)	Jiangchang (n = 38)	Huangtan (n = 40)
Se	0.0094 d	0.0095 d	0.0119 d	0.0185 c	0.1108 a	0.057 b
As	0.0297 bc	0.0223 c	0.0168 cd	0.0142 d	0.0608 a	0.0331 b
Li	0.0841 a	0.0561 b	0.0429 b	0.0157 c	0.0947 a	0.0549 b
Be	0.0015 c	0.001 c	0.0014 c	0.0026 b	0.0039 a	0.0019 bc
Ga	0.0172 bc	0.015 bc	0.0146 bc	0.0142 c	0.0325 a	0.0183 b
Ge	0.0014 b	0.0011 bc	0.0008 bc	0.0007 c	0.0026 a	0.0013 b
Rb	10.7397 c	28.079 a	15.749 b	37.86 a	10.8421 c	14.6513 b
Nb	0.0084 b	0.0063 bc	0.0045 c	0.0034 c	0.015 a	0.0089 b
Cs	0.0246 c	0.0849 a	0.0559 a	0.0684 a	0.0425 b	0.0592 a
Th	0.0186 b	0.0114 bc	0.0089 c	0.0092 c	0.025 a	0.0144 b
V	0.1295 bc	0.0988 c	0.0829 cd	0.0532 d	0.317 a	0.171 b
Cr	2.8583	1.5551	1.5277	4.4552	5.9361	2.5439
Co	0.0565 c	0.059 c	0.0829 b	0.1188 a	0.1177 a	0.1295 a
Ni	2.1287 c	1.5683 d	4.9838 a	3.9588 b	1.4105 d	1.3178 d
Cu	12.956 a	12.8797 a	13.1648 a	12.2383 a	9.969 b	8.9077 c
Zn	31.5147 b	29.9754 b	30.2844 b	49.9507 a	22.0355 c	23.4251 c
Sr	16.4219 bc	17.4596 b	15.3155 c	33.7568 a	6.3424 d	6.6427 d
Mo	1.211 b	0.8915 c	0.224 e	0.5245 d	1.2642 a	1.4522 a
Cd	0.1015 f	0.1793 e	0.3632 c	0.9844 a	0.4503 b	0.2258 d
Ba	2.764 c	1.1705 d	3.4722 b	20.2914 a	6.0221 b	5.1263 bc
Tl	0.0023 d	0.0023 d	0.0094 b	0.0423 a	0.006 c	0.0055 c
Al	45.6741 b	36.3351 c	34.0806 c	29.1579 c	105.9441 a	55.2508 b
Ca	6338.9766 b	8795.2547 a	5552.7848 c	4973.884 d	4149.4331 f	4187.0512 e
Fe	98.1835 bc	81.2943 c	88.5503 bc	104.7596 bc	176.6576 a	110.598 b
K	33,350.72 a	32,151.18 b	28,497.8 c	31,450.752 b	23,732.3274 d	24,014.445 d
Mg	3774.3658 b	4466.9468 a	3709.747 b	2830.2142 c	1927.8547 d	1922.862 d
Mn	27.7618 d	40.2097 c	87.8037 b	142.0024 a	26.4587 d	26.3564 d
Na	121.1388 b	177.2628 a	47.4063 c	53.0712 c	27.2487 e	39.783 d
La	0.0524 c	0.0442 c	0.0364 c	0.0773 b	0.0993 a	0.0529 c
Ce	0.0974 bc	0.0828 bc	0.0632 c	0.118 b	0.1862 a	0.0977 b
Nd	0.0412 bc	0.0353 bc	0.0257 c	0.031 c	0.0822 a	0.0431 b
Y	0.027 bc	0.0213 bc	0.0187 c	0.018 c	0.0571 a	0.029 b
Sc	0.0237 b	0.0215 bc	0.0243 b	0.0171 c	0.04 a	0.0208 e
Pr	0.0115 b	0.0099 bc	0.0073 c	0.0095 bc	0.0216 a	0.0113 b
Sm	0.0082 b	0.0068 bc	0.0051 bc	0.0047 c	0.0163 a	0.0084 b
Eu	0.0023 d	0.0017 f	0.0021 e	0.0067 a	0.0052 b	0.003 c
Gd	0.0071 bc	0.0057 bc	0.0047 c	0.0049 c	0.0151 a	0.0076 b
Tb	0.001 b	0.0008 bc	0.0007 c	0.0006 c	0.0021 a	0.001 b
Dy	0.0054 b	0.0044 bc	0.0036 c	0.003 c	0.0115 a	0.0054 b
Ho	0.0012 b	0.0009 bc	0.0007 c	0.0006 c	0.0021 a	0.0014 b
Er	0.0026 b	0.002 bc	0.0016 c	0.0013 c	0.0055 a	0.0026 b
Tm	0.0004 b	0.0003 bc	0.0003 bc	0.0002 c	0.0008 a	0.0004 b
Yb	0.0022 b	0.0018 bc	0.0016 bc	0.001 c	0.0048 a	0.0023 b
Lu	0.0003 bc	0.0003 bc	0.0002 c	0.0002 c	0.0007 a	0.0004 b
Total element	1,271,307.3527	1,375,285.3274	1,105,009.4770	993,712.6204	1,148,039.9990	1,216,859.6091
Macro elements	1,263,970.8340	1,367,719.3280	1,096,424.4080	982,698.0368	1,133,800.8284	1,206,565.6480
Micro elements	5003.7967	4977.8213	6518.8246	7822.7436	8988.1881	6824.2040
Trace elements	2332.1524	2587.6737	2065.8889	3191.4708	5249.4942	3468.9353

Note: The data represent the mean concentration value of 44 elements of HBJ samples. Different letters (a–f) in the same column indicate statistically significant differences (p < 0.05). A total of 44 were divided elements in to macro elements (Ca, Na, K, and Mg), micro elements (Zn, Fe, Cu, Ni, and Mn), trace elements (Se, As, Li, Be, Ga, Ge, Rb, Nb, Cs, Th, V, Cr, Co, Sr, Mo, Cd, Ba, Tl, Al, La, Ce, Nd, Y, Sc, Pr, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, and Lu); the classification is based on the amount of plant growth and development required.

. The results showed that the average content of all elements in the six production areas was in the range of µg/g, and the content of each element in different production areas was different. Among these 44 elements, the content of K was the highest, followed by Ca and Mg, with ranges of 23,732.3274–33,350.72 µg/g, 4149.4331–8795.2547 µg/g, and 1922.862–4466.9468 µg/g, respectively. According to the amount required for plant growth and development, the 44 multi-elements in the samples were divided into macro elements (Ca, Na, K, and Mg), micro elements (Zn, Fe, Cu, Ni, and Mn), and trace elements (Se, As, Li, Be, Ga, Ge, Rb, Nb, Cs, Th, V, Cr, Co, Sr, Mo, Cd, Ba, Tl, Al, La, Ce, Nd, Y, Sc, Pr, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, and Lu); these three sets of results show that the elements of the six origins are all different, and it can be seen that the contribution of multi-elements to the growth of HBJ in different origins is different. The results of one-way variance (p < 0.05) showed that there were significant differences in the 44 elements in six origins, which provided a certain basis for the identification of HBJ origins.

In order to understand the distribution of mineral elements in HBJ samples from different places, heatmaps were used to present the regional visual distribution (see Figure 3). The heatmap data analysis was performed using the OmicStudio tools at https://www.omicstudio.cn/tool (accessed on 3 December 2022). It can be seen from Figure 3 that, compared with Jiangchang and Huangtan, Yancheng, Rudong, Tongxiang, and Jinhua are rich in macro elements (especially Yancheng and Rudong); however, they are less rich in trace elements. This may be because Yancheng, Rudong, Tongxiang, and Jinhua are affected by the sea and saline alkali land along the coast. Among the four HBJ-producing areas of Yancheng, Rudong, Tongxiang, and Jinhua, some micro elements (Zn, Fe, Cu, Ni, and Mn) and trace elements in Jinhua are abnormally higher than those in the other three producing areas, which may be affected by precipitation and altitude (as shown in Figure 1 and Table 2, Jinhua has the highest altitude and the largest annual precipitation). This result is similar to previous results of one-way ANOVA, both of which indicate that there are certain differences in HBJ mineral elements in different places, and that mineral elements can be used for HBJ origin identification.

3.3. Chemometrics Analysis

3.3.1. PLSDA Model Analysis

In order to further identify the origin of HBJ, we established a PLSDA classification model to analyze the δ¹⁵N, δ²H, δ¹³C, and δ¹⁸O stable isotope ratios and the contents of 44 multi-elements in HBJ samples from different origins. PLSDA is a common type of supervised modeling, which is often used in research such as origin traceability and variety identification [37]. In the model, the cumulative contribution of latent variables (R²X, R²Y) and the goodness of fit (Q²) are important parameters for classification, and it is generally believed that these indicators should be greater than 0.5 to indicate that the model is feasible for classification, and the closer these are to 1, the better the classification performance is [16]. The PLSDA 3D score plot (Figure 4A) shows that, although R²X(cum) = 0.923, R²Y(cum) = 0.895, and Q²(cum) = 0.885 in the model, the Jinhua origin is clearly separated from other origins in the model 3D score plots, while Yancheng, Rudong, and Tongxiang are relatively concentrated, and the distinction is not very effective; the same can be found for Huangtan and Jiangchang, which may be due to their relatively close geographical locations (Figure 1). The most important variable (VIP) values are often considered as weights for judging the contribution of the classification model, and the criterion for identifying variables that are important to the model is often a VIP value > 1 [38]. The VIP values of Ni, Cs, Na, H, C, O, Ca, Rb, Mo, Cd, K, Ba, Se, Co, Tl, Mn, Mg, Sr, Cu, Eu, and Zn are all greater than 1, and Ni, Cs, Na, H, and C are considered to be the most important basis for distinguishing HBJ samples of different origins. The receiver operating characteristic curve (ROC) often reflects the relationship between the sensitivity and specificity of an analytical method and is a synthetic representation of the accuracy of the diagnostic method. As shown in Figure 4C, the ROC curve showed that K, Na, and Cu were the variables that had the greatest influence on this method, and the areas under the curve (AUC) were 0.960, 0.855, and 0.840, respectively. Taken together, stable isotopes and mineral elements play an important role in HBJ. The discrimination results of the PLSDA model showed a certain classification effect on HBJ from different origins; however, the classification effect is not very good.

3.3.2. OPLS-DA Model Analysis

OPLS-DA is a classification model developed in PLS-DA. Compared with PLSDA, OPLS-DA has a higher model resolution capability and effectiveness. In order to identify HBJ GI products (Tongxiang City, Zhejiang Province), based on stable isotope ratios and 44 multi-elements, an OPLS-DA model was established to distinguish Tongxiang Hangbaiju from samples from other origins. External validation was used to demonstrate model performance, and VIP values were used to determine isotopic and elemental identification markers. As shown in Figure 5, among the elements with VIP > 1, Mg, Ni, δ¹³C, Mo, and Cd are the main elements that contribute more to the model classification. This indicates that Mg, Ni, δ¹³C, Mo, and Cd elements can be used as the main predictors to distinguish the Tongxiang HBJ samples from other regional samples.

To further verify the distinguishing ability of Tongxiang HBJ products from other products, the external verification method of OPLS-DA was used to evaluate the model. All sample data from the six origins were randomly divided into training and prediction sets in a ratio of approximately 7:3 (Figure 5). As shown in the OPLS-DA score plots (Figure 5), the ordinate represents the within-group gap and the abscissa represents the between-group gap. The results showed that the GI products of HBJ (Tongxiang) could be clearly distinguished from other samples from other origins. This method can be used for the accurate traceability of HBJ.

In order to explore the differences between HBJ GI products (Tongxiang City, Zhejiang Province) and samples from other regions, we divided the samples into a Tongxiang production area group and other production area groups for OPLS-DA to establish a classification model. Likewise, all sample data were randomly divided into training and prediction sets in a ratio of approximately 7:3. The results are shown in Figure 6A; in the OPLS-DA model established based on the δ¹⁵N, δ²H, δ¹³C, and δ¹⁸O stable isotope ratios and multi-element concentrations, the samples from the Tongxiang origin group and other origin groups were completely distinguished. R²X(cum) = 0.848, R²Y(cum) = 0.909, and Q²(cum) = 0.897, which proves that the established model has a good variable explanation and prediction ability. The VIP value of the model (Figure 6B) shows that multiple elements (including Ni, Mo, δ¹³C, Cu, Ce, Nd, Pr, As, Ho, O, Tl, Sm, La, Cd, V, Al, Er, Dy, Y, Tb, Lu, Li, Ga, Yb, Mg, and Tm elements) are meaningful for the identification of HBJ grown in Tongxiang HBJ as seperate from those from other origins (VIP > 1), among which Ni, Mo, δ¹³C, Cu, and Ce elements contribute the most. The results show that a stable isotope ratio combined with multi-element analysis can accurately trace the origin of HBJ from different origins. Moreover, the classification effect was better than the PLSDA model. Therefore, this study provides a new method for the protection of HBJ GI products. This method can also be promoted to trace the origin of more TCMs.

4. Conclusions

The stable isotope ratios of N, C, H, and O, and 44 multi-elements, were used to identify and analyze samples of HBJ from different origins in China, and to further differentiate and identify the origins of GI products. The results of the one-way ANOVA showed that the δ¹⁵N, δ¹³C, δ²H, and δ¹⁸O values and multi-element concentrations could be used to identify the origin of HBJ. PLSDA and OPLS-DA classification models were established to trace the origin of HBJ. The results of the PLSDA model showed a certain identification effect on the origin of HBJ, although the effect was not obvious. However, the OPLS-DA model was established to compare samples from Tongxiang with samples from other origins, and the Tongxiang samples could be distinguished from samples from other origins. In order to further study the differences between the Tongxiang samples and other samples, the samples were divided into a Tongxiang production area group and a group of other production areas for the OPLS-DA to establish a classification model. The results show that many elements (including Ni, Mo, δ¹³C, Cu, Ce, Nd, Pr, As, Ho, δ¹⁸O, Tl, Sm, La, Cd, V, Al, Er, Dy, Y, Tb, Lu, Li, Ga, Yb, Mg, and Tm) were of great significance to the model classification (VIP > 1), among which the Ni, Mo, δ¹³C, Cu, and Ce elements contributed the most, and the Tongxiang origin group and other origin groups were effectively differentiated. Therefore, this study provides a simple, fast, and effective method to protect the interests of consumers and ensure the quality of HBJ’s GI products.