Assessing Seasonal Effects on Identification of Cultivation Methods of Short–Growth Cycle Brassica chinensis L. Using IRMS and NIRS

Seasonal (temporal) variations can influence the δ13C, δ2H, δ18O, and δ15N values and nutrient composition of organic (ORG), green (GRE), and conventional (CON) vegetables with a short growth cycle. Stable isotope ratio mass spectrometry (IRMS) and near-infrared spectroscopy (NIRS) combined with the partial least squares-discriminant analysis (PLS-DA) method were used to investigate seasonal effects on the identification of ORG, GRE, and CON Brassica chinensis L. samples (BCs). The results showed that δ15N values had significant differences among the three cultivation methods and that δ13C, δ2H, and δ18O values were significantly higher in winter and spring and lower in summer. The NIR spectra were relatively clustered across seasons. Neither IRMS-PLS-DA nor NIRS-PLS-DA could effectively identify all BC cultivation methods due to seasonal effects, while IRMS-NIRS-PLS-DA combined with Norris smoothing and derivative pretreatment had better predictive abilities, with an 89.80% accuracy for ORG and BCs, 88.89% for ORG and GRE BCs, and 75.00% for GRE and CON BCs. The IRMS-NIRS-PLS-DA provided an effective and robust method to identify BC cultivation methods, integrating multi-seasonal differences.


Introduction
Organic agriculture is based on environmentally friendly, product-safe, sustainable, and comprehensive agricultural practices, and is becoming a popular choice for global agricultural development.Furthermore, increasing awareness of health and environmental concerns among consumers is prompting them to pay more attention to organic products.Global organic market sales reached 134.8 billion euros in 2022 [1].Countries actively promote organic farming methods while adopting relatively neutral and market-oriented policies based on their own circumstances.For example, China was ranked the thirdhighest producer of organic products in 2022, with a global market share of 9.2%.China's government has also proposed to improve green agricultural standard and systems by strengthening the certification management of green food, organic products, and agricultural products with geographical indication status [2], sending strong signals that green agriculture is an important form of high-quality and sustainable development for Chinese Foods 2024, 13, 1165 2 of 15 agriculture.Green agriculture, defined by China National Standard GB/T 33761-2017 [3] and Agriculture Industry Standard NY/T 391-2021 [4], is a blend of organic and conventional agriculture, which primarily reduces the need for chemical fertilizers and pesticides.
As the societal demand for sustainable agriculture grows, the authenticity of organic and green foods has become an urgent issue due to profit-driven practices such as false advertising and misleading labels [5].Several effective analytical methods have been used to ensure their authenticity, primarily including stable isotope analysis [6], spectral analysis [7], and chemical analysis [8].In these methods, stable isotopes (δ 13 C, δ 15 N, δ 2 H, and δ 18 O) can objectively reflect climatic conditions (temperature, humidity, light intensity, and precipitation), soil composition, and agricultural input (fertilizers, pesticides, etc.) information for plant agricultural products [9].For instance, δ 13 C values are indicative of photosynthesis pathways (C3, C4, or CAM) utilized by plants but are also influenced by light intensity, humidity, and environmental pollutants (car emissions) [9].δ 15 N values offer information about plant nutrient sources and agricultural fertilization practices.Organic fertilizers (such as animal-derived manures and plant composts formed of non-leguminous products) typically exhibit relatively high δ 15 N values, distinguishing organic cultivation methods from green or conventional farming practices [9][10][11].δ 2 H and δ 18 O values are used to identify different irrigation water sources and are influenced by rainfall, local temperatures, latitude, altitude, and distance from the sea [9,12].Moreover, near infrared spectroscopy (NIRS) in spectral analysis mainly captures essential vibrational and rotational stretching details related to the hydrogen bonds (C-H, N-H, O-H, and S-H) of nutritional compositions of agricultural products and offers chemical-free, rapid, and non-destructive advantages in analyzing the composition of agricultural products [13].
Both stable isotope analysis and NIRS have been employed as powerful tools for the identification of organic, green, and conventional agricultural products based on their growing fertilizer types, climate environment, and nutritional component variations arising from different cultivation methods [11,14].However, plant isotopes and NIR spectra can be altered by interannual climatic variations and seasonal effects, potentially affecting the accuracy of identifying cultivation methods [15,16].This is particularly pronounced for crops with short growth or maturation periods (<60 days).Short-growth cycle leafy vegetables are more sensitive to changes in precipitation or temperature, causing variability in their chemical compositions and absorption characteristics of various elements [17].However, seasonal effects on the stable isotopes and NIR spectra of short-growth cycle leafy vegetables need to be further investigated, along with the identification of shortgrowth cycle leafy vegetables cultivation methods using stable isotopes and NIRS.
Therefore, this study aimed to investigate the influence of seasonal variations on the δ 13 C, δ 15 N, δ 2 H, and δ 18 O values and NIR spectra of short-growth cycle and year-round planted Brassica chinensis L. (BC) under various cultivation methods.Furthermore, IRMS or/and NIRS data combined with chemometrics were used to ensure the authenticity of BC cultivation methods despite seasonal variations.

