The miRNA-mRNA Regulatory Modules of Pinus massoniana Lamb. in Response to Drought Stress

Masson pine (Pinus massoniana Lamb.) is a major fast-growing woody tree species and pioneer species for afforestation in barren sites in southern China. However, the regulatory mechanism of gene expression in P. massoniana under drought remains unclear. To uncover candidate microRNAs, their expression profiles, and microRNA-mRNA interactions, small RNA-seq was used to investigate the transcriptome from seedling roots under drought and rewatering in P. massoniana. A total of 421 plant microRNAs were identified. Pairwise differential expression analysis between treatment and control groups unveiled 134, 156, and 96 differential expressed microRNAs at three stages. These constitute 248 unique microRNAs, which were subsequently categorized into six clusters based on their expression profiles. Degradome sequencing revealed that these 248 differentially expressed microRNAs targeted 2069 genes. Gene Ontology enrichment analysis suggested that these target genes were related to translational and posttranslational regulation, cell wall modification, and reactive oxygen species scavenging. miRNAs such as miR482, miR398, miR11571, miR396, miR166, miRN88, and miRN74, along with their target genes annotated as F-box/kelch-repeat protein, 60S ribosomal protein, copper-zinc superoxide dismutase, luminal-binding protein, S-adenosylmethionine synthase, and Early Responsive to Dehydration Stress may play critical roles in drought response. This study provides insights into microRNA responsive to drought and rewatering in Masson pine and advances the understanding of drought tolerance mechanisms in Pinus.


Introduction
Drought is one of the most significant natural environmental factors that affect plant growth, yield, and survival [1]. As sessile organisms, plants have evolved mature mechanisms to cope with drought, for example, stomatal regulation, protective solute accumulation, reactive oxygen species (ROS) detoxification, and cell wall stiffening [2]. Woody tree plants are often challenged by drought stress during their long lifespan [3]. Many largescale forest mortality events caused by drought have been documented [4,5]. Moreover, the ongoing global climate change is making drought events more frequent, longer-lasting, and more intense [6][7][8], resulting in severer loss. Therefore, investigation into the underlying mechanism of how woody tree plants respond to drought will help to improve their drought tolerance and maintain growth and productivity [9,10].
Roots play an important role in plant responses to drought. They are responsible for water uptake in the whole plant and are also the first organ to sense soil-borne water were aligned to the P. massoniana transcriptome from our previous study on the same set of samples [46]. The mapping rates ranged from 63.71% to 78.93% across the 21 libraries (Table S1). A total of 421 miRNAs were identified, mainly at 21 nt, followed by 22 nt and 20 nt, comprising 261 (62.00%), 138 (32.8%), and 22 (5.2%) counts, respectively ( Figure 1B). Among them, 290 and 131 were known and novel miRNAs, respectively. Within the 290 known miRNAs, 21,152, and 117 were 20 nt, 21 nt, and 22 nt long, respectively. These 290 known miRNAs belonged to 38 miRNA families, such as miR950 (46 members), miR946 (33 members), and miR482 (30 members, Figure 1C). Among the 131 novel miRNAs, 1, 109, and 21 were 20 nt, 21 nt, and 22 nt long, respectively. These 131 novel miRNAs represented 60 families, for instance, miRN17 (nine members), miRN54 (eight members), and miRN11 (seven members, Table S2). The nucleotide composition of the mature miRNA sequence was evaluated. Notably, these miRNAs exhibited a preference for a 5′-uridine residue ( Figure S1).

