The study area (Figures 1A,B) was a large area of Masson pine artificial pure forest located at an altitude of 820–1,000 m above sea level (a.s.l) within Taijiang County, Guizhou Province, China (26°50′51.85″N, 108°17′32.67″E; Figures 1A,B). The area has a humid subtropical monsoon climate. The annual extreme maximum and minimum temperatures of the Masson pine test area are 37.10°C and −10.30°C, respectively, the annual average temperature is 15.70°C, the frost-free period is 286 days, and the annual precipitation is 1,050–1,400 mm.

In the spring of 2014, an area within the barren mountains at the experimental site was burned, and a hole-shaped area of land was used for this study. Thirty-five half-sib Masson pine families were planted from left to right using green belt afforestation (Figure 1B). Fourteen of the families were of Guangxi province and named GX01 to GX14, and 20 families were of Guizhou province and were named GZ01 to GZ21. One-year-old seedlings from these 35 families were planted for afforestation. The planting holes measured 40 cm × 30 cm × 30 cm, and the plant row spacing was 2.0 m × 2.5 m. The same management method, weeding by hand and no fertilization, was used for all families.

In the spring of 2021, the needles of the 35 families were collected to analyze the different metabolites. The needles from four directions of the same tree canopy were pooled into one sample, three trees from each family were collected as samples, and the fresh samples were stored in deep freezer at −80°C until further metabolome analysis. The tree height (H, in m) was measured with a measuring rod. The DBH (in cm) was measured with a diameter tape. The tree height and DBH of the selected families were measured in 2021. Tree volume (V, in m³) was calculated based on the formula for the local living stand volume of the Masson pine. The changes in DBH and tree volume were relatively consistent among families, and the tree height of the families of Guizhou province was greater than those of Guangxi province on average. GX10, GZ16, and GZ13 had the largest average tree height; GX10, GZ16, and GX05 had the largest DBH; and GX10, GZ16, and GX05 had the largest average standing volume (Figure 1C).

The test sample, consisting of 50 mg of needles cut into small pieces, was placed in a 10-mL plastic centrifuge tube containing 30 μl of mixed internal standard (adipic acid 10.04 mg/ml, phenyl glucoside 8.04 mg/ml, positive valine 4.9 mg/ml, and methanol and water, 1:1) and 3 ml of extraction solvent mix (methanol:chloroform:water = 2.5:1:1). After homogenizing for 1 min, the sample was extracted, sonicated at 4°C for 20 min, and centrifuged at 3,000–5,000 rpm for 5 min. Afterwards, a 150-μl aliquot of the supernatant extract was added to a 1.5-ml centrifuge tube, dried with nitrogen, and set aside.

Subsequently, 40 μl of 25 mg/ml methoxyamine hydrochloride was added to the supernatant extract (the methanol-water phase), sealed with parafilm, vortexed-mixed for 5 min, and incubated in a constant temperature water bath at 40°C for 120 min (briefly vortexing once at 60 min for approximately 30 s). After cooling to 25°C, briefly centrifuging, and adding 200 μl of bis(trimethylsilyl)trifluoroacetamide (trimethylchlorosilane 1%), the tube was sealed with parafilm, vortex-mixed for 30 s, and reacted in a thermostatic water bath at 81°C for 120 min (briefly vortexing at 60 min for 10 s). Subsequently, after cooling the reaction mixture to room temperature, vortexing for 2 min, and centrifuging at 12,000 rpm for 5 min, the supernatant was then injected into a linear tube. The methanol–water phase was derivatized and analyzed by GC–MS.

GC-MS analysis was performed on an Agilent 7890A GC system coupled with an Agilent 5975C mass spectrometer detector, equipped with an HP-5 MS (60 m × 250 μm × 0.25 μm) capillary column (Agilent Technologies Inc., Santa Clara, CA, United States), using the following conditions: inlet temperature of 280°C, injection volume of 2 μl, split ratio of 20:1, and column flow rate of 1.0 ml/min. The heating program was run for a total of 72 min, and the following temperatures and times were used: 60°C for 2 min, ramping at 5°C/min to 230°C for 5 min, and 8°C/min to 290°C for 23.5 min.

The MS conditions were as follows: ion source temperature of 230°C, quadrupole temperature of 150°C, ionization energy of 70 eV, transmission line temperature of 280°C, full scan mass range of 45–600 aum, solvent delay of 11.90 min, and the acquisition mode was full-scan acquisition mode.