Sample Collection and Preparation
In this study, 175 BC samples comprising 63 organic (ORG) samples (defined by China National Standard GB/T 19630-2019) [18], 44 green (GRE) samples, and 68 conventional (CON) samples were collected between September 2020 and September 2021 from Shanghai vegetable farms.Samples collected from September 2020 to November 2020 were classified as autumn (mean temperature 19.6 • C), from December 2020 to February 2021 as winter (mean temperature 7.4 • C), from March 2021 to May 2021 as spring (mean temperature 16.9 • C), from June 2021 to August 2021 as summer (mean temperature 27.5 • C), and from September 2021 as autumn-repeat (autumn-re) (mean temperature 26.7 • C) (Table 1).Autumn-re samples were used as interannual samples to verify seasonal effects and identification models.BC samples were prepared as reported in our previous paper [16].Briefly, about 2.5 kg of each vegetable was collected, rinsed with deionized water to remove soil or dust, pulped using a blender, frozen at −18 • C for 6 h, and freeze-dried at −54 • C for at least 72 h.The dried samples were ground into a fine and uniform powder with a particle size less than 0.15 mm.The powder was stored in a desiccator to prevent the absorption of local atmospheric water and preserve the weakly reflected exchangeable and non-exchangeable δ 2 H and δ 18 O signatures from the different cultivation methods and seasonal effects.

Stable Isotope Analysis
The δ 13 C, δ 15 N, δ 2 H, and δ 18 O values of samples were determined using a Flash IRMS elemental analyzer (EA) interfaced to a DELTA V Advantage isotope ratio mass spectrometry system (IRMS, Thermo Fisher Scientific Inc., Bremen, Germany) using similar methods outlined in Liu et al. [19].In C/N mode, the oxidation and reduction furnace temperatures of the EA were set at 980 • C. High purity helium was used as the carrier gas with a flow rate of 180 mL/min.High purity CO 2 and N 2 were used as the reference gases with flow rates of 60 mL/min.About 1.6 mg samples were weighed into tin capsules for δ 13 C and δ 15 N analysis with an 80% dilution ratio of CO 2 produced by these samples during analysis.In H/O mode, the EA pyrolysis temperature was set at 1380 • C. High purity helium was used as the carrier gas with a flow rate of 100 mL/min.High purity CO and H 2 were used as the reference gases, also at a flow rate of 100 mL/min.About 0.3 mg samples were weighed into silver capsules for δ 2 H and δ 18 O analysis with a 40% dilution ratio of H 2 and a 60% dilution ratio of CO produced by these samples during analysis.Isotope ratios were calculated using the following Equation (1): where X represents δ 13 C, δ 15 N, δ 2 H, or δ 18 O; R sample denotes the abundance ratio of heavy isotope against light isotope, e.g., 13 C/ 12 C, 15 N/ 14

NIRS Analysis
NIR spectra were collected using a Nicolet iS50 Fourier transform near-infrared spectrometer (Thermo Fisher Scientific, Waltham, MA, USA) with an integrating sphere mode and an InGaAs detector.The spectral range was from 10,000 cm −1 to 4000 cm −1 , spectral resolution was set at 8 cm −1 , scan time was 32, and an internal blank was used as the reference for the measurements.The powdered samples were thoroughly mixed before each scan and then placed in a rotating sample cup and scanned three times.NIR spectral data were captured using OMNIC 9 software and stored in absorbance format.An averaged spectrum, generated from the three replicate analyses, contained 1557 variables and was used for both the calibration and validation sets.The laboratory temperature was maintained at a constant 25 • C throughout the analysis period.

