Chemometric authentication of Pu'er teas in terms of multielement stable isotope ratios analysis by EA-IRMS and ICP-MS.

In this work, the stable isotope ratios of carbon, nitrogen, hydrogen, oxygen, and mineral elements and their stoichiometric methods were examined as possible factors that could certify Chinese tea based on its production years. A total of 43 multi-element stable isotope ratios of Xiangzhujing Pu'er tea in five production years were determined through inductively coupled plasma mass spectrometry (ICP-MS) and elemental analyzer-isotope ratio mass spectrometry (EA-IRMS) methods. Two unsupervised learning techniques (principal component analysis and hierarchical clustering analysis) and three supervised learning techniques (partial least squares discriminant analysis [PLS-DA], back-propagation artificial neural network [BP-ANN], and linear discriminant analysis [LDA]) were used on the basis of 18 statistically significant multi-elemental stable isotope ratios to build authentication models for Pu'er tea. The clustering abilities of the two unsupervised learning methods were worse than those of the three supervised learning methods. The three supervised models correctly separated the corresponding production years of the samples. The authentication performance was obtained through BP-ANN and LDA, with 100% recognition and prediction abilities, which were better than those of PLS-DA. δD, δ13C, and 154Sm/152Sm were determined as the markers for the accurate authentication of Pu'er tea in different production years. The profiles of multi-element stable isotope ratios obtained via ICP-MS and EA-IRMS with chemometric methods could serve as potential and powerful factors for authenticating Chinese tea in different production years. This study contributed to the generalization of the use of multi-elemental stable isotope ratio fingerprinting as a promising tool for testing the authenticity of tea worldwide.

1. Introduction

With consumers’ growing interest in health, the demand for information about food authenticity, origin, growth methods, and processing techniques has also increased. For instance, the authentication of tea has been a great deal of interest. Tea is one of the three most popular nonalcoholic beverages worldwide, and it has many health benefits [1,2]. Teas can be aged to ensure the original quality and even improve the quality of tea based on basic tea packaging. For example, the aging of dark tea improves its quality. Properly stored tea does not deteriorate. After aging, the aroma of teas is more fragrant, and their taste is more mellow. Chinese dark tea (CDT) is a post-fermented tea and one of the six major teas in China. Pu’er tea is the most common CDTs that have been accepted and valued by an increasing number of consumers with a healthy lifestyle [3,4]. Given the different production years of Pu’er tea, tea quality within a certain period of time enhances with aging time under appropriate storage conditions [5-8]. Therefore, within a certain period, the older the Pu’er tea, the more expensive the price. Consumers are interested in the production time of Pu’er tea. The production year of teas should be authenticated with different analytical methods.

Most studies on the authentication of tea leaves in different production years have focused on the determination of chemical compositions (e.g., caffeine, catechins, and aroma components) and the examination of morphological characteristics (e.g., shape, size, and color) [9-15]. However, determination and examination methods often lack reproducibility mainly because of variable analytical capabilities, annual changes in the chemical composition of targets, and different environmental or culture conditions of samples [16,17].

In comparison with organic compounds, multiple elements have stable isotope ratios that can also reflect the different growth conditions of plants. Elements are superior to these compounds because they are less affected by processing and storage time, and their contents are stable [18]. Therefore, to accurately authenticate teas in different production years, we proposed an analytical method based on element analyzer-isotope ratio mass spectrometry (EA-IRMS) and inductively coupled plasma mass spectrometry (ICP-MS). In the proposed method, multi-element stable isotope ratios (mainly C, N, H, O, and 87Sr/86Sr) combined with stoichiometric models were used to accurately trace the origin of various teas [18-26]. In general, given the natural differences in physical, chemical, or microbial isotope fractionation processes, the measurement of the stable isotope ratios of C, N, H, and O was considered a promising tool for authenticating various agricultural, such as tea. For example, plant metabolic processes, such as C3, C4, and crassulacean acid metabolism photosynthesis, lead to differences in δ13C. Variations can also occur because of environmental factors, such as water availability, drought stress, nutrient availability, and anthropogenic effects. In plants, the effect of metabolic activity impact on δ13C is greater than that of environmental factors [27,28]. The stable N isotope composition of plants is mainly influenced by soil environment and local nitrogen fertilizer systems, such as type, brand, chemical form, strength, and application timing of different fertilizers [29-32]. Metabolic differences in agricultural can result in different δD and δ18O in sap or plant tissues, and δD and δ18O of such are closely related to climatic conditions, including precipitation, irrigation/surface water sources, height, and distance to the coast [29]. The stable isotope ratios of small mineral elements, such as Sr, are among the most reliable and widely used isotope fingerprints for the determination of tea, food or animal sources [18,33]. However, no study has used stable isotope ratios to authenticate teas in different production years.

Herein, we reported our preliminary investigation on the changes in the 43 multi-elemental stable isotope ratios of C, N, H, O, and various mineral elements in Pu’er tea in different production years (from 2014 to 2018) through EA-IRMS and ICP-MS for the first time. On the basis of multi-elemental stable isotope ratios, we conducted multivariate analyses and applied two unsupervised learning techniques (i.e., principal component analysis [PCA] and hierarchical clustering analysis [HCA]) and three supervised learning techniques (i.e., partial least squares discriminant analysis [PLS-DA], back-propagation artificial neural network [BP-ANN], and linear discriminant analysis [LDA]) to authenticate teas in different production years. We also selected the most important multi-elemental stable isotope ratios as markers and simultaneously compared the authentication performance of the models. We expected that the preliminary results would demonstrate the feasibility of our authentication technology for teas in different production years in China to prevent the false labeling of teas and provide a basis for applying the proposed method for other tea samples worldwide.

2. Materials and methods

2.1. Tea samples

A same brand of Pu’er tea from Xiangzhujing and produced by the same manufacturer in Yunnan Province (Fengqing, Lincang, Yunnan; latitude, 24° 35′; longitude, 100° 04′) in China in the harvesting seasons in 2014, 2015, 2016, 2017, and 2018 (n = 9 for each year) were collected in 2019. These Pu’er tea product samples were stored for 1–5 years. The grades of samples are PINGCHA, and the tea tree variety are all Fengqing large leaf species. Their authenticity, traceability, and equivalent production regimes were ensured by the Yunnan Academy of Agricultural Sciences. The collected tea samples were stored in a freezer at −40°C until multi-elemental stable isotope ratio analysis.