Metabolites were identified using the GC–MS spectral libraries (NIST08 and Willy08 libraries). The chemicals used were qualitative and quantitative with ions, the information obtained was organized, and the specific metabolites of the different Masson pine families were characterized. The metabolites were classified using standard methods, and referencing the associated literature (Zhang et al., 2020). The statistical analyses were performed using the R (version 4.0.0) software (Kanehisa and Goto, 2000; Chong and Xia, 2018).

The STAMP (version 2.1.3) software (Parks et al., 2014) with Welch’s t-test was used to determine the difference in metabolite profiles between the Masson pine families from the two geographical provenances at α = 0.05, and log₁₀ (x + 1) was used to convert the analysis data before analysis. This analysis was performed using the R (version 4.0.0) software (R Core Team, 2019), in particular, the “tidyverse” and “ggplot2” (Wickham, 2016) packages were used for the visualization of the results.

Principal component analysis (PCA): Unsupervised PCA was performed using the statistics function prcomp within the R (version 4.0.0.) software. The data were scaled to unit variance before unsupervised PCA (Ren et al., 2022).

Identification of differential metabolites: significantly differentially regulated metabolites between the Masson pine families from two provenances were determined by variable importance in project (VIP) ≥ 1 and absolute log₂FC (fold change) ≥ 1. The VIP values were extracted from the orthogonal partial least squares discriminant analysis (OPLS-DA) results, and were determined using the R package MetaboAnalystR: OPLSR. Anal in the R software. To avoid overfitting, a permutation test (200 permutations) was performed. The metabolic group data were analyzed according to the OPLS-DA model, the score chart of each group was drawn, and the differences between each group were further determined (Thévenot et al., 2015). The prediction parameters of the evaluation model included R2X, R2Y, and Q2, where R2X and R2Y represent the interpretation rate of the model to the X and Y matrices, respectively, and Q2 represents the predictive ability of the model. The closer these three indicators are to 1, the more stable and reliable the model is. When Q2 > 0.5, it can be considered as an effective model, and when Q2 > 0.9, it is an excellent model.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation and pathway enrichment analysis: the identified differential metabolites were annotated using the KEGG Compound database,¹ and the annotated metabolites were then mapped to the KEGG Pathway database.² Pathways with significantly differentially regulated metabolites were mapped and then input into the Metabolite Set Enrichment Analysis (MSEA) and enrichment was assessed using the MSEA Metabolic Pathway library.

In total, 116 metabolites were identified and grouped into 10 groups, and the main chemical groups were saccharides and alcohols, amino acids and derivatives, organic acids, amines, lipids, phosphoric acids, phenolic acids, nucleotides and derivatives, sterols, etc. There were 43 saccharides and alcohols, 23 amino acids and derivatives, and 22 organic acids (Figure 2A). These metabolites are predominantly the products or intermediates of Masson pine photosynthesis, which is the main physiological function of needles. The average content of organic acid was the highest (102,649 μg/g), followed by that of saccharides and alcohols (92,242 μg/g) and phosphoric acids (2,945 μg/g), while the content of other groups was lower than 2,000 μg/g (Figure 2B).

The contents of organic acids, amines, and others groups were significantly different between the Masson pine from Guangxi and Guizhou provinces (p < 0.05), and other groups showed no significant difference between Masson pine from the two provinces (p > 0.05). The content of organic acids of Masson pine from Guangxi province is significantly higher than that of Masson pine from Guizhou province (p < 0.05), but the content of other groups is lower than that of Masson pine from Guizhou province, among which the content of amines and others are significantly lower than that of Masson pine from Guizhou province (Figure 3).

PCA can provide a preliminary understanding of the overall metabolic differences between Masson pine families from various provenances and the variability between families within each province. To better determine the differences between the two provinces, we used both PCA and OPLS-DA methods. The results showed that there was an obvious separation between the Masson pine families from Guangxi and Guizhou provinces (Figure 4A), and three principal components explained 34.79% of the sample information in total (Figure 4B).