Statistical Analysis and Chemometrics Methods
One-way analysis of variance (ANOVA) was applied to evaluate and compare differences in the δ 13 C, δ 15 N, δ 2 H, and δ 18 O values of BCs attributed to seasonal effects using Matlab R2020a software (MathWorks, Natick, MA, USA).Boxplots created in Microsoft Office Excel 365 (Microsoft, Redmond, WA, USA) were employed to visually represent the differences between the four isotopes across different BC cultivation methods or seasons.
The spectral data underwent Norris smoothing and derivative (NSD) treatment prior to modeling, which aimed to reduce or eliminate random baseline shifts, light scattering, and noise interferences, ensuring that only useful information was incorporated into the spectral signal [20].NSD pretreatment consists of smoothing involving parameters 's' and 'g', where 's' represents the number of data in one segment and 'g' is the number of data in one gap, and a derivative containing the first or second derivative.Before building the models, the Kennard-Stone (KS) algorithm was used to divide the samples into calibration set (75%) and validation set (25%) [21].That is, the samples were selected one by one based on the furthest distance from each other using Euclidean distance, thus dispersing them across the multivariate space.Partial least squares-discriminant analysis (PLS-DA) [22] was utilized to identify the BC cultivation methods using Matlab R2020a software and SIMCA 14.1 software (Umetrics, Umeå, Sweden).K-fold cross validation (CV) was applied to determine the optimal number of latent variables (optLVs).The optLVs were used to build a calibration model of PLS-DA by sensitivity (SE), specificity (SP), area under curve (AUC), and classification accuracy for evaluation (see Supplementary Materials S1 for the equations of SE, SP, and classification accuracy) [23], and then a validation model was employed to predict the remaining samples (25%).The predictive ability of the model was assessed by an accuracy.

Overall and Seasonal Isotopes of Different BC Cultivation Methods
There were differences in the BC stable isotopes determined for different cultivation methods (Figure 1).ORG BC had the highest overall mean δ 15 N (10.50 ± 6.20‰) and δ 18 O (21.42 ± 2.85‰) values, coupled with the lowest overall mean δ 13 C (−29.17 also found three CON BC samples with δ 15 N values exceeding 20‰ (Figure 1, blue possibly due to the fact that these samples were from farms transitioning from co tional to organic cultivation.Slight δ 13 C, δ 2 H, and δ 18 O variations among BC culti methods were most likely due to the different fertilizer types and seasonal effects perature, light intensity, and precipitation), resulting in differences in carbon cycling tosynthetic efficiency, and water use efficiency [9,26,27].2).There were no significant differences in seasonal mea and δ 2 H values among the three cultivation methods of BC, possibly indicating th photosynthetic efficiency and water use efficiency of BC under different cultivation ods were similar in each season, given their similar growing locality in Shanghai.Se mean δ 15 N values of BC also varied similarly among the three cultivation methods, f ing the sequence ORG > GRE > CON.Significant differences in δ 15 N values were ob between ORG and CON BCs during winter, summer, and autumn-re (p < 0.05).CO exhibited significant differences compared to ORG and GRE BCs (p < 0.05) in sprin significant differences were observed in autumn, possibly due to fertilizer type var and the two outlier values (20.20‰ and 22.53‰) of CON BCs in autumn.ORG BC the highest seasonal mean δ 18 O value (23.37 ± 2.61‰) in spring and were significant ferent from CON BCs (21.65 ± 1.50‰) (p < 0.05).This possibly arose from the enhanc of soil permeability and water retention caused by organic fertilizer application a Commonly, the organic fertilizers used in the cultivation of ORG and GRE BCs are animal manures or plant composts which undergo denitrification, promoting the volatilization of the light stable isotope 14 N fraction, enhancing residual 15 N in the fertilizer [24] and resulting in higher ORG and GRE BCs δ 15 N values compared to CON BC.However, if the storage and fermentation time of organic fertilizer, especially for animal manures, is short and denitrification has not occurred, or if the organic fertilizer is plant-based (legume), the δ 15 N values of organic vegetables may not be so high [25].This effect might explain the low δ 15 N values (around 0‰) of some ORG BC samples (Figure 1).The study also found three CON BC samples with δ 15 N values exceeding 20‰ (Figure 1, blue circle), possibly due to the fact that these samples were from farms transitioning from conventional to organic cultivation.Slight δ 13 C, δ 2 H, and δ 18 O variations among BC cultivation methods were most likely due to the different fertilizer types and seasonal effects (temperature, light intensity, and precipitation), resulting in differences in carbon cycling, photosynthetic efficiency, and water use efficiency [9,26,27].
The δ 13 C, δ 15 N, δ 2 H, and δ 18 O values among ORG, GRE, and CON BC varied across different seasons (Table 2).There were no significant differences in seasonal mean δ 13 C and δ 2 H values among the three cultivation methods of BC, possibly indicating that the photosynthetic efficiency and water use efficiency of BC under different cultivation methods were similar in each season, given their similar growing locality in Shanghai.Seasonal mean δ 15 N values of BC also varied similarly among the three cultivation methods, following the sequence ORG > GRE > CON.Significant differences in δ 15 N values were observed between ORG and CON BCs during winter, summer, and autumn-re (p < 0.05).CON BCs exhibited significant differences compared to ORG and GRE BCs (p < 0.05) in spring.No significant differences were observed in autumn, possibly due to fertilizer type variations and the two outlier values (20.20‰ and 22.53‰) of CON BCs in autumn.ORG BCs had the highest seasonal mean δ 18 O value (23.37 ± 2.61‰) in spring and were significantly different from CON BCs (21.65 ± 1.50‰) (p < 0.05).This possibly arose from the enhancement of soil permeability and water retention caused by organic fertilizer application and the improvement of metabolic activity due to higher daily temperatures in late spring, leading to more positive 18 O enrichment in BC tissue from H 2 18 O and C 18 O 2 [9,27,28].The results confirmed that the δ 15 N values, mainly influenced by fertilizer type, characterized the three cultivation methods of BC.There were significant differences in some δ 15 N and δ 18 O values among different BC cultivation methods in a single season, indicating that different fertilizers were used for BCs grown under the same cultivation method.Moreover, BCs were also influenced by seasonal factors, such as temperature, light intensity, and precipitation, and required further investigation.