2.2. Sample preparation

Approximately 0.2 g of freeze-dried Pu’er tea samples ground with a ball mill (MM301, Retsch, Germany) and screened with a 100-mesh sieve was accurately weighed (accuracy of 0.0001 g) and then placed in a Teflon digestion vessel. The Teflon digestion vessel was soaked in 20% HNO3 overnight, cleaned with ultrapure water until no sour taste existed, and dried before use. The samples were then treated with a mixture of 5 mL of HNO3 (65% w/w, Merck) and 2 mL of H2O2 (30% w/w, Merck). Heat digestion was conducted in a microwave digestion apparatus (MARS 6 CEM,Matthews, USA) in accordance with the set microwave digestion procedure which digested for 90 min by increasing the power to 1400Wand the temperature to 180°C in a three stepwise fashion. The apparatus was cooled to room temperature to remove the samples after digestion. Subsequently, the lid of the tank in the fume hood was slowly opened to exhaust the gas, and the digestion tank was placed on an acid extraction heater (BHW, Botong, China) at 140°C for acid extraction. Afterward, the digestion tank and the lid were washed with a small amount of ultrapure water 3–4 times, and the washing liquid was combined in a Teflon digestion volumetric bottle at a fixed volume of 50 mL for standby application. The above steps were performed for mineral element stable isotope ratio analysis. The samples were ground to powder, passed through a 100-mesh sieve, and dried in an oven at 65°C for another 48 h for C,N,H,O stable isotope ratio analyses.

2.3. Mineral element stable isotope ratio analysis through ICP-MS

Thirty-nine mineral element stable isotope ratios (109Ag/107Ag,138Ba/137Ba,81Br/79Br,112Cd/111Cd,114 Cd/112Cd,114Cd/111Cd,53Cr/52Cr,72Ge/70Ge,74Ge/72 Ge,74Ge/70Ge,202Hg/200Hg,7Li/6Li,96Mo/95Mo,98 Mo/96Mo,98Mo/95Mo,60Ni/58Ni,207Pb/206Pb,208Pb/207 Pb,208Pb/206Pb,123Sb/121Sb,80Se/78Se,120Sn/118Sn,88 Sr/86Sr,47Ti/46Ti,48Ti/47Ti,48Ti/46Ti,205Tl/203Tl,66Zn/64 Zn,68Zn/66Zn,68Zn/64Zn,153Eu/151Eu,154Sm/152Sm,158 Gd/156Gd,160Gd/158Gd,160Gd/156Gd,164Dy/162Dy,168 Er/166Er,174Yb/172Yb,176Lu/175Lu) analysis of the prepared extracts was performed in a ratio mode by using an ICP-MS instrument (Perkin-Elmer NexION 300X, USA). The internal standards (Re, In, and Rh) and the tuning solution (Be, Ce, Fe, In, Li, Mg, Pb, and U) were used to correct the matrix effects and compensate for the possible variations in instrument performance during determination.

2.4. C,N,H,O stable isotope ratio analysis via EAIRMS

Approximately 4.5 mg (δ13C and δ15N) of Pu’er tea was packed into a 4 mm × 11 mm tin powder capsule, and 1.5 mg (δD and δ18O, balanced for 3 days in the dryer) of Pu’er tea was packed into a 3.5 mm × 5.0 mm silver powder capsule for C,N,H,O stable isotope ratio analysis. Automatic injection was conducted using a flash 2000 element analyzer (EA, ThermoFisher, USA). The teas were loaded into a dynamic fast burning furnace (Prepacke reactor, ThermoFisher, UK), and isotope ratio mass spectrometry (Delta V Advantage, Thermo-Fisher, USA) was introduced. The high-temperature combustion furnace δ13C and δ15N isotopes of the elemental analyzer were set at 980°C, with a carrier gas flow (He) of 230 mL/min. Pyrolysis (δD and δ18O) was conducted at a carrier gas flow rate (He) of 120 mL/min at 1450°C. The isotopic ratio was measured as follows: δsample (‰) = (Rsample/Rstandard – 1) × 1000, where the sample represents δ13C, δ15N, δD, or δ18O, and the R value of the isotope ratio represents 13C/12C, 15N/14N, D/1H, or 18O/16O analytical sample in accortdance with the International Atomic Energy Agency (IAEA, Vienna, Austria) standard. The accuracy determined through reproducibility of the analysis for δ13C, δ15N, δD, and δ18O were 0.15‰, 0.2‰, 1.5‰, and 0.3‰, respectively. δ13C, δ15N, δD, or δ18O (‰) have established reference standards in the international community, that is, Vienna Pee Dee Belemnite for C, N2 (air) for N, and Vienna Standard Mean Ocean Seawater (VSMOW) for H and O.

The main reference materials (IAEA, Vienna, Austria), including IAEA-N1 (NH4SO4, δ15N = 0.4 ± 0.2‰), USGS24 (graphite, δ13C = 16 ± 0.1‰), VSMOW (ocean water, δ2H = 0‰), SLAP (δ2H = − 428.0‰), IAEA601 (δ18Ovsmow = 23.3 ± 0.3‰), and IAEA602 (δ18Ovsmow = 71.4 ± 0.5‰), were used for the multipoint calibration of isotope ratios.

2.5. Statistical analysis

In this study, the 43 multi-elemental stable isotope ratios for each sample were separately analyzed for each analyte. In addition, triplicate measurements were made if the average standard deviation of the duplicates was outside the expected measurement error (±1‰). For multiple comparisons, Tukey’s post hoc test was used, all differences were considered statistically significant if p < 0.05 at 95% confidence interval.