OPLS-DA can discriminate the differences between the two groups and improve the effectiveness and analysis ability of the model. Therefore, we used R2X, R2Y, and Q2 of the OPLS-DA model to verify its reliability. For the OPLS-DA model obtained from the data of the Masson pine families from the two provinces, these three parameters were 0.267, 0.878 and 0.75, respectively (Figure 5A). In addition, we used random permutation and combination experiments to evaluate the accuracy of the OPLS-DA model. The R2 and Q2 values on the left side were less than the initial values, indicating that the model was reliable (Thévenot et al., 2015). The OPLS-DA scores showed a clear separation between the Masson pine families from the two provinces (Figure 5B). The metabolites with VIP ≥ 1 were determined by OPLS-DA S-plot (Figure 5C). The top 10 upregulated and downregulated metabolites of GZ vs. GX are shown in Figure 5D.

With the threshold value of VIP > 1.0 and absolute log₂FC ≥ 1, we identified six differential metabolites between Masson pine families from Guizhou and Guangxi provinces, including three upregulated and three downregulated metabolites (Figure 6A). According to the VIP value, the differential metabolites were caffeic acid, L-ascorbic acid, gentiobiose, xylitol, d-pinitol, and β-sitosterol (Figure 6B). These differential metabolites play an important role in distinguishing the provenance of Masson pine families.

To examine the relationship between different metabolites and metabolic pathways among Masson pine families from different provenances, we searched the Human Metabolome Database (HMDB) or KEGG database for four significantly differentially regulated metabolites, namely caffeic acid, L-ascorbic acid, xylitol, and β-sitosterol. We used the KEGG database to annotate and display the differentially regulated metabolites (Kanehisa and Goto, 2000), and metabolism and environmental information processing were the two main metabolic pathways enriched in the differentially regulated metabolites. Among them, metabolism includes metabolic pathways, biosynthesis of secondary metabolites, steroid biosynthesis, phenylpropanoid biosynthesis, pentose and glucuronate interconversions, glutathione metabolism, biosynthesis of cofactors, ascorbate and aldarate metabolism, while environmental information processing includes ABC transporters (Figure 7A). These metabolites enriched in these pathways may be the most important metabolites determining potential provenance differences. In addition, the KEGG pathways enrichment analysis of the differentially regulated metabolites showed that the pathways most significantly enriched in differentially regulated metabolites between the Masson pine families from the two provenances were steroid biosynthesis, phenylpropanoid biosynthesis, glutathione metabolism, pentose and glucuronate interconversions (Figure 7B).

The study of plant metabolic components in different Masson pine tree families is of great significance for the identification of different Masson pine families from various provinces. Genetic information, tree height, volume, branch growth, nitrogen and phosphorus nutrient efficiency differ considerably among the various Masson pine families (Hu et al., 2016; Ye et al., 2020; Ni et al., 2021; Shi et al., 2021). However, little attention has been paid to the metabolic differences of different Masson pine families from various provinces.

All Masson pine families have the same characteristics metabolite profile, and the main chemical groups were saccharides and alcohols, amino acids and derivatives, organic acids, amines, lipids, phosphoric acids, phenolic acids, nucleotides and derivatives, sterols, and others. Remarkably, the results of this study showed that the Masson pine families from two provinces showed different metabolite characteristics. The content of organic acids of Masson pine families from the Guangxi province was higher than that of Masson pine families from the Guizhou province, but the other groups were lower than those of Masson pine families from the Guizhou province. The Masson pine families from the Guizhou and Guangxi provinces could be clearly separated by PCA based on differentially regulated metabolites. These results suggested that some key differentially regulated metabolites were formed during the domestication of Masson pine trees from different provinces, indicating that they played an important role in such process. Therefore, our study once again verified that the profile of metabolites in all pines was determined by the interaction of genetic background and the abiotic and biotic factors, and the content of metabolites in the Masson pine families was affected by the geographical environment.

