Genome-Wide Identification and Analysis of Small Nucleolar RNAs and Their Roles in Regulating Latex Regeneration in the Rubber Tree (Hevea brasiliensis)

Small nucleolar RNAs (snoRNAs) are a class of conserved nuclear RNAs that play important roles in the modification of ribosomal RNAs (rRNAs) in plants. In rubber trees, rRNAs are run off with latex flow during tapping and need to be regenerated for maintaining the functions of the laticifer cells. SnoRNAs are expected to play essential roles in the regeneration of rRNAs. However, snoRNAs in the rubber tree have not been sufficiently characterized thus far. In this study, we performed nuclear RNA sequencing (RNA-seq) to identify snoRNAs globally and investigate their roles in latex regeneration. We identified a total of 3,626 snoRNAs by computational prediction with nuclear RNA-seq data. Among these snoRNAs, 50 were highly expressed in latex; furthermore, the results of reverse transcription polymerase chain reaction (RT-PCR) showed the abundant expression of 31 of these snoRNAs in latex. The correlation between snoRNA expression and adjusted total solid content (TSC/C) identified 13 positively yield-correlated snoRNAs. To improve the understanding of latex regeneration in rubber trees, we developed a novel insulated tapping system (ITS), which only measures the latex regenerated in specific laticifers. Using this system, a laticifer-abundant snoRNA, HbsnoR28, was found to be highly correlated with latex regeneration. To the best of our knowledge, this is the first report to globally identify snoRNAs that might be involved in latex regeneration regulation and provide new clues for unraveling the mechanisms underlying the regulation of latex regeneration.

However, several previous studies have shown that many snoRNAs known as "orphan snoRNAs" have typical characteristics of snoRNAs but lack complementarities to snRNA, rRNA, tRNA, or other reported RNAs (Bachellerie et al., 2002). In addition to canonical modification functions, some snoRNAs were involved in the process of pre-mRNA alternative splicing. For example, compelling evidence has been gathered for the SNORD115 family of snoRNAs, which were found to be deleted in individuals with Prader-Willi syndrome (PWS), to direct alternative splicing of the serotonin receptor subtype (5-HT2CR) pre-mRNA (Kishore and Stamm, 2006;Doe et al., 2009;Kishore et al., 2010). In addition, snoRNAs can be further processed into small ncRNA (sno-miRNAs), which are almost similar in function to the miRNAs (Lee et al., 1993;Ender et al., 2008;Kawaji et al., 2008;Taft et al., 2009;Burroughs et al., 2011). Understanding the emerging relationship between snoRNAs and RNA silencing may reveal more about the role of RNA in gene regulation (Lui and Lowe, 2013). Additionally, snoRNAs regulate chromatin assembly and are involved in cellular response to stress (Nicoloso et al., 1996;Chow and Heard, 2009;Michel et al., 2011;Schubert et al., 2012). Thus, these findings suggest that the potential roles of snoRNAs should be extensively investigated.
To date, snoRNAs have not been extensively studied in plants, and only a few studies have reported snoRNAs in Arabidopsis thaliana, Oryza sativa (rice), and other plants. To date, 295 and 544 snoRNA genes were discovered in A. thaliana and rice, respectively (Liang-Hu et al., 2001;Chen et al., 2003;Liu et al., 2012). Recently, the emergence of high-throughput next-generation sequencing (NGS) technology has facilitated the rapid identification of ncRNAs, such as snoRNAs, and the determination of their expression profiles at a genome-wide level in different species. Global identification and analysis of snoRNAs in rice showed 433 snoRNA candidates, including 125 previously un-annotated snoRNAs based on their structural features and common motifs (Liu et al., 2012). In maize, high-throughput sequencing was performed to identify a total of 169 ncRNAs, including 70 snoRNAs. Moreover, target site analysis revealed that 22 snoRNAs can guide to 38 sites of 2 -O-methylation and pseudouridylation modification of rRNAs and snRNAs (Li et al., 2018). A combination of high-throughput sequencing and massive analysis of complementary DNA (cDNA) ends (MACE) was used to identify the complete tomato pollen sncRNAome, namely, miRNAs, snoRNAs, and tRNAs, and the regulatory network of these small RNAs was integrated into heat stress response in different stages of the development of tomato pollens (Bokszczanin et al., 2015).
The rubber tree (Hevea brasiliensis) is a tropical perennial crop originating in the Amazon basin and is the most economically important tropic crop in many tropical countries that provides the largest renewable natural rubber (NR) for the industry (Lau et al., 2016). The rubber tree latex is the cytoplasm of laticifers, and the latex exudes when the bark is tapped. A major part of the latex consists of rubber particles, lutoids, and other cytoplasmic components; in addition, it contains a small amount of laticifer nuclei (Shearman et al., 2014). High turgor pressure in the laticifers results in the release of the cytoplasm from the laticifers (Dickenson, 1965). However, biotic and abiotic stresses, namely, drought, salinity, temperature, mechanical wounding, and pathogen attacks, are the main factors affecting the production of latex (Jaspers and Kangasjärvi, 2010). During tapping, rRNA and proteins are run off from the laticifers and need to be re-synthesized before the next tapping to ensure the continuous yield of latex (McMullen, 1962). rRNAs are the fundamental component of ribosomes, which are the machinery for protein synthesis (Yamashita et al., 2016). Therefore, rRNA maturation, during which proteins required to maintain the biological activity in laticifers are synthesized, is considered to be the first step in latex regeneration. rRNAs are first transcribed as pre-rRNAs, and only after 2 -O-methylation or pseudouridylation modifications they are processed into mature rRNAs (Tupý, 1988). Therefore, snoRNAs are expected to play important roles in latex regeneration. To date, a few snoRNAs have been identified in H. brasiliensis, and their functional roles in H. brasiliensis are largely unknown (Gébelin et al., 2012); therefore, deep sequencing should be performed to identify and analyze the snoRNAs and will contribute to a better understanding of the mechanism underlying the regulation of their expression (Chow et al., 2007).
In this study, we constructed and sequenced a non-coding RNA population ranging from 50 to 300 nt from the rubber tree leaves and latex samples for global screening and identifying potential functional snoRNAs in the rubber tree. We identified a total of 3,626 snoRNAs by this extensive analysis, and showed that 50 of them were highly expressed in latex; the expression of 31 out of the 50 snoRNAs in latex was confirmed by reverse transcription polymerase chain reaction (RT-PCR). The results of the analysis of the correlation between snoRNA expression and latex yield showed that 13 laticifer-abundant snoRNAs were highly associated with latex yield and total solid content (TSC). HbsnoR28 was highly correlated with latex regeneration in the depletion test. Thus, our results improve the understanding of latex metabolism and regeneration in the rubber tree.