One-way ANOVA, LDA, and BP-ANN were applied using SPSS Statistics version 23.0 (IBM, USA). PCA and PLS-DA were conducted using SIMCA version 13.0 (Umetrics, Umeå, Sweden). Before PCA, LDA, PLS-DA, and BP-ANN were carried out, all variables were “autoscaled.” HCA is used to visualize the correlation between variables and samples. PCA is an unsupervised pattern recognition technique that reduces dimensionality and provides the local view and trend of spatial data [34]. Three supervised learning techniques were applied to establish the authenticity prediction model of Pu’er tea. The dataset with statistically significant variables was assessed through PLS-DA, BP-ANN, and LDA to construct the model. PLS-DA is a regression discrimination tool used to predict a set of variables and identify them as function classes from numerous independent variables [35]. As a supervised pattern recognition technique, BP-ANN has been widely used in stoichiometry. This technique usually consists of three layers, namely, input, hidden, and output layers. The input layer is characterized by rows for classification as an input to the layer; in the output layer, each class in the dataset has an output node [36]. LDA tests the affinity of each sample to the previously defined group by minimizing intragroup variance and maximizing intergroup variance [37].

3. Results and discussion

3.1. Variation in δ13C, δ15N, δD, and δ18O

We initially examined the different δ13C, δ15N, δD, and δ18O profiles over the Pu’er tea in terms of five production years (from 2014 to 2018; Table 1). The distribution and change in δ13C, δ15N, δD, and δ18O in Pu’er tea in different years are shown in Fig. 1. δ13C, δ15N, δD, and δ18O significantly differed in most production years (p < 0.001). The boxes correspond to the interquartile range containing the middle 50% of data, whereas the whiskers indicate the highest and lowest values at 95% and 5% over the entire data range, respectively. The squares inside the boxes represent the mean values, whereas the lines across each box and the filled circles on the box/whisker charts indicate the median and outlier values, respectively. × symbol represents 99% and 1% of the whole data range, and − indicates maximum and minimum values. VPDB, Vienna Pee Dee Belemnite; VSMOW, Vienna Standard Mean Ocean Water. The distribution patterns of these values showed a normal Gaussian distribution even with a few outlier values.

Table 1

Stable isotope ratios and the mean comparison results of the one-way ANOVA of Pu’er tea samples in 2014–2018.

Variables	2014 samples (n = 9)			2015 samples (n = 9)			2016 samples (n = 9)			2017 samples (n = 9)			2018 samples (n = 9)			F	P-Value
	mean ± SD	min	max	mean ± SD	min	max	mean ± SD	min	max	mean ± SD	min	max	mean	min	max