Plants synthesize hundreds of thousands of metabolites, and different species and genera produce different types of primary and secondary metabolites (Allwood et al., 2008). Both primary and secondary metabolites play important roles in plant stress resistance (Zaynab et al., 2018; Zhang et al., 2020). In this study, only the methanol–water phase was derivatized and analyzed by GC-MS, as a result mainly primary metabolites were identified in this study. A total of 116 metabolites were identified, which were mainly primary metabolites. Primary metabolites are present in nearly every part of the plant species and are synthesized through the same or nearly the same biochemical pathways (Canarini et al., 2019), and these metabolites are critical for plant growth, development, and reproduction (Abdelrahman et al., 2019). These metabolites are produced during the growth phase and function as signaling molecules that trigger defense responses through signal transduction and pathogen-recognition pathways (Kachroo and Robin, 2013), and their induction patterns are similar to those of secondary metabolites (Devoto et al., 2005). For example, proline plays an important role in plant defense mechanisms and stress responses, and biotic and abiotic stresses usually lead to proline accumulation in plants, which serves as a potential indicator of plant stress resistance (Maroli et al., 2016). The needles of Masson pine trees are useful for metabolite research as they are a storage organ and contain metabolites with functional traits. Masson pine terpenoids contribute to drought resistance (Quan and Ding, 2017), and they defend the plant against herbivores (Chen et al., 2019). We identified many kinds of amino acids and organic acids in Masson pine families, and these components are associated with plant growth traits and stress resistance.

This study showed that the content of many metabolites was significantly different between Masson pine families from two provinces, indicating that they participated in the process of domestication and breeding of Masson pine families and may be involved in some metabolic regulation. We found that six deferentially regulated metabolites contributed the most to the metabolic differences of Masson pine families of GZ vs. GX provinces, and these metabolites are the key metabolites for the ecological adaptation of plants (Fernández de Simón et al., 2017; Mitiouchkina et al., 2020; Sánchez-Hidalgo et al., 2021; He et al., 2022). In the genus Pinus, there is great potential for incorporating phytochemical metabolism into conifer breeding programs (Nantongo et al., 2021), and several studies have also directly or indirectly linked primary metabolites to biotic stress responses (Nantongo et al., 2022). Chemical analysis of Pinus has shown that volatile substances in needles differ considerably between species (Ioannou et al., 2014). The results of this study indicate that the d-pinitol in the Masson pine families from Guizhou province is associated with a high accumulation of metabolites during the breeding process, and this provides a direction for the future development of high value-added natural Masson pine products. The inositol ether d-pinitol is a naturally occurring compound found in Pinaceae that has positive effects on human hypoglycemia and cardiovascular protection (Moreira et al., 2018), stimulates insect spawning (Honda et al., 2012), and is associated with hypoglycemic, cardiovascular, antiviral, and larvicidal activities. Xylitol is a naturally occurring alcohol found in most plant materials, including many fruits and vegetables (Cheng et al., 2009), which may be extracted from the hardwood species containing xylan (Ahuja et al., 2020). β-Sitosterol is a sterol found in almost all plants, which as the main component of many plants and vegetables has a variety of biological effects (Liu et al., 2019). β-Sitosterol is also found in some medicinal plants and is widely used in the pharmaceutical industry (Babu and Jayaraman, 2020).

Caffeic acid is found in all plants as it is a key intermediate in the biosynthesis of lignin (Reuter et al., 2020). Previous studies have shown that an increase in caffeic acid content in trees will lead to an increase in lignin content (Bubna et al., 2011; Zhao et al., 2022). Caffeic acid is a common phenolic compound widely distributed in plants, which has broad antibacterial and antiviral activities, as well as a variety of pharmacological properties, such as anti-inflammatory, anti-anxiety, and anti-depressive activities (Deguchi and Ito, 2020; Dukie et al., 2020; Dziedzinski et al., 2020; Grabska-Kobylecka et al., 2020). L-ascorbate acid and β-gentiobiose are commonly used antioxidants, which serve as important physiological indicators of plant resistance to stress and have important significance for plant identification and breeding (Takahashi et al., 2014; Yang et al., 2018; Shah et al., 2022). L-ascorbic acid is an important metabolite that is crucial to plant physiological and biochemical processes, which besides being involved in the synthesis of ethylene, gibberellin, anthocyanin and other molecules involved in various metabolic pathways in plants, it can also improve the anti-stress and antioxidant capacity of crops (Anjum et al., 2014; Noctor et al., 2014; Akram et al., 2017). β-Gentiobiose is a rare disaccharide that accumulates abundantly in plants and is involved in various aspects of plant development, such as fruit ripening and dormancy release of bud and growth dormancy (Dumville and Fry, 2003; Takahashi et al., 2018). In this study, the Masson pine families with high contents of L-ascorbate acid and β-gentiobiose showed stronger stress resistance.