Preparation of the Plant Material
Five-year-old mature rubber trees (H. brasiliensis cv. Reyan7-33-97) were maintained at the experimental plantation of the Rubber Research Institute, Chinese Academy of Tropical Agricultural Sciences (Danzhou, Hainan, China). The leaves from the first umbrella were harvested, and the latex was collected using the conventional S/2 d/3 (half spiral cut tapped at the third daily frequency) tapping system. Samples from 10 tapped trees were combined into one biological replicate. Each experiment was conducted using at least three replicates. Total latex production was measured using the conventional approach of S/2 d/3 and was weighed as TSC after air drying at 60 • C.

Isolation of Nuclear RNA and Construction of Complementary DNA Libraries of Non-coding RNAs
The leaves were cut into 0.5-to 1-cm-wide strips. The samples were treated with cellulase (R-10; Solarbio R , Beijing, China) to produce protoplasts (Tang et al., 2007). The protoplasts were washed five times with a washing buffer (154 mM NaCl, 125 mM CaCl 2 , and 5 mM KCl) and preserved at −80 • C. The nuclei were isolated using a nuclear isolation buffer [NIB; 10 mM N-2-hydroxyethylpiperazine-N-2-ethane sulfonic acid (HEPSE), 0.8 M sucrose, 5 mM MgCl 2 , 5 mM ethylenediamine tetraacetic acid (EDTA), pH 7.6, and then 0.07% β-mercaptoethanol after autoclaving]. The nuclei of laticifer cells present in the latex were centrifuged at 16,000 rpm for 2 h under 4 • C using an ultracentrifuge (Sovall TM LYNX 4000; Thermo Fisher Scientific TM , United States). We collected only the precipitate containing most of the lutoid and a small amount of laticifer nuclei. The pellet was washed two times by violent vortexing and resuspending in the NIB. Then, the TRIzol method (Invitrogen TM , United States) was used to isolate the nuclear RNAs according to the protocol of the manufacturer. The concentrations of samples were measured using an ultra-microspectrophotometer (NanoDrop2000; Thermo Fisher Scientific TM , United States). Furthermore, rRNAs were isolated using Ribo-Zero rRNA Removal Kits (Plant Leaf, MRZPL116; Illumina R , United States) according to the standard protocol.
The sequencing libraries were constructed as previously described (Deng et al., 2006). Briefly, RNA fractions of 50-500 nt in length were isolated by polyacrylamide gel electrophoresis (PAGE) and sequentially ligated to the 5 and 3 end adaptors. Then, the adapter-ligated RNAs were reverse transcribed using primers with partial complementarity to the adaptors. The cDNAs with both end adaptors were subsequently sequenced using HiSeq2500 (Illumina, United States).