δ13C (‰)	−25.14 ± 0.14a	−25.25	−24.78	−25.06 ± 0.12a	−25.20	−24.85	−25.08 ± 0.15a	−25.32	−24.84	−25.65 ± 0.14b	−25.88	−25.51	−25.69 ± 0.13b	−25.88	−25.48	49.03	0.000
δ15N (‰)	2.01 ± 0.31bc	1.32	2.40	1.76 ± 0.36c	1.27	2.26	2.33 ± 0.52b	1.32	3.03	3.04 ± 0.34a	2.51	3.44	2.93 ± 0.19a	2.51	3.10	21.69	0.000
δD (‰)	−32.37 ± 0.79b	−33.31	−31.34	−37.76 ± 0.82c	−38.92	−36.57	−31.96 ± 0.66b	−32.91	−31.24	−37.04 ± 1.20c	−38.49	−34.60	−26.80 ± 1.21a	−28.99	325.18	190.40	0.000
δ18O (‰)	28.52 ± 0.46b	27.55	29.19	26.32 ± 0.35d	25.79	26.71	27.67 ± 0.41c	27.23	28.47	25.65 ± 0.27e	25.15	26.02	29.13 ± 0.41a	28.29	29.50	127.99	0.000
109Ag/107Ag	0.93 ± 0.14	0.76	1.14	0.99 ± 0.06	0.90	1.07	1.02 ± 0.11	0.87	1.22	0.96 ± 0.13	0.79	1.15	0.98 ± 0.09	0.80	1.07	0.85	0.500
138Ba/137Ba	1.05 ± 0.02ab	1.03	1.09	1.01 ± 0.07b	0.89	1.07	1.02 ± 0.07ab	0.91	1.08	1.06 ± 0.01a	1.04	1.09	0.91 ± 0.01c	0.89	1.09	14.01	0.000
81Br/79Br	0.97 ± 0.20	0.68	1.32	1.02 ± 0.24	0.78	1.41	1.03 ± 0.19	0.59	1.19	0.98 ± 0.14	0.67	1.13	1.02 ± 0.18	0.72	1.23	0.19	0.940
112Cd/111Cd	1.00 ± 0.15	0.68	1.23	0.97 ± 0.07	0.86	1.11	0.98 ± 0.13	0.77	1.22	0.97 ± 0.07	0.86	1.08	0.98 ± 0.14	0.74	1.14	0.14	0.965
114Cd/112Cd	1.08 ± 0.18	0.99	1.53	1.01 ± 0.08	0.91	1.17	0.99 ± 0.08	0.89	1.15	0.99 ± 0.04	0.89	1.03	0.94 ± 0.09	0.82	1.08	2.30	0.075
114Cd/111Cd	1.07 ± 0.10a	0.95	1.27	0.97 ± 0.07ab	0.90	1.08	0.97 ± 0.14ab	0.74	1.19	0.95 ± 0.05b	0.87	1.02	0.92 ± 0.13b	0.76	1.07	2.64	0.048
53Cr/52Cr	0.91 ± 0.03a	0.87	0.96	0.84 ± 0.02c	0.82	0.90	0.85 ± 0.04c	0.82	0.95	0.82 ± 0.01d	0.80	0.84	0.88 ± 0.03b	0.83	0.92	15.94	0.000
72Ge/70Ge	1.02 ± 0.13a	0.78	1.24	0.85 ± 0.18ab	0.62	1.09	0.80 ± 0.23b	0.60	1.33	0.75 ± 0.14b	0.60	0.99	0.87 ± 0.20ab	0.49	1.09	2.82	0.037
74Ge/72Ge	1.55 ± 0.32	0.98	1.99	1.35 ± 0.22	1.06	1.79	1.52 ± 0.22	1.24	1.92	1.49 ± 0.36	1.03	2.15	1.56 ± 0.33	1.09	2.25	0.74	0.572
74Ge/70Ge	1.58 ± 0.40a	1.05	2.17	1.13 ± 0.21b	0.82	1.39	1.20 ± 0.0.33b	0.79	1.85	1.12 ± 0.30b	0.70	1.71	1.33 ± 0.31ab	0.71	1.79	3.35	0.019
202Hg/200Hg	1.00 ± 0.10	0.85	1.11	1.01 ± 0.12	0.87	1.20	1.01 ± 0.10	0.82	1.17	1.01 ± 0.12	0.89	1.24	0.99 ± 0.12	0.87	1.17	0.04	0.997
7Li/6Li	1.04 ± 0.05	0.96	1.12	1.01 ± 0.06	0.90	1.09	1.04 ± 0.06	0.97	1.15	1.05 ± 0.09	0.99	1.23	0.99 ± 0.04	0.93	1.08	1.63	0.186
96Mo/95Mo	1.03 ± 0.09a	0.90	1.18	0.96 ± 0.17ab	0.76	1.35	0.90 ± 0.07bc	0.82	1.01	0.74 ± 0.07d	0.66	0.85	0.82 ± 0.06cd	0.74	0.92	11.25	0.000
98Mo/96Mo	0.95 ± 0.07a	0.84	1.06	0.76 ± 0.15b	0.57	1.11	0.68 ± 0.09bc	0.48	0.81	0.66 ± 0.10c	0.52	0.84	0.70 ± 0.05bc	0.62	0.76	13.18	0.000
98Mo/95Mo	0.98 ± 0.11a	0.82	1.12	0.71 ± 0.08b	0.62	0.84	0.61 ± 0.06c	0.48	0.70	0.49 ± 0.10d	0.38	0.64	0.57 ± 0.04c	0.46	0.61	45.89	0.000
60Ni/58Ni	1.03 ± 0.02	1.00	1.05	1.03 ± 0.02	1.00	1.06	1.03 ± 0.02	1.00	1.06	1.02 ± 0.01	1.00	1.04	1.02 ± 0.02	0.99	1.04	1.78	0.151
207Pb/206Pb	1.00 ± 0.02	0.97	1.03	1.00 ± 0.02	0.96	1.02	1.03 ± 0.02	0.97	1.05	1.02 ± 0.03	0.97	1.06	1.00 ± 0.02	0.96	1.04	2.59	0.051
208Pb/207Pb	1.01 ± 0.03a	0.95	1.03	1.00 ± 0.02ab	0.96	1.03	0.97 ± 0.02c	0.94	1.00	0.98 ± 0.02bc	0.96	1.01	0.98 ± 0.02abc	0.94	1.03	4.57	0.004
208Pb/206Pb	1.01 ± 0.02	0.98	1.03	1.00 ± 0.02	0.97	1.04	0.99 ± 0.01	0.97	1.02	1.00 ± 0.03	0.94	1.04	0.99 ± 0.02	0.96	1.01	1.29	0.292
123Sb/121Sb	1.00 ± 0.04	0.92	1.05	1.02 ± 0.07	0.94	1.16	0.96 ± 0.05	0.89	1.04	0.97 ± 0.08	0.86	1.12	0.95 ± 0.08	0.87	1.11	1.61	0.190
80Se/78Se	1.02 ± 0.04	0.95	1.06	1.02 ± 0.03	0.97	1.06	0.99 ± 0.05	0.93	1.06	0.98 ± 0.03	0.94	1.04	0.99 ± 0.03	0.96	1.05	2.04	0.107
120Sn/118Sn	1.02 ± 0.06	0.93	1.10	0.96 ± 0.13	0.76	1.18	0.98 ± 0.09	0.86	1.08	1.01 ± 0.08	0.88	1.11	1.02 ± 0.11	0.86	1.23	0.75	0.566
88Sr/86Sr	0.87 ± 0.01	0.85	0.88	0.86 ± 0.01	0.85	0.88	0.86 ± 0.01	0.85	0.89	0.87 ± 0.01	0.86	0.88	0.87 ± 0.01	0.85	0.90	1.87	0.135
47Ti/46Ti	0.85 ± 0.10	0.64	0.99	0.95 ± 0.09	0.83	1.08	0.90 ± 0.12	0.76	1.07	1.01 ± 0.17	0.78	1.23	0.96 ± 0.22	0.39	1.08	1.45	0.237
48Ti/47Ti	1.08 ± 0.11a	0.86	1.21	1.05 ± 0.06a	0.94	1.13	1.09 ± 0.20a	0.90	1.58	0.91 ± 0.06b	0.82	0.98	0.87 ± 0.06b	0.78	0.94	7.19	0.000
48Ti/46Ti	0.91 ± 0.06	0.78	1.00	0.99 ± 0.06	0.89	1.04	0.97 ± 0.13	0.81	1.25	0.91 ± 0.12	0.75	1.08	0.84 ± 0.19	0.35	1.01	2.00	0.114
205Tl/203Tl	1.04 ± 0.09	0.93	1.19	1.01 ± 0.12	0.86	1.27	1.07 ± 0.21	0.79	1.47	0.97 ± 0.08	0.86	1.08	0.95 ± 0.11	0.84	1.18	1.29	0.291
66Zn/64Zn	1.01 ± 0.02	0.98	1.03	1.01 ± 0.02	0.98	1.03	1.00 ± 0.02	0.97	1.04	1.00 ± 0.02	0.98	1.05	1.00 ± 0.02	0.98	1.03	0.80	0.534
68Zn/66Zn	0.99 ± 0.02	0.98	1.03	1.00 ± 0.01	0.99	1.02	1.01 ± 0.02	0.98	1.03	1.00 ± 0.02	0.96	1.01	0.99 ± 0.02	0.96	1.03	1.37	0.262
68Zn/64Zn	1.01 ± 0.01	0.98	1.02	1.01 ± 0.01	1.00	1.02	1.00 ± 0.01	0.98	1.02	1.00 ± 0.01	0.98	1.01	0.99 ± 0.01	0.97	1.01	2.48	0.059
153Eu/151Eu	1.30 ± 0.06	1.22	1.42	1.27 ± 0.07	1.19	1.40	1.30 ± 0.09	1.22	1.51	1.25 ± 0.06	1.14	1.33	1.23 ± 0.04	1.18	1.31	1.80	0.149
154Sm/152Sm	1.07 ± 0.03a	1.02	1.12	0.91 ± 0.03b	0.87	0.95	0.85 ± 0.06c	0.79	0.94	0.93 ± 0.05b	0.87	1.01	0.91 ± 0.04b	0.84	0.99	28.32	0.000
158Gd/156Gd	0.60 ± 0.04a	0.54	0.67	0.56 ± 0.04ab	0.51	0.63	0.54 ± 0.04b	0.48	0.60	0.58 ± 0.04a	0.52	0.64	0.60 ± 0.05a	0.55	0.67	3.98	0.008
160Gd/158Gd	0.92 ± 0.03	0.88	0.95	0.90 ± 0.03	0.85	0.94	0.88 ± 0.05	0.79	0.95	0.92 ± 0.04	0.87	0.97	0.92 ± 0.04	0.86	0.99	1.82	0.144
160Gd/156Gd	0.55 ± 0.05a	0.48	0.61	0.51 ± 0.03ab	0.47	0.55	0.47 ± 0.06b	0.38	0.55	0.53 ± 0.05a	0.46	0.58	0.55 ± 0.03a	0.51	0.63	4.89	0.003
164Dy/162Dy	0.94 ± 0.02	0.92	0.97	0.96 ± 0.04	0.88	0.99	0.92 ± 0.04	0.86	0.97	0.94 ± 0.03	0.90	0.97	0.93 ± 0.05	0.88	1.03	1.60	0.193
168Er/166Er	0.96 ± 0.04	0.90	1.02	0.99 ± 0.04	0.93	1.05	0.97 ± 0.06	0.92	1.08	1.00 ± 0.02	0.96	1.02	0.96 ± 0.06	0.90	1.05	1.23	0.312
174Yb/172Yb	0.97 ± 0.04b	0.92	1.05	0.97 ± 0.02b	0.92	0.99	0.98 ± 0.04b	0.92	1.04	1.01 ± 0.03a	0.98	1.06	1.00 ± 0.04ab	0.96	1.06	3.23	0.022
176Lu/175Lu	1.10 ± 0.15	0.83	1.26	1.03 ± 0.12	0.84	1.27	1.09 ± 0.09	0.97	1.30	0.95 ± 0.09	0.82	1.07	1.05 ± 0.16	0.86	1.34	2.10	0.098