Expression Profile of DEMs
Gene expression profiles can help to predict gene function [47]. Therefore, the DEMs were categorized into six clusters using K-means clustering analysis according to expression profiles ( Figure 3A, Table S4). Cluster 1 and cluster 2 contained DEMs with peak expressions in D1 versus C1 and D3 versus C3 (drought conditions, Figure 3A,B, Table S4). Cluster 1 primarily featured miRNAs from miR946 and miR11425 families and cluster 2 exhibited a prevalence of miRNAs from miR3710 and miRN11 families. In Cluster 3, DEMs peaked in expression during D3 versus C3 (drought condition) and they were mainly from miR950 and miRN9 families. Clusters 4 and 5 contained DEMs whose expressions were peak in D2 versus C2 (post-rewatering). In cluster 4, DEMs were predominately from miR166 and miR11487 families ( Figure 3C, Table S4). In cluster 5, the majority of DEMs were from miRN17 and miR11487 families. In cluster 6, DEMs culminated

Expression Profile of DEMs
Gene expression profiles can help to predict gene function [47]. Therefore, the DEMs were categorized into six clusters using K-means clustering analysis according to expression profiles ( Figure 3A, Table S4). Cluster 1 and cluster 2 contained DEMs with peak expressions in D1 versus C1 and D3 versus C3 (drought conditions, Figure 3A,B, Table S4). Cluster 1 primarily featured miRNAs from miR946 and miR11425 families and cluster 2 exhibited a prevalence of miRNAs from miR3710 and miRN11 families. In Cluster 3, DEMs peaked in expression during D3 versus C3 (drought condition) and they were mainly from miR950 and miRN9 families. Clusters 4 and 5 contained DEMs whose expressions were peak in D2 versus C2 (post-rewatering). In cluster 4, DEMs were predominately from miR166 and miR11487 families ( Figure 3C, Table S4). In cluster 5, the majority of DEMs were from miRN17 and miR11487 families. In cluster 6, DEMs culminated expression in D1 versus C1 (drought condition). The majority of DEMs in this cluster belonged to miR482 and miRN54 families. expression in D1 versus C1 (drought condition). The majority of DEMs in this cluster belonged to miR482 and miRN54 families.