Data Analysis and de novo Transcriptome Assembly
Raw reads were trimmed using Trim_Adaptor script 1 with the minimum clip length fixed at six nucleotides. The clean paired-end reads were assembled using the Trinity program (Trinityrnaseq_r20131110) with the following parameters: seqType, fastaq; JM, 500G; minimum k-mer value, 1; CPU, 32 (Grabherr et al., 2011). The clean reads were mapped onto the reference genome of the rubber tree (Tang et al., 2016) with Tophat2 v2.1.0. to calculate the mapping location using default parameters. For differential expression calculation, RSEM

Identification of Small Nucleolar RNAs and Prediction of Modification Sites in the Ribosomal RNAs Sequence
Small nucleolar ribonucleic acid candidates were identified using snoScan 0.9b (Lowe and Eddy, 1999) and snoSeeker v1.1 (Yang et al., 2006) on the basis of the assembled transcriptome. On the basis of the known modification sites in A. thaliana, O. sativa, and other model organisms, we predicted the putative modification sites on the rRNA sequence according to consensus rRNA sequence between the rubber tree and the model plant organisms. Guide snoRNA candidates were identified by FIGURE 1 | Ribonucleic acid (RNA) sequencing of the nuclei prepared from (A) latex and (B) leaves of the rubber tree. The reads were mapped onto the rubber tree genome and (C) percentages of the mapped reads on each region were calculated. As the read is 150 nt, one reads may map onto several regions.

Identification of Putative Small Nucleolar RNA Target Sites
Methylation sites in the rRNAs were computationally predicted using the CDSeeker program by scanning for complementary regions of the snoRNAs from corresponding 5.8S, 18S, and 28S sequences in H. brasiliensis. The first nucleotide of the D or D' box was included in the guide sequence, because this position may participate in the formed duplex. One mismatch, two GU base pairs, and no bulges were allowed in the duplex between the guide and the rRNA region containing the target sequence. Typically, the H/ACA box snoRNAs are distinguished by the presence of an ACA motif at the 3 -end and a common hairpin-hingehairpin-tail secondary structure with the H box (ANANNA) in the hinge region.

Analysis of the Expression of Small Nucleolar RNAs
Total RNAs were extracted from the leaves and latex using the Trizol method. We reverse-transcribed 2 µg of total RNA to cDNA using a SuperScript III reverse transcriptase kit (Invitrogen, United States) according to the protocol of the manufacturer. All the expression analyses were conducted with three biological replicates. HbI5I was used as the internal reference for gene expression analyses. Semi-quantitative PCR was performed using the following conditions: 95 • C for 5 min, followed by 30 cycles of denaturing at 95 • C for 15 s, annealing at 60 • C for 30 s, and elongation at 72 • C for 30 s. About 10 µl of the PCR product was analyzed by agarose gel electrophoresis. Quantitative RT-PCR (qRT-PCR) was performed to detect the expression level of each snoRNA in the hybrid population. Specific primers for respective snoRNAs (Supplementary Table 4) were designed using Primer 5, and all of the amplified fragments were sequenced for target verification. The qRT-PCR reaction was performed on the LightCycler 2.0 System (Roche Diagnostics, Penzberg, Germany) using the SYBR Green premix kit (TaKaRa, Dalian, China) according to the instructions of the manufacturer. The relative abundance of transcripts was calculated using the Light Cycler Relative Quantification Software. Spearman's correlation was performed with adjusted TSC (TSC/C) and qRT-PCR results using R scripts.