Values with different superscripts differ significantly in terms of Pu’er tea in each year (p < 0.05).

Fig. 1

Box/whisker charts of C N H O stable isotope ratios of Pu’er tea in five production years (from 2014 to 2018). The charts show the distribution of all data summarizing the variation in (A) δ13CVPDB, (B) δ15NAIR, (C) δDVSMOW, and (D) δ18OVSMOW in the Pu’er tea samples in terms of the batches of different blend raw materials.

The carbon isotope ratio δ13C low variation ranged from − 25.88‰ to − 24.78‰.~2‰ variation in δ13C in produce grown in the same area due to slight variations in nutrient and water levels available for plant growth [27]. According to previous studies [18-26], δ13C of tea from China varied between −24‰ and −28‰, which is similar to the δ13C (−25.88‰ to −24.78‰) observed for the Pu’er tea in this study. These δ13C distributions also conform to the known δ13C range (−30‰ to −22‰) of C3 plants [27]. Xiangzhujing Pu’er tea in 2015 had lest in negative δ13C mean, and Pu’er tea in 2018 had most in negative δ13C mean (p < 0.001; Fig. 1A and Table 1). The extremely lower δ13C may be related to the increase in artificial CO2 emissions in 2018 [38].

δ15N low variation (~2‰) ranged between 1.27‰ and 3.44‰, with a highest δ15N mean for Pu’er tea in 2017 and a lowest δ15N mean for Pu’er tea in 2015 (p < 0.001; Fig. 1B and Table 1). In general, δ15N of organic fertilizers (1‰–37‰, representatively >5‰) is higher than that of synthetic fertilizers (−4‰–4‰) [39]. Therefore, synthetic fertilizers might be applied to tea leaves in this study. Previous studies [18-26] showed that δ15N of Chinese tea varies between 0‰ and 8‰. These results were similar to our findings; that is, the range of δ15N was 1.27‰–3.44‰. Nitrate levels in irrigation water, N isotope fractionation in soil, and fertilizer availability also affect δ15N of various agricultural [32]. Therefore, low δ15N of Pu’er tea in 2015 predicted the difference in the results of N isotopic fractionation and fertilizer availability, although we did not examine δ15N of organic fertilizers or soil cultivation. As such, the physical, chemical, and microbial properties of soil were affected rather than the specific use of chemical fertilizers by nitrogen groups.

δD of Pu’er tea ranged from −38.92‰ to −25.18‰. Pu’er tea in 2018 had a characteristic profile of δD and the mean δD of Pu’er tea in 2015 and 2017 was lower than that of other Pu’er teas (p < 0.001; Fig. 1C and Table 1). Similarly, Pu’er tea in 2018 had a characteristic profile of δ18O and the mean δ18O of Pu’er teas in 2017 was lower than that of other Pu’er teas (p < 0.001; Fig. 1D and Table 1). Increased precipitation with heavy water isotopes (i.e., or ) results in heavier H and O isotopes in precipitation and groundwater [31]. δD and δ18O in Pu’er tea differed probably because the level was affected by climate in a particular year and topographic features. However, the most important parameter in determining δD and δ18O in Pu’er tea remains unclear. Further studies should be conducted to clarify the comprehensive effects of topography and climate parameters on the variation in δD and δ18O of Pu’er tea in different years.

In conclusion, as δ18O permits to discriminate among the five harvesting years by itself (p < 0.001). From 2014 to 2018 in Lincang, Yunnan, there has a great difference in annual climate, especially annual precipitation each year [40], which cause fractionation of this elements in a particular time frame. It probably could interpret the underlying mechanism of authentication of Pu’er teas using δ18O.

