Pinus massoniana Lamb. is an important coniferous tree species in ecological environment construction and sustainable forestry development. The function of gene gradual change and coexpression modules of needle and root parts of P. massoniana under continuous drought stress is unclear. The physiological and transcriptional expression profiles of P. massoniana seedlings from 1a half-sibling progeny during drought stress were measured and analyzed. As a result, under continuous drought conditions, needle peroxidase (POD) activity and proline content continued to increase. The malondialdehyde (MDA) content in roots continuously increased, and the root activity continuously decreased. The needles of P. massoniana seedlings may respond to drought mainly through regulating abscisic acid (ABA) and jasmonic acid (JA) hormone-related pathways. Roots may provide plant growth through fatty acid β-oxidative decomposition, and peroxisomes may contribute to the production of ROS, resulting in the upregulation of the antioxidant defense system. P. massoniana roots and needles may implement the same antioxidant mechanism through the glutathione metabolic pathway
Pinus massoniana Lamb. (Fam.: Pinus; Gen.: Pinus) is widely distributed in China and plays a pivotal role in ecological environment construction and sustainable forestry production. Drought is the most widespread stress factor affecting plant productivity [1]; tree radial growth is mainly limited by dry conditions from May to October, with drought-driven tree mortality [2]. Climate predictions have forecast that hot and dry events will become more frequent in the coming decades [3][4]. Molecular responses to drought stress have been extensively studied in broadleaf species, but studies of needle-leaf species have been limited [5]. Studying the response of needle-leaf plants under drought plays an important role in screening high-quality drought-resistant germplasm resources and revealing the mechanism of drought resistance.
In recent years, the study of drought in P. massoniana has focused on the photosynthetic characteristics, plasma membrane structure, osmotic adjustment, endogenous enzyme activity, and endogenous hormone content of seedlings under drought stress; most studies tend to be basic theoretical research [6]. Drought significantly reduces the photosynthetic rate, increment of stem basal diameter and plant height, and biomass accumulation [7]. P. massoniana enhances its drought resistance by altering the ratio of monoterpenes and sesquiterpenes: the relative contents of monoterpenes decreases with the intensity of drought stress, and the relative content of sesquiterpenes increases [8][9]. The drought resistance of P. massoniana is controlled by genetic factors; the differences in drought resistance between different families may be related to the genetic differences formed by long-term domestication in different geographical environments [10]. The response to stress is through the regulation of numerous complex and interrelated metabolic networks. The transcriptome is essential for interpreting the functional elements of the genome, revealing the molecular components of tissues and cells and providing an understanding of development [11]. The underground part mainly shows the growth strategy of “increasing” the root system of P. massoniana, which grows in an obvious manner. Transcription factors (TFs) play an important role in the adaptation of P. massoniana needles to drought stress [12][13]. However, there is no research on the continuous gradual change of genes and the role and function of needles and roots under continuous stress conditions. In this study, a high-throughput sequencing method is used to identify the gene expression differences between the needles and roots of a half-sibling progeny of P. massoniana under continuous drought. This study provides basic data for screening drought resistance factors and understanding drought resistance mechanisms.
Through physiological measurements of the process of drought stress (Figure 1), soil absolute moisture content (AMC) gradually decreased (Figure 1a); the plant needles’ malondialdehyde (MDA) (ML) activity (Figure 1b) and the content of needle proline (Pro) (Figure 1c) also continued to rise. Needle peroxidase (POD) activity (Figure 1d) reached its peak in 6 d, after which the needle content continued to decline. Root MDA (MR) activity (Figure 1e) and root activity (RV) (Figure 1f) continued to increase over time. Pro and ML showed a significant positive correlation (p = 0.94), while RV and Pro showed a significant negative correlation (p = −0.95) (Figure 1g). Principal component analysis (PCA) analysis (Figure 1h) found that the difference between the 6-d group and other groups was mainly reflected in POD activity.

Figure 1. Physiological changes of P. massoniana under persistent drought stress. (a) Absolute water content in soil (AMC); (b) malondialdehyde (MDA) content in needles; (c) needle proline (Pro) content; (d) needle peroxidase (POD) activity; (e) MDA content in roots; (f) root activity; (g) correlation analysis of various biochemical indicators; (h) PCA analysis. Note: In Figure 1a–f, we used the “0 d” group as a reference to compare the mean and added p-value and significance markers (“*” indicates that there is a difference compared with the “0 d” group; * p < 0.05, ** p < 0.01, *** p < 0.001, **** p ≤ 0.0001 and “ns” means no difference). In Figure 1g, the number represents the magnitude of the correlation, and the color represents the positive and negative correlations.
By analyzing the needle and root coexpression module of P. massoniana during drought, the coexpressed gene set collection was determined, and the analysis revealed several main subexpression profiles with similar expressions, called coexpression modules [14]. We generated a consensus network to identify conserved modules that are needle set-specific and root set-specific in P. massoniana.
A total of 11 modules (24 to 1208 genes) were identified in the needle module, and a total of 14 modules (37 to 932 genes) were identified in the root module (Figure 2a). In the needle module, the genes of the black and magenta submodules showed a downward trend with drought time. GO enrichment showed that flavonoid biosynthesis (GO: 0009813), the flavonoid metabolism process (GO: 0009812), and other biological processes were enriched in the megenta submodule. The genes of the turquoise submodule in needles showed an increasing trend with prolonged drought time (Figure 2b). KEGG pathway enrichment analysis found that glutathione metabolism (ko00480), MAPK signaling pathway-plant (ko04016), flavonoid biosynthesis (ko00941), and other pathways were enriched. In the root module, the turquoise submodule and the needle submodule shared the most genes, up to 762. GO enrichment of these 762 genes showed that phospholipid–hydroperoxideglutathione peroxidase activity, glutathione peroxidase activity (glutathione peroxidase activity), and other molecular functions were enriched.