Latex Regeneration Using an Insulated Tapping System
We designed an insulated tapping system (ITS) in this study (Figures 6A,B). Briefly, a 50-cm-tall bark was insulated by girding the trunk at a height of 50 and 100 cm. The girding was carefully performed to cut the yellow layer, where, the laticifers  are abundant while causing no damage to the inner watery layer, which transports nutrients. In the ITS, the tapping will exude the latex only from the dedicated 50-cm bark, and the transportation of nutrients will not be affected.

Nuclear RNA Sequencing and Data Analysis
To identify snoRNAs globally in the rubber tree, we constructed two cDNA libraries using the total RNAs extracted from the nuclei of the latex ( Figure 1A) and leaves (Figure 1B) of the rubber tree. Then, RNA-seq was performed and 4-Gb clean paired-end reads data were obtained for each sample.
A total of 27 million raw reads were obtained for the leaf samples and 21 million for the latex, with the Q30 of each sample higher than 90% (Supplementary Table 1). After filtering adapter sequences and discarding low-quality reads, the reads were mapped to the draft rubber tree genome (Liu et al., 2019). The mapping results indicated that the majority of sequences were in the exon and intergenic regions (Figure 1C), demonstrating that the RNAs extracted from the crude nuclei contained abundant amounts of mRNA and pre-mRNA. The nucleic transcriptome was then de novo assembled using the Trinity program (Grabherr et al., 2011). To reduce redundancy, only isoforms that had the highest expression level within each subcomponent were selected. Then, contigs that match with the rubber tree structural genes and EST sequences were removed. Finally, a total of 229,669 contigs representing unique transcripts were obtained, with an average length of 1,499 nt and N50 of 1,188 nt. These contigs were used as the ncRNA transcript repertoire from the rubber tree latex and leaf for further analysis.

Computational Identification of Small Nucleolar RNAs
To identify snoRNA candidates, the rRNA modification sites of 2 -O-methylation and pseudouridylation were first predicted  Figure 1). These rRNA modifications were further used as the input for C/D and H/ACA subfamily snoRNA computation according to the pipeline shown in Figure 2.
For snoRNA computation, the ncRNA transcript repertoire was analyzed with the snoScan (Lowe and Eddy, 1999) and snoSeeker programs (Yang et al., 2006). A total of 22,125 snoRNA candidates were predicted, namely, 11,269 and 10,856 identified using snoScan and snoSeeker, respectively. Then, these candidate snoRNAs were subjected to a series of filtering steps (Figure 2) to remove any possible contaminations and the sequences belonging to coding genes. Subsequently, 3,626 snoRNA candidates were retained, in which 1,712 were guided C/D snoRNAs, 86 guided H/ACA snoRNAs, and 1,828 were orphan snoRNAs ( Table 2).

Identification of Latex-Abundant Small Nucleolar RNAs
To elucidate the potential regulatory roles of snoRNAs in latex regeneration, we analyzed the expression of snoRNAs in the latex and leaves. The expression profiles of 3,626 candidate snoRNAs were calculated and the differentially expressed (DE) snoRNAs were identified. Because of differences in the abundance of total reads between the leaves and latex samples, the pre-rRNA, which only exists in the nucleus, was used as the reference for normalization of snoRNA expression. Our results showed that 50 snoRNAs were highly expressed in latex ( Figure 3A). Among the latex-abundant snoRNAs, 22 belong to the C/D subfamily and 28 belong to the H/ACA subfamily ( Table 3).
Semi-quantitative RT-PCR was performed to verify the expression levels of laticifer-abundant expressed snoRNAs identified from RNA-seq data analysis. The pre-rRNA transcript was used as the reference gene because mature rRNAs are not present in the nucleus. The amplification fragment of pre-rRNAs was designed to span the internal transcribed spacer 1 (ITS1). Among the 50 latex-abundant snoRNAs, 31, including 12 C/D snoRNAs (12 guides) ( Figure 3B) and 19 H/ACA snoRNAs (1 guide, 18 orphans) ( Figure 3C) were highly expressed in laticifer compared with leaves ( Table 2). The RT-PCR results further confirmed that these snoRNA were present abundantly in the latex (Figure 3).