3.2. Variation in mineral element stable isotope ratios

The homogeneity and normal distribution of Pu’er tea were evaluated through one-way ANOVA with univariate variance. The mean and standard deviation of 39 different mineral element stable isotope ratios are summarized in Table 1. Some mineral element stable isotope ratios of Pu’er tea in five production years were significantly different (p < 0.05; Table 1). One-way ANOVA showed that 14 of the 39 mineral element stable isotope ratios in Pu’er tea (138Ba/137Ba, 114Cd/111Cd, 72Ge/70Ge, 74Ge/70Ge, 53Cr/52Cr, 98Mo/96Mo, 96Mo/95Mo, 98Mo/95Mo, 208Pb/207Pb, 48Ti/47Ti, 154Sm/152Sm, 158Gd/156Gd, 160Gd/156Gd, and 174Yb/172Yb) significantly differed in the five production years. As indicated by Duncan’s multiple comparison, Pu’er tea in 2014 had a characteristic mineral element stable isotope ratio fingerprint. In this study, 53Cr/52Cr, 98Mo/96Mo, 98Mo/95Mo, and 154Sm/152Sm ratios were higher in Pu’er tea in 2014 than in other Pu’er tea. Pu’er tea in 2015 had a characteristic mineral element stable isotope ratio fingerprint of 48Ti/47Ti. 114Cd/111Cd were higher in Pu’er tea in 2016 than in other Pu’er tea. Pu’er tea in 2017 had a characteristic mineral element stable isotope ratio profile of 138Ba/137Ba and 174Yb/172Yb.And Pu’er tea in 2018 had a characteristic mineral element stable isotope ratio profile of 158Gd/156Gd and 160Gd/156Gd.This result might be mainly due to the use of different chemical fertilizers from 2014 to 2018, which caused tea to absorb different minerals from soil each year, and further studies should be conducted.

However, any difference could not be reliably distinguished with a single variable. Large differences were observed between the minimum multi-element stable isotope ratios and the maximum multi-element stable isotope ratios in the samples, which differed in the same year. Stoichiometric method that could accurately authenticate tea should be used to evaluate the clustering trend of the samples based on the 18 multi-element stable isotope ratios (p < 0.05) in five production years.

3.3. Exploratory statistical analysis

3.3.1. HCA

The levels of these 18 significant multi-element stable isotope ratios in Pu’er tea in five production years were subjected to HCA and PCA algorithms. HCA showed a cluster analysis tree of 18 multi-element stable isotope ratios (p < 0.05) in the samples (Fig. 2). HCA described three variable clusters. Group 1 contained 208Pb/207Pb, 98Mo/96Mo, 96Mo/95Mo, 98Mo/95Mo, δ13C, and 48Ti/47Ti. The stable isotope ratios of Pb, Mo, and Ti could be in a cluster. Group 2 included δD, δ18O, 72Ge/70Ge, 74Ge/70Ge, and 53Cr/52Cr. The stable isotope ratio of Cr and Ge could be in a cluster. δD and δ18O that were influenced by irrigation water could also be in a cluster. Group 3 comprised 154Sm/152Sm, 158Gd/156Gd, 160Gd/156Gd, 138Ba/137Ba, 114Cd/111Cd, δ15N, and 174Yb/172Yb. The stable isotope ratios of rare earth elements (Yb, Gd, and Sm) could also be in a cluster. In this study, we found that the stable isotope ratios of some mineral elements (Pb, Mo, and Ti; Cr and Ge) were highly correlated to form a cluster. The correlation between the isotope ratios of these elements was attributed to the similarity of soil behavior in the genetic process from parent rock formation to diagenetic mineral weathering and leaching. However, this finding should be further verified. The stable isotope ratios that formed a cluster behaved similarly to Pu’er tea.

Fig. 2

Hierarchical clustering analysis (HCA) between variables and samples. Correlation coefficient is indicated by the intensity of colors as shown by the color scale. (2014, 2015, 2016, 2017, and 2018 indicate Pu’er tea in 2014, 2015′, 2016′, 2017′, and 2018′, respectively).

HCA revealed six sample clusters. Years 2014, 2015, 2017, and 2018 formed clear clusters, while 2016 is divided into two clusters (one individual and second as a part of 2017). HCA was conducted to accurately authenticate Pu’er tea in different production years by choosing 18 multi-element stable isotope ratios, which could not be specified.

3.3.2. PCA

PCA reduces the number of variables used in data description and is the most commonly used method for calculating components, such as linear potential variables. The unsupervised principal component variance contribution rate and the load of 18 principal component variables used for the certification of Pu’er tea in the five studied years are shown in Fig. 3. In this study, the first three principal components (PCs) extracted in accordance with the Kaiser criterion represented PC1 (29.16%), PC2 (18.55%), and PC3 (13.38%; Fig. 3). The load-scattering diagrams of 18 variables showed that PC1 was strongly influenced by δ15N and 174Yb/172Yb. The key components of PC2 included 160Gd/156Gd and 158Gd/156Gd. PC3 was basically contributed by δ18O and δD (Fig. 3B, D).

Fig. 3

PCA plots of multi-element stable isotope ratios showing the PC1, PC2, and PC3 effects to authenticate Pu’er tea in five production years (from 2014 to 2018). The ellipse on the score plots represents the 95% confidence region for Hotelling’s T2. Score plot of PC1 and PC2 (A), loading plot of PC1 and PC2 from PCA (B), Score plot of PC1, PC2 and PC3 (C), loading plot of PC1, PC2 and PC3 from PCA (D).