Seasonal Isotopes for Each BC Cultivation Method
The stable isotopes of different BC cultivation methods showed distinct variations based on seasonal time series (Figure 2).Mean δ 13 C values for each season of organic (δ 13 C Organic ), green (δ 13 C Green ), and conventional (δ 13 C Conventional ) BCs were higher in winter.A significant difference between winter and summer (p < 0.05) was noted, due to seasonal variations in temperature and light intensity influencing photosynthetic efficiency [9,27].In summer, BCs exhibited vigorous photosynthesis, preferentially absorbing a higher proportion of lighter 12 CO 2 from the atmosphere, leading to lower δ 13 C values in BC tissue.Mean δ 15 N values for different BC cultivation methods showed no significant seasonal differences, as each cultivation method had specific fertilizer treatments.Each seasonal mean δ 2 H and δ 18 O value of different BC cultivation methods exhibited a similar trend across seasons (Figure 2).In winter and spring, the three cultivation methods of BCs had more positive mean δ 2 H and δ 18 O values compared to summer and autumn.There were significant δ 2 H and δ 18 O value differences between winter and/or spring and autumn, as well as between summer and autumn-re (p < 0.05).This trend was contrary to the usual expectations for precipitation, where δ 2 H and δ 18 O values are more positive in summer and more negative in winter [12,29].However, it was consistent with the seasonal δ 2 H and δ 18 O variation in local atmospheric precipitation and the irrigation water source from the Yangtze River [30], which are significantly influenced by the distinct monsoon system and local topography around Shanghai [31][32][33].In summer, precipitation in Shanghai, originating from the ocean, undergoes isotopic fractionation due to evaporation and condensation during a long transport process, leading to 2 H and 18 O depletion.In winter, the vapor from nearby water bodies, serving as a primary source of humidity, moisture content, and precipitation for Shanghai, typically shows relatively higher isotopic ratios [31,32].
These results suggest that seasonal changes can affect the δ 13 C and δ 15 N values of short-growth cycle BC under each cultivation method due to variable temperatures, light intensity, and fertilizer types.Additionally, seasonal effects can impact the δ 2 H and δ 18 O values through changes in precipitation and/or irrigation water sources.Therefore, when identifying cultivation methods for short-growth cycle vegetables, it is essential to investigate these seasonal effects to ensure that the developed model exhibits high applicability and accuracy across different seasons.
Foods 2024, 13, x FOR PEER REVIEW 7 of 16 negative in winter [12,29].However, it was consistent with the seasonal δ 2 H and δ 18 O variation in local atmospheric precipitation and the irrigation water source from the Yangtze River [30], which are significantly influenced by the distinct monsoon system and local topography around Shanghai [31][32][33].In summer, precipitation in Shanghai, originating from the ocean, undergoes isotopic fractionation due to evaporation and condensation during a long transport process, leading to 2 H and 18 O depletion.In winter, the vapor from nearby water bodies, serving as a primary source of humidity, moisture content, and precipitation for Shanghai, typically shows relatively higher isotopic ratios [31,32].These results suggest that seasonal changes can affect the δ 13 C and δ 15 N values of short-growth cycle BC under each cultivation method due to variable temperatures, light intensity, and fertilizer types.Additionally, seasonal effects can impact the δ 2 H and δ 18 O values through changes in precipitation and/or irrigation water sources.Therefore, when identifying cultivation methods for short-growth cycle vegetables, it is essential to investigate these seasonal effects to ensure that the developed model exhibits high applicability and accuracy across different seasons.