Prediction of the Target Sites of Latex-Abundant Small Nucleolar RNAs
To elucidate the functions of snoRNAs in the rubber tree, we predicted their potential target sites according to sequence recognition rules. Previous studies have indicated that snoRNAs have a stable secondary structure and that their modification sites (2 -O-methylation sites and pseudouridylation sites) are usually conserved in eukaryotes, such as A. thaliana, O. sativa, and other model organisms (Watkins and Bohnsack, 2012;Dupuis-Sandoval et al., 2015). The 13 latex-abundant guide snoRNAs (12, C/D subfamily member and 1, H/ACA subfamily member) were found to target 31 specific modification sites (Figure 4). For 18S rRNAs, there are 12 modification sites that are targeted by seven-latex abundant snoRNAs, whereas for 28S rRNAs, the numbers are 19 and 10. Except HbsnoR23, which was predicted to target the pseudouridylation of U2790 of 28S rRNA, all the other 12 snoRNAs belong to the C/D subfamily and direct the 2 -O-methylation modification of 30 sites on 18S and 28S rRNAs ( Table 4).

Identification of Small Nucleolar RNAs Involved in Regulating Latex Regeneration
To identify snoRNAs that are related to latex yield potential, a correlation analysis was conducted between the latex yield and the expression of snoRNAs in an F1 population of Reyan66-2 × PR107. We randomly selected a total of 50 F1 individual trees, and the adjusted yields were ranked. In the rubber tree, a larger girth means a longer tapping cut and usually produces more latex, because more laticifers are present in the bark with a larger circumference. Therefore, we estimated the potential of latex production by normalizing the TSC by dividing with the trunk circumference (TSC/C), which will avoid the estimation bias caused by differences in trunk size. A Spearman correlation was performed to assess the relationship between the expression of snoRNAs ( Figure 5A) and the ranked TSC/C ( Figure 5B).
The results showed that 13 of the 26 laticifer-abundant snoRNAs were positively correlated with latex yield potential (TSC/C) (R 2 > 0.2, Figure 5C and Supplementary Table 3). Among the snoRNAs that were correlated with latex yield, HbsnoR1, HbsnoR23, and HbsnoR46 showed high correlation coefficients (R 2 > 0.4, Figures 5D-F). HbsnoR23 and HbsnoR1 belong to the guide H/ACA subfamily, while HbsnoR46 is an orphan H/ACA snoRNA. This result suggested that orphan H/ACA snoRNAs might participate in the regulation of latex production through pathways other than guiding rRNA modifications. The rubber tree laticifers are connected as a network structure that enables the latex to flow for a long distance when the bark is tapped. Many factors influence the movement of the latex in the laticifer network, such as laticifer turgor and plugging index. To minimize the replenish effect of long-distance movement of the latex, a 50-cm-tall bark was insulated by girding the trunk at a height of 50 and 100 cm. The girding was carefully conducted to cut the yellow layer, where the laticifers are abundant while causing no damage to the inner watery layer, in which the nutrients are transported (Figures 6A,B). Thus, this ITS will exude the latex only from the dedicated 50-cm bark, and the transportation of the nutrients will not be affected.
First, we evaluated the minimum time required for latex regeneration in the rubber tree. We selected five individual trees for tapping using the ITS. The trees were continuously tapped at intervals of 6, 12, 24, 48, and 72 h, and the TSC was calculated. The results showed that latex can be completely regenerated as quickly as 6 h ( Figure 6C). To further investigate the dynamics of latex regeneration, the laticifers in the ITS were isolated by continuous tapping until no latex exuded. After 12 h of regeneration, the trees were tapped at 30-min intervals, and the TSC was calculated. The results indicated that the TSC at each tapping dropped sharply in the first six tappings and then reached a stable state ( Figure 6D). After 3 h, the TSC of each tapping was maintained at a relatively stable level, suggesting that it was regenerated in the previous 30 min. We, therefore, considered the TSC at 3 h as the latex regeneration ability and investigated the correlation with snoRNA expression at this time point. We investigated the expression of 13 snoRNAs positively correlated with latex yield by Q-PCR and calculated their correlation with adjusted TSC (TSC/C). The results showed that one snoRNA (HbsnoR28) displayed a high correlation coefficient (R 2 = 0.6039) with TSC/C ( Table 5). The examination of HbsnoR28 expression showed that its expression was downregulated and highly consistent with TSC/C during continuous tapping, suggesting the regulatory roles of this snoRNA in latex regeneration ( Figure 6E).