PC1 and PC2 were the main components. The Pu’er tea could be preliminarily divided into five classes corresponding to 2014 to 2018 from the score plots for PC2 versus PC1 (Fig. 3A). Pu’er tea in 2014 showed negative scores in PC1 (the second and third quadrants), and this finding was easily differentiated from Pu’er tea in other years. A minimal overlap was found between Pu’er tea in 2015 and 2016, resulting in negative scores in PC2 (the third and fourth quadrants). Pu’er tea in 2017 with positive scores in PC1 and PC2 (the first quadrant) and Pu’er tea in 2018 with positive scores in PC1 (the first and fourth quadrants) could be easily differentiated from Pu’er tea in other years. The score plot for PC1 versus PC2 versus PC3 could show an improved differentiation among the scores corresponding to the different Pu’er tea in the five studied production years (Fig. 3C). However, different production years partially overlapped, and the PCA classification was not strong. In summary, the fractional scattering figure of Pu’er tea could be preliminarily divided into five classes based on the five studied production years under the orthogonal coordinate system of PC1, PC2, and PC3 on the basis of their 18 multi-element stable isotope ratio fingerprints.

The application of these two unsupervised methods of pattern recognition (HCA and PCA) revealed the natural grouping of the Pu’er tea in the five studied production years in the original data matrix, indicating a tendency to group samples. Remarkably, HCA and PCA showed that Pu’er tea in the five studied production years slightly overlapped. Although HCA and PCA could provide a visual picture of how the samples were clustered, they could not present information about the quality of clustering and confidence in such clustering. Therefore, other supervised stoichiometric tools were further explored to authenticate Pu’er tea.

3.4. Classification and predictive modeling for sample authentication

PLS-DA, BP-ANN, and LDA models were established on the basis of the determination of 18 element stable isotope ratios in five Pu’er tea in different production years (from 2014 to 2018; p < 0.05) to use the multi-element stable isotope ratio for the authentication of these.

3.4.1. PLS-DA

PLS-DA is a projection method that maximizes the separation of observation groups by rotating PCA, which was an appropriate way to authenticate Pu’er tea in the five studied production years (Fig. 4). The PLS-DA score plot (Fig. 4A) was clear clustering, which revealed that the two highest-ranking R2X accounted for 47.1% of the total variance. The first R2X, accounting for 29.0% of this variance, separated the multi-element stable isotope ratio feature of the Pu’er tea in the five studied production years, which was 13.3% of this variance in the second R2X. In PLS R2X 1, the largest contributor was the highest stable isotope ratio of Mo in Pu’er tea in 2014 and indicated as eigenvectors of 0.418, 0.339, and 0.316 for 98Mo/95Mo, 98Mo/96Mo, and 96Mo/95Mo, respectively. In R2X 2, the largest contributors were δ13C, 138Ba/137Ba, and 48Ti/47Ti. The eigenvectors of δ13C, 138Ba/137Ba, and 48Ti/47Ti were 0.340, 0.304, and 0.227, respectively (Fig. 4B). The variable importance in projection (VIP) value explains the contribution of variables in the projection; a VIP value of >1 is usually used to identify variables essential for modeling [41]. Six variables, namely, 138Ba/137Ba (VIP: 1.365), δD (VIP: 1.353), δ18O (VIP:1.264), δ13C (VIP: 1.200), δ15N (VIP: 1.146), and 154Sm/152Sm (VIP: 1.133), were observed as significant multi-element stable isotope ratio features for the authentication of Pu’er tea (Fig. 4C).

Fig. 4

Comparison of PLS-DA results derived from 18 multi-element stable isotope ratios of Pu’er tea in five production years (from 2014 to 2018). Score plot (A), loading plot from PLS-DA (B), and variable importance in projection (VIP) values obtained from the PLS-DA model (C). The ellipse on the score plots represents 95% confidence region for Hotelling’s T2. (1, 2, 3, 4, and 5 indicate Pu’er tea in 2014, 2015, 2016, 2017, and 2018, respectively).

Validation is critical to ensure the reliability of the developed PLS model because of the overfitting of PLS to the model, and data usually result in a split class score [42]. Although no threshold is set to compare criteria or determine importance, quality assessment statistics (Q2) is usually the result of cross-validation [43]. In general, Q2 greater than a 0.5 threshold indicates good predictability, whereas poor Q2 corresponds to noisy data [44]. Our PLS-DA model displayed a good performance (R2X = 0.765, R2Y = 0.843, and Q2 = 0.666) and explained 76.5% and 84.3% of the variations in X and Y, respectively, with a predictive Q2 ability of 66.6%.

To our knowledge, the applicability of PLS-DA for tea authentication has never been reported. At present, multi-element stable isotope ratio analysis combined with PLS-DA provides a reliable Pu’er tea authentication in China. The PLS-DA score plot showed a cluster from Pu’er tea in different production years. However, a small overlap, which also existed in PCA, was observed between Pu’er tea in 2015 and 2016. This result could be attributed to their computational principles. PCA and PLS-DA can be applied to form components that capture most information in explanatory variables by maximizing the variance [35].

3.4.2. BP-ANN

A BP-ANN model was studied to further improve the authentication accuracy. In this study, a three-layer BP-ANN was established and applied to authenticate the Pu’er tea in different production years. The input and output layers had 18 and 5 neurons, respectively. The number of hidden layers was 1, and the number of neurons was 10. The recognition and prediction abilities of the model were evaluated.

In this study, 70% of the samples were randomly selected from Pu’er tea as a training set and cross-validation to evaluate the recognition ability. The remaining 30% of the samples were selected as the test set for the external verification test to evaluate the prediction ability. In the model training process, the five groups of Pu’er tea representing the corresponding production years could be successfully divided. In Table 2, the recognition ability of overall accuracy was 100.0% in all Pu’er tea. The remaining 30% of Pu’er tea were classified using the established model for external verification to further test the prediction ability of this model. Table 2 shows that 30% of the samples could be correctly predicted into the five groups. The overall prediction accuracy was 100.0%, indicating the good applicability of the proposed BP-ANN model.

Table 2

Model training and prediction results of the BP-ANN model.

	Pu’er tea in 2014	Pu’er tea in 2015	Pu’er tea in 2016	Pu’er tea in 2017	Pu’er tea in 2018	accuracy (%)