PLS-DA Isotope Models to Identify BC Cultivation Methods
Among the overall stable isotopes of different BC cultivation methods, only the δ 15 N values were significantly different among ORG, GRE, and CON BCs (Figure 1).The three CON BC outliers with high δ 15 N values (Figure 1, blue circle) came from farms transitioning from conventional to organic cultivation and could not represent the typical values of CON BC.After excluding these outliers, 80.95% of ORG BC had δ 15 N values above 6‰, while 43.18% of GRE BC and 81.54% of CON BC were below 6‰.The ORG and CON BCs could be well distinguished, but GRE BC overlapped both ORG and CON BCs due to the use of both fertilization methods.There were no significant differences in the overall mean δ 13 C, δ 2 H, and δ 18 O values of BC among different cultivation methods, as well as their seasonal mean δ 13 C and δ 2 H values.However, stable isotopes of individual BCs still exhibited variations within cultivation method classes due to seasonal variations in temperature, light intensity, and precipitation.A supervised PLS-DA method (IRMS-PLS-DA) was used to investigate these differences and attempt to improve the identification accuracy for ORG, GRE, and CON BCs (Tables 3 and S1), especially for GRE BC.Moreover, the three outliers were included in the modeling data in order to ensure the universality of PLS-DA identification models.

PLS-DA Isotope Models to Identify BC Cultivation Methods
Among the overall stable isotopes of different BC cultivation methods, only the δ 15 N values were significantly different among ORG, GRE, and CON BCs (Figure 1).The three CON BC outliers with high δ 15 N values (Figure 1, blue circle) came from farms transitioning from conventional to organic cultivation and could not represent the typical values of CON BC.After excluding these outliers, 80.95% of ORG BC had δ 15 N values above 6‰, while 43.18% of GRE BC and 81.54% of CON BC were below 6‰.The ORG and CON BCs could be well distinguished, but GRE BC overlapped both ORG and CON BCs due to the use of both fertilization methods.There were no significant differences in the overall mean δ 13 C, δ 2 H, and δ 18 O values of BC among different cultivation methods, as well as their seasonal mean δ 13 C and δ 2 H values.However, stable isotopes of individual BCs still exhibited variations within cultivation method classes due to seasonal variations in temperature, light intensity, and precipitation.A supervised PLS-DA method (IRMS-PLS-DA) was used to investigate these differences and attempt to improve the identification accuracy for ORG, GRE, and CON BCs (Table 3 and Table S1), especially for GRE BC.Moreover, the three outliers were included in the modeling data in order to ensure the universality of PLS-DA identification models.
The first two principal component score plots of PLS-DA (Figure 3) revealed that most ORG BCs could be effectively distinguished from CON BCs, but GRE BCs exhibited significant overlap with both ORG and CON BCs.Therefore, the discriminant model of the three cultivation methods was established in pairs.The calibration model accuracy for ORG and CON BCs was 77.55%, and the validation accuracy was 75.76%, with 24 ORG BCs being misclassified as CON BCs and six CON BCs as ORG BCs, possibly due to variations in fertilizer types and the comprehensive effects of seasonal variations and fertilizer types.For instance, some ORG BCs (n = 11) with δ 15 N values below 5‰ might have utilized plant-based (legume) organic fertilizers or the inadequate fermentation of manures [25], while some transitioning CON BCs only utilized organic manures.Moreover, some misclassified ORG BCs with δ 15 N values higher than 5‰ might have been influenced by seasonal variations in temperature, light intensity, and precipitation [9,26,27], and some misclassified CON BCs might be attributed to increased environmental protection Foods 2024, 13, 1165 9 of 15 by farmers, resulting in higher organic fertilizer use [34].The highest misclassification rate occurred in summer (n = 10), probably due to higher levels of photosynthesis, water evaporation, and transpiration decreasing the differences between ORG and CON BCs during the summer [9,27] (Figure 2 and Table 3).The variable importance in projection (VIP) order used in the PLS-DA model was δ 15 N > δ 13 C > δ 18 O > δ 2 H, consistent with the stable isotope variations observed between ORG and CON BCs (Figure 1). between ORG and GRE BCs (Figure 1).Specifically, ORG BCs showed the highest overall mean δ 15 N values (10.50 ± 6.20‰) and the lowest overall mean  The calibration model of the GRE and CON BCs PLS-DA achieved an accuracy of 73.81%, with a validation accuracy of 53.57%.The number of misclassifications for GRE BCs (n = 30) was six times higher than that for CON BCs (n = 5), and misclassifications remained most pronounced during the summer, primarily due to the distinctive fertilization strategies employed by GRE BCs coupled with seasonal effects on the physiological and biochemical reactions of BCs [24,26].The VIP order of the PLS-DA model was δ 15 N > δ 2 H > 1 > δ 18 O > δ 13 C, suggesting that fertilizer type and the irrigation water source influenced by seasonal precipitation were important modeling variables [8], and aligning well with the overall stable isotope variations between GRE and CON BCs (Figure 1).Specifically, GRE BCs exhibited significantly higher δ 15 N values and the highest δ 2 H values compared to CON BCs.
The results indicated that, despite incorporating the differences in four stable iso- The accuracy of the PLS-DA for the calibration model of ORG and GRE BCs was 71.25%, and for the validation model it was 51.85%, with 30 GRE BCs being misclassified as ORG BCs and six ORG BCs as GRE BCs, possibly attributed to smaller differences in δ 15 N values between ORG and GRE BCs (Figure 1), as ORG BCs used only organic fertilizers while GRE BCs could use both organic and chemical fertilizers (GB/T 33761-2017 and NY/T 391-2021).Summer BCs (n = 12) still demonstrated the highest misclassification rate, which may also be due to a combination of different fertilizer effects and seasonal variations [24,26].The VIP order of the ORG and GRE BCs PLS-DA model was δ 15 N > δ 2 H > 1 > δ 13 C > δ 18 O, highlighting the significance of δ 15 N and δ 2 H values in distinguishing between GRE and ORG BCs.This order aligned with the overall stable isotope variations observed between ORG and GRE BCs (Figure 1).Specifically, ORG BCs showed the highest overall mean δ 15 N values (10.50 ± 6.20‰) and the lowest overall mean δ 2 H values (−80.46 ± 9.68‰), while GRE BCs demonstrated the highest overall mean δ 2 H values (−78.61 ± 9.60‰).
The calibration model of the GRE and CON BCs PLS-DA achieved an accuracy of 73.81%, with a validation accuracy of 53.57%.The number of misclassifications for GRE BCs (n = 30) was six times higher than that for CON BCs (n = 5), and misclassifications remained most pronounced during the summer, primarily due to the distinctive fertilization strategies employed by GRE BCs coupled with seasonal effects on the physiological and biochemical reactions of BCs [24,26].The VIP order of the PLS-DA model was δ 15 N > δ 2 H > 1 > δ 18 O > δ 13 C, suggesting that fertilizer type and the irrigation water source influenced by seasonal precipitation were important modeling variables [8], and aligning well with the overall stable isotope variations between GRE and CON BCs (Figure 1).Specifically, GRE BCs exhibited significantly higher δ 15 N values and the highest δ 2 H values compared to CON BCs.
The results indicated that, despite incorporating the differences in four stable isotopes of individual BCs, the BC cultivation methods' identification rates using PLS-DA models were not superior to those achieved by using only δ 15 N values, due to the influence of variations in fertilizer types and seasonal effects.Furthermore, the accuracy of identifying GRE BCs still needed to be improved.Overall, the study showed that BCs exhibited the highest misclassification rate in summer, resulting from the combination of high temperatures, intense sunlight, and frequent precipitation, leading to vigorous growth and the blurring of stable isotope differences among ORG, GRE, and CON BCs.We confirmed that seasonal factors had important impacts on the stable isotopes of BCs grown using different cultivation methods, consequently impacting the ability to clearly distinguish different BCs' cultivation methods by isotopes alone.