Figure 2. Root and needle coexpression network analysis of the P. massoniana drought process. (a) Overlap between needle set-specific and root set-specific modules in P. massoniana; (b) different module gene expression patterns in needles and roots. Note: In subfigure (a), each row in the table corresponds to one needle set-specific module, and each column corresponds to one root set-specific module. The number in the labeled heatmap indicates the genes counted in the intersection between two parts of P. massoniana. Coloring of the table indicates significant overlap, evaluated using Fisher’s exact test. In subfigure (b), the left side of the double-axis represents the needles, the right side represents the roots, the text includes the number of modules and module names, and the gene trends in the module are processed by taking the logarithm of the gene expression in the module, fitted in a linear simulation.
Large-scale biogeographical shifts in forest tree distributions are predicted in response to altered precipitation and temperature regimes associated with climate change [15]. Seed sources with higher plasticity perform better in core habitat conditions [16]. Paying attention to the response of forest species to drought will help us use molecular markers and field experiments to screen for drought-resistant core germplasm resources, which are necessary to improve forest productivity under the conditions of climate change. Transcription differences between different tissue parts/organs are essential for a comprehensive understanding of the stress response of the entire plant. By measuring the physiological indexes of needles and roots during drought of the half-sibling progeny of P. massoniana, the gene expression profiles during the drought process were constructed. This study provides further insights into the process of drought response and tolerance in P. massoniana.
In response to increased H2O2 levels due to stress, plants have evolved different enzymatic and nonenzymatic mechanisms, including free radical scavengers such as superoxide dismutase, catalase, and peroxidase, as well as ascorbic acid-glutathione circulating enzymes [17][18]. During drought, total NSC (TNSC, soluble sugars plus starch) values decrease, particularly in leaves of P. sylvestris L. [19]. MDA and Pro, SOD activity, root shoot ratio, and root mass ratio increased significantly during drought treatment in P. massoniana [20]. Under drought conditions, MDA in needles (Figure 1b) and the content of proline continued to rise (Figure 1c), the MDA in roots (Figure 1e) continued to rise, and root vitality (Figure 1f) continued to decline.
As a sun plant, P. massoniana has higher light compensation and saturation points and a stronger ability to adapt to strong light. Under moderate (field-water-holding capacity 35%–45%) and severe (20%–30%) drought stress, net photosynthetic rate and transpiration rate decreased significantly [21]. The reduction of photosynthesis and chlorophyll contents under drought stress is the most direct manifestation of plants under drought stress. KEGG pathway enrichment analysis showed that plant hormone signal transduction (ko04075) was always enriched in the different teams of needles. STEM analysis found that photosynthesis (ko00195), photosynthesis-antenna proteins (ko00196), porphyrin and chlorophyll metabolism (ko00860), and other metabolic pathways were enriched in the needle Profile9 module (Figure 3a.1). The modules’ inner genes showed the expression trend 0, −1, −2, −3, −4. ADH1, as one of the marker genes for drought or ABA response, is induced by high ABA levels in Arabidopsis under drought stress [22][23]. Overexpression of AtADH1 increases the transcription levels of multiple stress-related genes and the accumulation of soluble sugars and callose deposition and enhances the abiotic and biological stress resistance of Arabidopsis [24]. During drought stress induced by PEG in Arabidopsis thaliana, ADH1 and two other ADH genes responded to drought in needles and roots together [25]. Under drought stress, 24 ADH1 were upregulated in the needle Profile41 module, and 12 ADH1 were upregulated in the needle Profile41 module. In cell signaling, protein kinases (PKs) and phosphatases (PPs) are key enzymes in the regulation mechanism of protein reversible phosphorylation [26]. Plant PP2Cs participate in various signal cascades, such as ABA and salicylic acid (SA)–ABA crosstalk [27]. OPR is involved in the biosynthesis of JA. The key process in the synthetic pathway is the conversion of cis-12-oxo-docosadienoic acid (OPDA) to 12-oxo-phytoenoic acid (OPC-8:0), catalyzed by OPR; the reaction takes place on the peroxisomes [28][29][30]. The function of OPR may be conserved in higher plants, and rice OsOPR7 can make up for the phenotype of Arabidopsis opr3 mutants [31]. In the Profile41 module in needles (Figure 3a.2), 8 ADH1, 6 OPRs, and PP2C in the ABA pathway, MYC2 in the JA pathway, and other related genes were enriched, indicating that the needles of P. massoniana may mainly regulate ABA and JA-hormone-synthesis-related pathways in response to drought.