Model training
Pu’er tea in 2014	5	0	0	0	0	100.0
Pu’er tea in 2015	0	6	0	0	0	100.0
Pu’er tea in 2016	0	0	7	0	0	100.0
Pu’er tea in 2017	0	0	0	7	0	100.0
Pu’er tea in 2018	0	0	0	0	3	100.0
Recognition ability (%)						100.0
Cross-validation
Pu’er tea in 2014	5	0	0	0	0	100.0
Pu’er tea in 2015	0	6	0	0	0	100.0
Pu’er tea in 2016	0	0	7	0	0	100.0
Pu’er tea in 2017	0	0	0	7	0	100.0
Pu’er tea in 2018	0	0	0	0	3	100.0
Recognition ability (%)						100.0
Model test sets
Pu’er tea in 2014	4	0.0	0.0	0.0	0	100.0
Pu’er tea in 2015	0.0	3	0.0	0	0	100.0
Pu’er tea in 2016	0.0	0.0	2	0.0	0	100.0
Pu’er tea in 2017	0.0	0.0	0	2	0	100.0
Pu’er tea in 2018	0	0	0	0	6	100.0
prediction ability (%)						100.0

In terms of the importance and standardization importance of variables in the BP-ANN (Table S1), the six most important variables as markers in the BP-ANN model were as follows: δ18O (importance: 0.101, standardization importance: 100%), δD (importance: 0.088, standardization importance: 87.40%), 154Sm/152Sm (importance: 0.082, standardization importance: 81.70%), δ13C (importance: 0.077, standardization importance: 76.70%), 138Ba/137Ba (importance: 0.073, standardization importance: 72.30%), and 48Ti/47Ti (importance: 0.068, standardization importance: 67.30%). The performance of the BP-ANN model was more accurate than that of the PLS-DA method.

3.4.3. LDA

The LDA model was also used to further improve the authentication accuracy. The recognition and prediction abilities of the LDA model were also evaluated. The Pu’er tea were divided into training sets, cross-validation (70% of samples) and test sets (30% of samples) for external validation. The recognition capability of the model was examined with the training set, whereas the prediction ability was verified through cross-validation and external validation. The results showed that two statistically significant discriminant Fisher functions were formed: Wilks’ lambda = 0.000, X2 = 207.895, df = 20, p < 0.001 for the first Fisher function, and Wilks’ lambda = 0.015, X2 = 108.972, df = 12, p < 0.001 for the second Fisher function. A significant Wilks’ lambda indicated that the discriminant function was the basis for differentiating populations. The testing of the uniformity of variability (BoxMindex = 111.636, F = 1.039, p = 0.396) was not significant at 95% confidence level, suggesting that the variability of the samples in each production year was consistent. The first discriminant Fisher function accounted for 77.4% of the total variance, and the second discriminant Fisher function accounted for 14.2%. Both accounted for 91.6% of the total variance, showing a high level. After forward stepwise variable selection based on Wilk’s lambda criterion, δD, δ13C, δ15N, 98Mo/95Mo, and 154Sm/152Sm were five variables that factored into the simplified model.

Discriminant analysis showed that the Pu’er tea in different production years were well authenticated (Fig. 5). The first discriminant function showed an evident authentication function for Pu’er tea in 2015 and 2018 compared with that in the other years. The second discriminant function had an evident identification function for Pu’er tea in 2014 and 2017 compared with that in the other years. Recognition ability was expressed as the percentage of Pu’er tea in different production years correctly classified during model training, with an overall accuracy of 100.0%, indicating an excellent result. The prediction ability was expressed as the percentage of Pu’er tea in different production years correctly classified using a typical cross-validation procedure and an external validation procedure. The overall prediction accuracy was 100.0% for cross-validation and external validation methods. This result suggested a satisfactory value, especially for this method (Table S2). The multi-element stable isotope ratios (independent variables) constituting the first and second discriminant functions are shown in Table S3. δD strongly influenced the first discriminant function. This finding showed that Pu’er tea in 2015 was the most positive, whereas the four other stable element isotope ratios mainly contributed to the second discriminant functions.

Fig. 5

Authentication of Pu’er tea in five production years (from 2014 to 2018) based on the 18 multi-elemental stable isotope ratios (p < 0.05).

The performance of separation between the groups in BP-ANN and LDA was better than that of PLS-DA possibly because of different training algorithms. BP-ANN training involves multiple perception levels, whereas the goal of LDA is to project high-dimensional data into a low-dimensional space and obtain the maximum class discrimination through the ratio of the maximum interclass and intraclass distances [36,37]. Among the variables selected through different methods, the three models shared δD, δ13C, and 154Sm/152Sm as markers of Pu’er tea in different production years.

4. Conclusions

The fingerprints of Pu’er tea in five production years (from 2014 to 2018) consisted of 43 multi-elemental stable isotope ratios and were obtained through EA-IRMS and ICP-MS. These techniques are reliable tools for generating chemical information-rich multi-elemental stable isotope ratio fingerprints. This study confirmed that 18 statistically significant multi-elemental stable isotope ratio signatures could be applied as fingerprints to authenticate Pu’er tea. This study also demonstrated significant differences in some of the multi-elemental stable isotope ratios and Pu’er tea in different production years by projecting the multi-elemental stable isotope ratio data on the computing platform of two unsupervised and three supervised learning techniques. The clustering ability of the two unsupervised learning methods were worse than that of the three supervised learning methods. The supervised models proposed using the LDA, PLS-DA, and BP-ANN algorithms showed that the authentication performance of BP-ANN and LDA was better than that of PLS-DA, with a recognition ability of 100% and a prediction ability of 100%. δD, δ13C, and 154Sm/152Sm were the markers for enabling the accurate authentication of Pu’er tea in different production years.

In conclusion, we developed a feasible method for multi-elemental stable isotope ratio analysis and improved the process of identifying the authenticity of tea by using a stoichiometric model. Our study could enhance the authentication procedures for existing tea and facilitate the detection of fake brand labels. The established prediction model had the advantages of high accuracy, strong robustness, human factor independence, and so on. As such, it had a certain application prospect in food classification and group distribution, especially in food source traceability and adulteration research. This work presented the great potential of multi-elemental stable isotope ratio fingerprinting for evaluating the authenticity of foods, such as tea. However, our results must be interpreted cautiously, as a limited number of Pu’er tea samples were tested. Future studies should evaluate a larger number of Pu’er tea samples from different years to establish global isotope markers for reliable authentication on a wider scale.