NIR Spectra to Identify BC Cultivation Methods
The rapid identification of agricultural product cultivation methods using NIRS mainly depends on distinctive nutritional composition signals in the spectra from different cultivation methods.Hydrogen-containing BC nutritional components primarily consist of dietary fiber, small amounts of sugars, proteins, and fats.Raw spectra of the ORG, GRE, and CON BCs had a similar spectral shape with significant wavenumber peaks observed at 4010 cm −1 , 4250 cm −1 , 4330 cm −1 , 4670 cm −1 , 5050 cm −1 , 5170 cm −1 , 5780 cm −1 , 6350 cm −1 , 6780 cm −1 , and 8370 cm −1 , corresponding to characteristic groups of the main nutrients in BC (Figure 4a) [23,35].The peak at 4010 cm −1 corresponds to the combined frequency of C-H stretching and C-C stretching, indicating the presence of cellulose; both the peaks at 4250 cm −1 and 4330 cm −1 represent a second-order frequency doubling of C-H bending vibration in C-H groups, potentially indicating the presence of polysaccharides and lipids, respectively; the peak at 4670 cm −1 denotes the combination frequency of C-H stretching and C-H deformation, suggesting the presence of fats; the peak of 5050 cm −1 may indicate the combination frequency of N-H antisymmetric stretching and N-H in-plane bending in CONH 2 groups of amide II, hinting at the presence of proteins; the peak at 5170 cm −1 represents the combination frequency of O-H stretching and HOH deformation in OH and HOH groups, indicating the existence of polysaccharides; the peak at 5780 cm −1 indicates the first-order frequency doubling of C-H stretching vibration in methylene groups of hydrocarbon structures; the peaks at 6350 cm −1 and 6780 cm −1 are the first-order frequency doubling of N-H stretching vibration in amide groups, indicating the presence of proteins; and the peak at 8370 cm −1 might be the second-order frequency doubling of C-H stretching vibration in methyl groups [23,35].Deviations in peak positions occurred for individual BCs due to the influences of cultivation methods and seasons.5780 cm −1 indicates the first-order frequency doubling of C-H stretching vibration in methylene groups of hydrocarbon structures; the peaks at 6350 cm −1 and 6780 cm −1 are the firstorder frequency doubling of N-H stretching vibration in amide groups, indicating the presence of proteins; and the peak at 8370 cm −1 might be the second-order frequency doubling of C-H stretching vibration in methyl groups [23,35].Deviations in peak positions occurred for individual BCs due to the influences of cultivation methods and seasons.Overall, the BC raw spectra overlapped significantly (Figure 4a), as did the spectra of the three cultivation methods under each season (Figure S1).However, BC spectra tended to cluster under different seasons (Figure 4b), suggesting that seasons played an important role in the accumulation of nutrients in BC.The raw spectra of BCs displayed baseline drifts, band overlapping, and weak characteristic peaks (Figure 4a,b), making it challenging to directly distinguish the three cultivation methods based on their spectra.Therefore, NSD was used to reduce interferences and enhance feature signals in the NIR spectra (Figure 4c,d) [20], and the PLS-DA was utilized to build identification models for BC cultivation methods (Table 3).
Based on the BC raw spectra, the calibration model for ORG and CON BCs achieved an accuracy of 87.76%, and the validation accuracy was 78.79%.It was hard to visually distinguish between different BC cultivation methods from the first two principal component score plots of PLS-DA (Figure S2a) due to the NIR spectral data having 1557 variables and a high number of optLVs in modeling (generally more than 10, Table S2).More ORG BCs (n = 13) were wrongly classified compared to CON BCs (n = 7).In addition, summer (n = 7) and winter BCs (n = 5) accounted for a higher proportion of misclassified samples, suggesting possible similarities in nutritional components between ORG and CON BCs during these seasons.The NSD (5,5,2) pretreatment improved the model accuracies, reaching 91.84% for the calibration model and 81.82% for the validation model (Tables 3 and S2).Misclassified samples still occurred more frequently in winter (n = 7) and summer (n = 4), possibly due to slower or faster growth in winter or summer, respectively, leading to similar accumulation rates of nutritional components and lowering cultivation methods differences.The PLS-DA model for ORG and GRE BCs achieved 100% calibration model accuracy and 62.96% validation model accuracy, and 60% (6/10) misclassification samples occurred in summer (n = 4) and winter (n = 2).NSD(9,9,2) preprocessing improved the validation accuracy to 70.37%.The highest number of misclassifications occurred in spring (n = 4), possibly indicating that NSD optimized BCs spectra during summer while reducing spectral information differences in spring.The PLS-DA model for GRE and CON BCs achieved a calibration accuracy of 96.43% and a validation accuracy of 71.43%.Overall, 81.82% (9/11) of misclassified BCs were from summer (n = 7) and winter (n = 2).However, the NSD preprocessing did not improve the predictive abilities of the model, with a 67.86% accuracy for the validation model, suggesting that useful spectral information might be removed when reducing disturbing signals.
Therefore, the optimal PLS-DA model for ORG and CON BCs showed a favorable predictive performance with an accuracy of 81.82% due to differences in their nutritional composition, while the predictive accuracies for ORG and GRE BCs (70.37%) and GRE and CON BCs (75.00%) did not achieve such good results, mainly due to the special cultivation requirements of GRE BCs.In general, higher BC misclassification rates occurred in summer and winter, possibly attributed to the dynamic physiological and biochemical responses of BC during these seasons.These responses were influenced by temperature and light intensity, leading to similar nutrient compositions and, thus, similar NIR spectral signals.
The results confirm the importance of investigating seasonal effects on NIR spectra to build a higher-accuracy and more widely adaptable model for identifying different BC cultivation methods.