DISCUSSION
Compared with studies on snoRNA in animals, those on snoRNAs in plants have not been extensively performed. Previously, 295 and 544 snoRNAs have been identified in plants such as A. thaliana (Kim et al., 2010) and O. sativa (Liu et al., 2012), respectively. H. brasiliensis is one of the most important industrial crops; however, a snoRNA has not been extensively studied in H. brasiliensis thus far. To date, only five snoRNAs have been annotated in a microRNA transcriptome from the leaf, bark, and root tissues (Gébelin et al., 2012). NcRNA sequencing techniques are powerful tools for the identification of snoRNAs. The RNA sequencing of the H. brasiliensis ncRNA library produced 229,569 sequences, and 3,626 of them were putative snoRNA candidates. After a further DE expression analysis and RT-PCR validation, 31  of them were identified to be latex-abundant snoRNAs. To the best of our knowledge, the global investigation of laticifer snoRNAs in the rubber tree have been reported for the first time in this study, and our results may be very important in understanding the mechanism underlying latex metabolism and regeneration. The functions of canonical snoRNAs are guiding the modification of rRNAs, snRNAs, tRNAs, or known stable RNA molecules. The most prevalent modifications found in rRNA molecules are the methylation of the ribose moiety at the 2 -hydroxyl group (2 -O-methylation) and the conversion of uridine to pseudouridine ( ). In this study, 31 putative modification sites were targeted by 13 laticifer-abundant guide snoRNAs in H. brasiliensis. Modification sites on the rRNA are mostly concentrated on the functionally important regions of the ribosome according to the 3D structural map, which would influence the activities of the ribosome. The putative modification sites on rRNA in H. brasiliensis may provide a clue for understanding the mechanism of how snoRNAs direct the modification of rRNA and regulate latex regeneration. In addition, 18 laticifer-abundant orphan snoRNAs were identified. Currently, no target RNA molecules are reported; however, it could be suggested that they might be involved in the regulation of physical processes, such as RNA silencing, rRNA chaperones, alternative splicing, and chromatin packaging.
To date, a limited number of studies have investigated snoRNA functions associated with physiological processes in plants. Loss of three snoRNAs (U33, U34, and U35a) in the rpL13a locus was sufficient to induce resistance to oxidative and lipotoxic stresses in vitro and prevented the spread of oxidative stress in vivo (Michel et al., 2011). In the conventional tapping system, large amounts of proteins and rRNAs are lost with the latex exudate and need to be replenished before the next tapping. rRNA was the first to be synthesized and is required for peptide translation. The regeneration of RNA will directly affect latex regeneration. Our results showed that 13 snoRNAs were positively correlated with latex production. Among these positively correlated snoRNAs, HbsnoR28 (H/ACA subfamily) showed a high correlation (R 2 = 0.6039) with latex regeneration potential in our ITS. This 138-nt-long orphan H/ACA snoRNA is novel and displays no homology with known snoRNAs in other species. Characterization of the functions of snoRNAs will be of great significance to examine latex metabolism and regeneration.
In this study, we developed a novel ITS, which separates the laticifers in dedicated regions from the laticifer network of the whole trees. The purpose of the ITS is to eliminate the long distance transportation of the latex in the laticifer network, which greatly influences the estimation of latex regeneration potential because many factors influence latex flow during the tapping process. This newly developed tapping system may cause some damage to the tree, and cannot be implemented in the latex harvest process for farmers. However, it has great advantages for studying the functions of laticifers and latex metabolism and regeneration, as the ITS simplified factors that influence latex biosynthesis and regeneration. We suggest that this system may have a great application significance for fundamental studies on laticifer biology and physiology.

CONCLUSION
In summary, non-coding RNA and transcriptome libraries were constructed in the rubber tree leaf and latex. We performed highthroughput sequencing to identify a large number of snoRNAs in the rubber tree, and 31 snoRNAs were confirmed to be highly expressed in latex. Further functional analysis of these snoRNAs will be useful for a better understanding of the complexity and dynamics of latex regeneration in the rubber tree. Thus, our findings might provide a valuable foundation for exploring the snoRNA-mediated regulatory mechanism involved in latex metabolism and regeneration.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: National Center for Biotechnology Information (NCBI) BioProject database under accession number PRJNA742874.

AUTHOR CONTRIBUTIONS
HC coordinated the project and supervised the data generation and analysis. XC and ZD conducted most of the experiments and wrote the draft of the manuscript. XZ and XH maintained the rubber tree population and collected the tapping data. DY, WW, ZA, and QL partially performed the field experiments and discussed the manuscript. HH helped to design the whole project. All authors reviewed and approved the final manuscript.