Target Gene Prediction via Degradome Sequencing
Globally, 16,705,004 sequences were obtained from degradome sequencing, and 16,697,634 clean reads sequences were obtained after fastp quality control. Q20 and Q30 sequences accounted for 98.4% and 43.8% of the clean read sequences, respectively (Table  S5). Within the set of 421 miRNAs, there were 419 that targeted 3582 mRNA, forming 15,522 miRNA-mRNA modules (Table S6). Of these, 2942 modules were 'category 0', 1626 modules were 'category 1', and 6995 modules were 'category 2'. In the subset of 248 DEMs, 247 targeted 2069 genes, of which 1772 were annotated by SwissProt (Table S7). Out of 1772 annotated genes, 155 isoforms were annotated as transcription factor genes. The top five transcription factors with the highest number of genes were squamosa promoterbinding-like protein (SPL, 28 genes), dehydrin (26 genes), growth-regulating factor (23 genes), homeobox-leucine zipper protein (15 genes), transcription factor GAMYB (15 genes), and NAC domain-containing protein (nine genes). On the other hand, 1617 isoforms were annotated as functional protein genes. The top five categories of functional protein genes, ranked by the number of genes, were disease resistance protein (245 genes), TMV resistance protein (103 genes), probable disease resistance protein (60 genes), 60S ribosomal protein (29 genes), and 40S ribosomal protein (19 genes).
Further, GO term enrichment analysis was performed on DEM target genes in each cluster to infer the mediatory role of DEMs during drought and rehydration. DEMs in cluster 1 were linked to GO terms such as "oxidoreductase activity using superoxide radicals as acceptor" and "superoxide dismutase activity" ( Figure 3C, Table S8). The DEMs in cluster 2 were associated with "pollen tube growth", "cell tip growth", and "developmental cell growth" ( Figure 3C, Table S8). DEMs in cluster 3 were related to "histone exchange" and "methionine adenosyltransferase activity" ( Figure 3C, Table S8). DEMs in cluster 5 were correlated with "large ribosomal subunit" and "cytosolic ribosome" (Figure

Target Gene Prediction via Degradome Sequencing
Globally, 16,705,004 sequences were obtained from degradome sequencing, and 16,697,634 clean reads sequences were obtained after fastp quality control. Q20 and Q30 sequences accounted for 98.4% and 43.8% of the clean read sequences, respectively (Table S5). Within the set of 421 miRNAs, there were 419 that targeted 3582 mRNA, forming 15,522 miRNA-mRNA modules (Table S6). Of these, 2942 modules were 'category 0', 1626 modules were 'category 1', and 6995 modules were 'category 2'. In the subset of 248 DEMs, 247 targeted 2069 genes, of which 1772 were annotated by SwissProt (Table S7). Out of 1772 annotated genes, 155 isoforms were annotated as transcription factor genes. The top five transcription factors with the highest number of genes were squamosa promoter-bindinglike protein (SPL, 28 genes), dehydrin (26 genes), growth-regulating factor (23 genes), homeobox-leucine zipper protein (15 genes), transcription factor GAMYB (15 genes), and NAC domain-containing protein (nine genes). On the other hand, 1617 isoforms were annotated as functional protein genes. The top five categories of functional protein genes, ranked by the number of genes, were disease resistance protein (245 genes), TMV resistance protein (103 genes), probable disease resistance protein (60 genes), 60S ribosomal protein (29 genes), and 40S ribosomal protein (19 genes).
Further, GO term enrichment analysis was performed on DEM target genes in each cluster to infer the mediatory role of DEMs during drought and rehydration. DEMs in cluster 1 were linked to GO terms such as "oxidoreductase activity using superoxide radicals as acceptor" and "superoxide dismutase activity" ( Figure 3C, Table S8). The DEMs in cluster 2 were associated with "pollen tube growth", "cell tip growth", and "developmental cell growth" ( Figure 3C, Table S8). DEMs in cluster 3 were related to "histone exchange" and "methionine adenosyltransferase activity" ( Figure 3C, Table S8). DEMs in cluster 5 were correlated with "large ribosomal subunit" and "cytosolic ribosome" ( Figure 3B, Table S8). DEMs in cluster 6 were linked to "cul3-RING ubiquitin ligase complex" and "ubiquitin ligase complex" ( Figure 3C, Table S8).

The Negatively Correlated miRNA-mRNA Modules
Taking advantage of the available transcriptomics data for the same samples from our previous study [46], negatively correlated DEM-target modules were identified through Pearson's correlation analysis. This analysis involved assessing the normalized DEM expressions (average RPM of three biological replicates) and target expressions (average FPKM of three biological replicates). Given the complexity of the network, which hindered effective visualization, a concise mini-miRNA network was generated using three criteria: (i) inclusion of degradome signals with quality categories zero, one, or two at the cleavage sites of the targets; (ii) selection of differentially expressed genes (DEGs) from the transcriptome as targets; (iii) and a DEM-target expression coefficient of r < −0.80. As a result, 100 negatively correlated DEM-DEG modules were identified, containing 30 DEMs and 30 DEGs (Table S9). For a more intuitive presentation, log 2 (fold change) corresponding to these DEMs and DEGs was used to visualize the expression profiles ( Figure 4A).   (DEMs) and differentially expressed genes (DEGs) under drought and rehydration. The heatmap on the left presents the expression of the DEMs, and the heat map on the right presents the expression of the DEGs. The red and blue dots represent DEMs and DEGs, respectively, and the solid black lines in between represent the degradome-validated targeting relationships between the DEMs and the DEGs. T-plots and miRNA-mRNA alignments represent pma-miR156a-5p cleaves isoform_118780 (B), pma-miR1312c-3p cleaves isoform_166904 (C), pma-miR166l-3p cleaves isoform_24257 (D). The red dots and triangles represent the cleavage positions on the target genes.

Discussion
Studies have demonstrated that miRNAs regulate drought response in crop plant species [22,48,49] as well as forest tree species [50]. Pine trees respond to drought through a complex process that involves the expression reprogramming of multiple genes [51,52]. A study conducted on P. pinaster showed that a significant number of miRNAs were involved in the drought response of pine trees [37]. Roots play a crucial role in a plant's response to drought [53][54][55]. However, the role of miRNAs in the roots of Masson pine trees in response to drought and rewatering has not been reported yet. Here, we employed high-throughput small RNA sequencing to identify miRNAs and their expression profiles from 21 libraries under control and stress conditions to study the effect of drought and rehydration in seedling root of P. massoniana. A total of 421 miRNAs were identified, among which 248 miRNAs exhibited differential expression under drought stress and rewatering. Through degradome sequencing, potential miRNA-mRNA regulatory modules were predicted.

Features of P. massoniana miRNA Population
In this study, each library generated a minimum of 19,969,315 clean reads, providing sufficient sequencing depth for subsequent analysis. The reads exhibited high quality, with Q20 base ratio > 99.50%, and the Q30 base ratio was >92.51%. The most abundant read lengths were 21 nt, followed by 20 nt and 22 nt, consistent with previous finding in P. massoniana [40,43] and Pinus tabuliformis [56]. Here, the 421 miRNAs were 20, 21, and 22 nt long, which is a characteristic commonly found in plant miRNAs [57]. Sequences of 21 and 22 nt were the most prevalent in both known and novel miRNAs ( Figure 1B). Moreover, the mature miRNA had a strong bias toward a 5 -uridine residue ( Figure S1), consistent with previous observation in pine [44,52]. These factors, namely the length and the preference of 5'-uridine residue, can influence the partitioning of miRNAs onto specific AGO proteins [58][59][60], thereby affecting the functions performed by these AGO proteins [57]. AGO1 is considered the most essential protein in the miRNA pathway [61,62] and it prefers to bind 21 nt miRNAs with a 5'-uridine residue [57,58].

miRNA Modules Mediate Translational Regulation in Drought Response
Ribosomal proteins (RPs) are components of ribosomes and perform multiple roles in biological processes, such as ribosome biogenesis, protein synthesis, cell growth, development, and abiotic stress response [69]. For example, RP genes could be induced by water deficit in rice root [70,71] and enhance the expression of two RP genes, respectively, as both resulted in improved drought tolerance in rice [70]. Moreover, knockdown of 60S ribosomal protein L14-2 resulted in reduced tolerance to drought stress in cotton [72]. Here, DEMs from five miRNA families: miR482, miR11524, miRN88, miR11476, and miR396 were found to target fifteen isoforms annotated as ribosomal protein (Table S8). Notably, the expression of pma-miRN88a-5p and pma-miRN88b-5p were negatively correlated with their target gene, isoform_34691 (RL222), by drought in D1 versus C1 and D3 versus C3 ( Figure 4A). These findings suggest that these miRNAs may regulate the ribosome and play a role in the drought response of P. massoniana.
Environmental stress activates unfolded protein response (UPR) in the endoplasmic reticulum (ER), a highly conserved response in plants [73,74]. Persistent UPR can cause programmed cell death, so UPR is under tight control [75,76]. ER-resident luminal-binding protein (BiP), a central UPR regulator [77,78], aids protein folding, re-establishing ER homeostasis [79][80][81]. BiP plays a vital role in drought tolerance. Overexpressing BiP improved drought tolerance in plants [82][83][84][85]. However, overproduced BiP proteins suppressed the expression of BiP, indicating a negative feedback mechanism of the UPR, whereby the cell may reduce nonessential BiP transcripts when functional BiP proteins are sufficient in protein folding during ER stress [86]. Similar observations have been made in yeast and mammalian cells, where the overexpression of functional BiP protein mitigates UPR [87,88]. In our study, pma-miR11571-5p was found to target four BiP genes and was upregulated by drought stress in D1 versus C1 and D3 versus C3 (Tables S4 and S8). The results indicate that pma-miR11571-5p was involved in the negative feedback regulation by suppressing the BiP's expression, thus maintaining the homeostasis of ER under drought stress.
Rapid responses to environment perturbation are vital for plants due to their sessile lifestyle [89]. Such responses, like signal transduction and cell cycle control, require prompt elimination of certain proteins, such as misfolded proteins or various normal short-lived regulators [89,90]. The ubiquitin proteasome system (UPS) mediates a major pathway responsible for protein degradation [91]. The UPS is initiated with a conserved cascade reaction involving E1, E2, and E3 enzymes, leading to the attachment of ubiquitin to specific proteins [89]. The most common E3 ligase in plants is the Skp1-Cullin-F-box (SCF) protein complex, which recognizes specific substrates through the binding interaction between SKP1-like ASK and F-box proteins [92][93][94]. Among the F-box protein family, F-box kelch proteins (FBKs) represent one of the largest subfamilies [93,95]. Previous studies in Arabidopsis [96] and wheat [97][98][99] have demonstrated that FBKs binds to ASK proteins and that the overexpression of FBK enhanced drought tolerance in plants. In our study, miR482 family members were found to target, with high confidence (category = 0), six genes annotated as FBK (Table S8). miR482 has also been shown to target F-box genes in lychee [100] and strawberry [101]. Therefore, our findings suggest that members of miR482 family played a role in drought response in P. massoniana.

miRNA Modules Mediate Cell Wall Modification in Drought Response
S-adenosylmethionine synthase (SAMS) catalyzes the synthesis of S-adenosylmethionine (SAM) from methionine and adenosine triphosphate (ATP) [102]. SAM is involved in multiple transmethylation reactions, including those related to lignin biosynthesis [103,104]. Increasing deposition of lignin in cell walls may be one of the mechanisms by which cells respond to drought [105][106][107]. Methylation of lignin precursors is a critical step in lignin synthesis, with SAM acting as the primary methyl group donor [108,109]. Previous studies on Pinus banksiana [110], peanut [111], soybean [112], and cucumber [113] have shown that the expression levels of SAMS protein and/or transcript in roots were responsive to drought stress. In this study, members from miR396 and miRN74 were found to target eight isoforms annotated as SAMS (Table S8). These genes were enriched in two GO pathways, "methionine adenosyltransferase activity" and "S-adenosylmethionine biosynthetic process" (Figure 3D, Table S8). Their findings suggest that miR396 and miRN74 was involved in lignin biosynthesis during drought in P. massoniana by mediating the expression of SAMS.
HOX32 encodes a transcription factor that belongs to the HD-ZIP III group [114]. In our study, isoform_242457, annotated as HOX32, was found to be targeted by pma-miR166l-3p ( Figure 4A, Table S9). A previous study on rice has demonstrated that miR166 targets OsHOX32, and knockdown of miR166, or overexpression of OsHOX32 led to a reduction in lignin content in cell wall [114]. Many studies have shown that drought increased lignin accumulation in the roots [105][106][107]. However, since lignin biosynthesis consumes a high and irreversible input of carbon sources, its deposition is tightly regulated through transcriptional, posttranscriptional, and posttranslational processes [114][115][116]. In our study, the pma-miR166l-3p were downregulated, while the expression of its target gene HOX32 were upregulated in D2 versus C2 (rewatering) and D3 versus C3 (drought) stages. These results suggest that the pma-miR166l-3p:HOX32 module may be involved in fine-tuning of lignin biosynthesis under drought and rehydration. HOX32 was slightly downregulated in D1 versus C1 (drought). This may be due to other factors regulating its expression [117].

miRNA Modules Mediate ROS Scavenging in Drought Response
ROS are free radicals of oxygen. They may have both beneficial and harmful effects [118]. During drought, ROS can accumulate excessively in cells, leading to oxidative stress [119]. To counteract this, cells possess various enzymes, such as superoxide dismutase (SOD), which provide the first line of defense against oxidative stress [120]. The copper-zinc superoxide dismutase (CSD) is the most common type of SOD [121]. Previous studies have shown that miR398 regulates drought tolerance in plants by targeting CSD genes [122][123][124]. Further, miR398 was found to be downregulated in the roots by drought stress in pea [68] and legume [123]. Overexpression of miR398 reduced the expression of CSD and impaired plant drought tolerance [124,125], while the knockdown of miR398 increased CSD expression and enhanced plant drought tolerance [124]. In this study, pma-miR398a-3p, pma-miR398b-3p, and pma-miR398c-3p targeted three genes annotated as CSD ( Figure 4A, Table S8). The expression of these pma-miR398 was downregulated in D1 versus C1 and D3 versus C3 (drought). These findings suggest that the downregulation of pma-miR398 make root cells to produce more SOD enzymes, enabling them to scavenge peroxides induced by drought stress.

Putative miRNA-Mediated Regulatory Network
Based on DEM-target correlations, a schematic model was proposed for miRNAmediated regulatory network during drought and rehydration ( Figure 4A, Table S9). In this network, pma-miR166 targeted the most genes among the negatively correlated miRNA-mRNA modules ( Figure 4A). In addition to HOX32 as previous mentioned, pma-miRNA166 members were found to target Early Responsive to Dehydration stress (ERD) genes such as ERD10. ERD10 is a dehydrin protein [126], which is a potent chaperon that activates other protective proteins or acts as a plasma membrane stabilizer to protect cells [127,128]. The expression of ERD10 can be rapidly increased by dehydration [126] and erd10 mutants show reduced drought tolerance [129]. Here, pma-miR166a-5p and six other members were downregulated, while ERD10 was highly upregulated by drought in D1 versus C1 and D3 versus C3. This result implies that reduced pma-miR166 expression promoted ERD protein expression in response to drought. Noticeably, pma-miRN89-5p, which was substantially downregulated by drought in D1 versus C1 and D3 versus C3, ( Figure 4A) targeted another ERD10 gene. This result indicates that pma-miRN89-5p plays a role in drought response. Moreover, pma-miR156a-5p and pma-miR482b-3p targeted GRP genes, including cell wall-associated GRP (isoform_292633 and isoform_33236) and RNA-binding GRP (isoform_32791, isoform_46586 and isoform_77096; Figure 4A). Cell wall-associated GRPs have been reported to be involved in cell elongation [130] and root size control [131]. As for RNA-binding GRP, they may play a role in RNA stabilization, processing, and transport according to a previous report [132]. miR396 plays a pivotal role in regulating plant architecture through its mediation of gibberellin (GA) signaling [133,134]. GA signaling is subjected to regulation by hormone transporters, such as NITRATE TRANS-PORTER1/PEPTIDE TRANSPORTER (NPF) [135,136]. NPF3 overexpression dramatically inhibited root growth [135]. Moreover, the suppression of GA activity resulted in improved plant drought tolerance [137,138]. In this study, pma-miR396 were found to target one gene annotated as Protein NRT1-PTR FAMILY 5.3 and PTR4 (Table S9), which is encoded by NPF5.3 gene ( Figure 4A). And, here, NPF5.3 was downregulated by drought treatment ( Figure 4A). This result indicated that the miR396:NPF5.3 module was involved in regulating GA signaling during drought in P. massoniana. pma-miR156a was found to target SPL12 ( Figure 4A). The miR156:SPL module exists in multiple species and was involved in drought response and growth regulation [139]. Thus, the pma-miR156a-5p:SPL12 module may play a role in drought response in P. massoniana. In addition to targeting RL222, pma-miRN88a-5p and pma-miR88b-5p also targeted UBI1P, suggesting a role in regulating E1 activating enzyme in UPS. pma-miRN86-3p targeted LPR2. LPR was reported to mediate the response of root meristems during phosphate availability [140]. pma-miR1312b-3p and pma-miR1312c-3p were found to target two genes (isoform_ 166,904 and isoform_95799, Figure 4A), which were annotated as a target of AvrB operation (TAO1, Figure 4A, Table S9). TAO1 was reported to play a key role in signaling during the response to pathogens [141]. pma-miR1312 might be involved in drought response at the signaling level. In addition, miR482, miR950, and miRN90 targeted three genes, isoform_41247, isoform_264908, and isoform_76607, respectively, which currently have no SwissProt annotation. They may be novel genes in response to drought stress.

Plant Materials and Stress Treatments
The materials originated from P. massoniana seedlings of a full-sibling family (code: 19-309) derived from two high resin-yielding parents, GZ001 (female) and GZ549B (male). Approximately one month after germination, young seedlings were randomly allocated into three treatment groups and four control groups with each group being repre-sented by three biological replicates. Subsequently, the seedlings were transplanted into non-woven pots (4 cm in diameter and 8 cm in depth) containing a mixture of coconut husk, loess soil, and peat in a ratio of 6:3:1 (v/v/v). During the drought experiment, all the sibling seedlings were cultivated under constant conditions of 26 • C, 16 h/8 h light/dark cycle, 60% relative humidity, and 80 µmol m −2 s −1 photon flux inside an RXZ-1000A-LED growth chamber (Ningbo Jiangnan Instrument Factory, Ningbo, Zhejiang, China). The experiment involved three treatments, namely, D1, which entailed withholding irrigation for seven days upon needle wilting; D2, which involved withholding irrigation for seven days and 17 h followed by rewatering and another seven-hour period of no irrigation; and D3, which entailed withholding irrigation for eight full days. The four controls C0, C1, C2, and C3 corresponded to the time points of treatment start, D1, D2, and D3, respectively, with regular watering every other day. To minimize the effects of circadian rhythm, sampling was conducted for half an hour at the same time each day. Approximately 0.5 cm from the root tip was collected and snap-frozen in liquid nitrogen. A total of 21 samples were obtained for small RNA sequencing, including three treatment groups and four control groups with three biological replicates each.

Small RNA Sequencing
RNA extraction was performed using the Invitrogen TRIzol ® Plus RNA Purification Kit (Thermo Fisher Scientific Inc., Waltham, MA, USA). The quality and quantity of RNA samples were evaluated using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific Inc.) and an Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). sRNA-seq libraries were constructed and sequenced at Beijing Genomics Institute (BGI, Shenzhen). Briefly, 1µg total RNA from each sample was separated using polyacrylamide gel electrophoresis (PAGE) and the 18-30 nt stripes were selected and recycled. The 3 and 5 adaptors were ligated to recycled RNA. The adapter-ligated RNA was reversetranscribed and the cDNA product was amplified by PCR for 16 cycles. Amplified products were purified, and then quantified using a Qubit fluorometer (ThermoFisher Scientific, Massachusetts, USA). Subsequently, the purified products were used to produce singlestrand circle DNA as the final cDNA library. A total of 21 libraries, namely, C0_1, C0_2, C0_3, C1_1, C1_2, C1_3, C2_1, C2_2, C2_3, C3_1, C3_2, C3_3, D1_1, D1_2, D1_3, D2_1, D2_2, D2_3, D3_1, D3_2, D3_3, were constructed and sequenced on a BGISEQ-500 platform (BGI, Shenzhen, Guangdong, China) using 50 bp single-end read chemistry.