Combined IRMS and NIRS to Identify BC Cultivation Methods
Individual IRMS or NIRS PLS-DA models could not effectively identify BCs under different cultivation methods from the same geological origin due to the influence of fertilizer types and seasonal variations.The stable isotopic and NIR spectral differences in ORG, GRE, and CON BCs were comprehensively evaluated to build the identification models of IRMS combined with NIRS (IRMS-NIRS-PLS-DA) (Tables 3, S3 and S4).
The IRMS-NIRS-PLS-DA model for ORG and CON BCs showed a higher predictive accuracy of 87.88% compared to the optimal IRMS-PLS-DA (75.76%) or NIRS-PLS-DA model (81.82%).Furthermore, the NSD(5,5,2) pretreatment optimized the calibration model, further improving the accuracy to 89.80%, although the accuracy of the new validation model remained unchanged, indicating the combined PLS-DA model exhibited a more robust and predictive ability.However, the first two principal component score plots showed an overlap between ORG and CON BCs (Figure S2b).The number of misclassifications in the optimal NSD-PLS-DA model for ORG BCs (n = 9) was higher than for CON BCs (n = 5), possibly due to fertilizer types and seasonal effects [8,26].The misclassified BCs mainly came from winter (n = 6) and summer (n = 4), possibly attributed to relative weak (winter) or strong (summer) photosynthesis, water evaporation, and transpiration, decreasing the differences in stable isotopes and nutritional compositions between ORG and CON BCs during these two seasons (Figure 2 and Table 3).The δ 15 N values still remained the most important variable for identifying between ORG and CON BCs according to the VIP of the NSD-PLS-DA model.
The IRMS-NIRS-PLS-DA model significantly improved the predictive performance for distinguishing between ORG and GRE BCs with an 81.48% accuracy, surpassing the ± 1.36‰) and δ 2 H (−80.46 ± 9.68‰) values.Conversely, CON BC exhibited the highest overall mean δ 13 C (−28.78 ± 1.48‰) and the lowest overall mean δ 15 N (3.58 ± 5.33‰) and δ 18 O (20.73 ± 2.29‰) values, and GRE BC had the highest overall mean δ 2 H (−78.61 ± 9.60‰) values.Only the δ 15 N values showed a significant difference (p < 0.05) resulting from all the study samples collecting from Shanghai farms within a radius of around 40 km and experiencing similar climatic influences throughout the year.ORG, GRE, and CON BCs have different fertilization requirements.Specifically, ORG BC exclusively uses organic fertilizers (GB/T 19630-2019) [18], GRE BC is permitted to use chemical fertilizers in appropriate amounts (GB/T 33761-2017 and NY/T 391-2021) [3,4], and CON BC does not have significant restrictions on the use of chemical fertilizers.