Figure 4. Cluster analysis of gene expression trends during drought stress in P. massoniana. (a.1) and (a.2): Two STEM significance modules (needle Profile9 and needle Profile41) of the needles during drought stress; (b.1) and (b.2): Two short time-series expression miner (STEM) significance modules (root Profile9 and root Profile41) of the roots during drought stress. Note: We used STEM software to perform trend clustering analysis and divided the trends into 5 stages according to the time before and after input and normalized the union of all the combined DEGs from the five stages. The left side of each subfigure is the significant module selected in the STEM analysis; the right side is the box plot of the corresponding gene expression after taking the logarithm.
The major ROS regulatory enzyme systems in plant peroxisomes include catalase, all ascorbic acid-glutathione cycle components, and superoxide dismutase (SOD) [32]. Peroxisomes contribute to several metabolic processes, such as β-oxidation of fatty acids, biosynthesis of ether phospholipids, and metabolism of ROS. β-oxidation of fatty acids and detoxification of ROS are generally accepted as being fundamental functions of peroxisomes [33][34]. β-oxidation consists of four enzymatic steps: ACX oxidation, MFP2 hydration and dehydration, and 3-ketoacyl-CoA thiolase (KAT) thiolysis. β-oxidation degrades fatty acids only in the peroxisomes in plants and their derivatives [35]. Peroxisomal abundance can be used as a cellular sensor for drought and heat stress responses; severe drought stress induces peroxisomal proliferation and affects the size of the peroxisomes [36][37]. GOLS is a key rate-limiting enzyme for raffinose synthesis [38]. Overexpression of AtGOLS2 in rice significantly improved the drought tolerance of transgenic rice; at the same time, the yield of transgenic rice lines under drought conditions was significantly higher than that of the control [39]. In the root Profie41 module (Figure 3b.2), ko00071, ko00592, and ko01212 pathways were enriched, and genes such as MFP2, ADH1, ADH7A1, and ACOX were enriched. In the ko00052 pathway, genes such as GOLS were enriched, and related differential genes presented gene expression trends of (0,1,2,3,4). KEGG pathway enrichment analysis found that phenylpropanoid biosynthesis (ko00940) was enriched in the root Profile9 module (Figure 3b.1); most (25/31) of the genes belong to peroxidase (K00430). β-oxidation is believed to enhance the production of stress-induced ROS, which adversely affects plant survival under adverse environmental conditions [40]. Under drought stress, plants tend to increase ROS levels, which, in turn, leads to an upregulation of the antioxidant defense system [41]. This shows that the declining peroxidase-related genes may act as inhibitor agents, but the related regulatory mechanisms are currently unclear. On the other hand, fatty acids are catabolized by β-oxidation to acetyl-CoA, producing sugar, glyoxylic acid, and sugar for seedling growth and development without photosynthesis xenobiosis [42]. This shows that under continuous drought conditions of P. massoniana, the roots may provide plant growth through fatty acid β-oxidative decomposition, and P. massoniana peroxisomes may contribute to the production of ROS, resulting in the upregulation of the antioxidant defense system.
P. pinaster root, stem, mature, and aerial organs, under nontarget metabolomics analysis, were found to have a large number of flavonoids, detected in aerial organs, mainly in root glutathione pathway induction. It is believed that aerial organs and roots activate different antioxidant mechanisms [43]. Transcriptome clustering analysis of P. halepensis found that ROS were cleared by an ascorbic acid–glutathione cycle, fatty acid and cell wall biosynthesis, stomatal activity, and biosynthesis of flavonoids and terpenoids [5]. The drought-tolerant seed sources of P. halepensis showed increased levels of glutathione, methionine, and cysteine [44]. The P. massoniana glutathione peroxidase gene, PmGPX6, is highly expressed in the roots. Overexpression of PmGPX6 in Arabidopsis and wild-type plants have little difference in phenotype and root length under normal water conditions, but under drought stress, in transgenic plants, the root system is longer [45]. Phospholipid hydroperoxide glutathione peroxidase, peroxidase, and other activities are related to reducing the accumulation of reactive oxygen in stress-induced environments [46]. The coexpression module of phospholipid–hydroperoxideglutathione peroxidase activity, glutathione peroxidase activity, and other molecular functions in the roots and needles of P. massoniana were enriched. It shows that the glutathione metabolic pathway between the roots and needles of P. massoniana may implement the same mechanism of activating antioxidants.