De Novo Prediction and Annotation of miRNAs
The adaptor for sRNA sequencing was obtained using dnapi.py [142]. Cutadapt [143] was used to perform quality control (-j 0 -a adaptor -quality-base 33 -m 18 -M 30 -O 4discard-untrimmed -q 20 -max-n 0). Reads that contained more than one N base (also known as ambiguous base) and base with quality value less than 20 were discarded. Then, the obtained clean reads were mapped against the Rfam (v12.1) database, with one base mismatch permitted, to remove rRNA, tRNA, snRNA, and snoRNAs. Such reads with counts less than 10 were also discarded. Then, reads were analyzed by miREvo [144] and miRDeep-P2 [145] for identifying known miRNA as well as predicting novel miRNAs. Briefly, reads with more than 15 mapping sites in the reference transcriptome were filtered out. The length of the potential miRNA precursor should not have exceeded 300 nt. The range of the length of mature miRNA and miRNA* should have been 20-24 nt. miRNA/miRNA* duplexes should not have contained a large loop and should contain up to five mismatches.

Degradome Library Construction and Target Gene Prediction
Equal amounts of RNA from both control (C0, C1, C2, C3) and drought (D1, D2, D3) treatments were pooled to construct one degradome library. The libraries were constructed following the method proposed by Fang et al. [151] and the sequencing was performed on the Illumina NovaSeq 6000 platform with 50 bp single-end chemistry at Genedenovo Biotechnology Co., Ltd. (Guangzhou, Guangdong, China). The degradome sequences were subjected to quality control using fastp (v 0.23.2) [152]. The resulting high-quality reads were then used to predict target genes using Cleaveland4 (parameters: -t -c 2) [153] with a full-length transcriptome (SRA accession number: PRJNA667166) as reference sequence from our previous study [45].