Figure 1 .
Figure 1.Stable isotope values of BC samples from different cultivation methods where "×" sents mean values; "-" represents median values; and "°" represents suspected outliers.D lowercase letters among different cultivation methods indicate a significant difference at the level.The δ 13 C, δ 15 N, δ 2 H, and δ 18 O values among ORG, GRE, and CON BC varied different seasons (Table2).There were no significant differences in seasonal mea and δ 2 H values among the three cultivation methods of BC, possibly indicating th photosynthetic efficiency and water use efficiency of BC under different cultivation ods were similar in each season, given their similar growing locality in Shanghai.Se mean δ15 N values of BC also varied similarly among the three cultivation methods, f ing the sequence ORG > GRE > CON.Significant differences in δ 15 N values were ob between ORG and CON BCs during winter, summer, and autumn-re (p < 0.05).CO exhibited significant differences compared to ORG and GRE BCs (p < 0.05) in sprin significant differences were observed in autumn, possibly due to fertilizer type var and the two outlier values (20.20‰ and 22.53‰) of CON BCs in autumn.ORG BC the highest seasonal mean δ 18 O value (23.37 ± 2.61‰) in spring and were significant ferent from CON BCs (21.65 ± 1.50‰) (p < 0.05).This possibly arose from the enhanc of soil permeability and water retention caused by organic fertilizer application a

Figure 1 .
Figure 1.Stable isotope values of BC samples from different cultivation methods where "×" represents mean values; "-" represents median values; and " • " represents suspected outliers.Different lowercase letters among different cultivation methods indicate a significant difference at the p < 0.05 level.

Figure 2 .
Figure 2. Stable isotopes of different BC cultivation methods based on seasonal time series where "×" represents mean values; "-" represents median values; "°" represents suspected outliers.Different lowercase letters among different seasons indicate a significant difference at the p < 0.05 level.

Figure 2 .
Figure 2. Stable isotopes of different BC cultivation methods based on seasonal time series where "×" represents mean values; "-" represents median values; " • " represents suspected outliers.Different lowercase letters among different seasons indicate a significant difference at the p < 0.05 level.

Figure 3 .
Figure 3.The first two principal component score plots of PLS-DA for different BC cultivation methods using IRMS.

Figure 3 .
Figure 3.The first two principal component score plots of PLS-DA for different BC cultivation methods using IRMS.

Table 1 .
Number of BC samples collected for each season and cultivation method.

Table 2 .
Stable isotope values of BC under different cultivation methods and seasons.
Different letters within a row indicate a significant difference for each cultivation method (p < 0.05).ORG: organic; GRE: green; CON: conventional.

Table 3 .
Model methods and accuracies of PLS-DA models for BC cultivation methods using IRMS and/or NIRS.