Validation of miRNA Expression via qRT-PCR
For the validation of miRNA expression obtained from sRNA-seq, the abundance of ten mature miRNAs was quantified via qRT-PCR. cDNAs of each of the 21 samples were synthesized using a Mir-X™ miRNA First Strand Synthesis Kit (TaKaRa, Dalian, China). Using the ACT1 gene as an internal control, the cDNAs were verified by qRT-PCR with the corresponding mature sequences of the miRNAs used as forward primers (Table S10), and the universal miR 3' primer, included in the Mir-X™ miRNA First Strand Synthesis Kit, as reverse primer. Reactions with three replicates for each of the samples were performed on ViiA™ 7 Real-Time PCR Systems (Applied Biosystems, Waltham, MA, USA). The relative expression data were calculated using 2 −∆∆CT method [154]. The qRT-PCR results are shown in Figure S2.

Conclusions
In this study, 21 sRNA-seq libraries were sequenced to identify miRNAs in P. massoniana seedling roots under drought and rewatering conditions. A total of 421 miRNAs were identified and among them, 248 miRNAs were differentially expressed. The Gene Ontology enrichment analysis of predicted target genes of the differentially expressed miRNAs indicated their participation in drought response via various mechanisms such as translational and posttranslational regulation, cell wall modification, and ROS scavenging. miRNAs, such as miR482, miR398, miR11571, miR396, miR166, miRN88, and miRN74, along with target genes, such as those encoding F-box/kelch-repeat protein, 60S ribosomal protein, copper-zinc superoxide dismutase, luminal-binding protein, Sadenosylmethionine synthase, and Early Responsive to Dehydration Stress, could potentially have vital functions in responding to drought conditions. miRNA-mRNA modules, such as pma-miR396b-3p:NPF5.3, pma-miR156a-5p:SPL12 and pma-miR1312b-3p:TAO1, could also play important roles in drought responses in Pinus. This study presents a valuable resource for further molecular investigation on complex regulatory network of gene expression and uncovering new players functioning in drought tolerance in Pinus.  Data Availability Statement: sRNA-seq and degradome sequencing data in this paper have been deposited in the Genome Sequence Archive [155] in National Genomics Data Center [156], China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (GSA: CRA012245) and are publicly accessible at https://ngdc.cncb.ac.cn/gsa (accessed on 15 August 2023).