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.

Introduction

Small nucleolar RNAs (snoRNAs), which range from 50 to 300 nt in size (Kiss, 2002), are a distinct class of non-coding RNAs (ncRNAs) that play important roles in ribosomal RNA (rRNA) maturation by guiding ribose 2′-O-methylation or pseudouridylation modification of the precursor rRNA (Tollervey and Kiss, 1997; Bachellerie et al., 2002; Watkins and Bohnsack, 2012). Generally, there are two subclasses of snoRNAs: C/D box snoRNAs and H/ACA box snoRNAs (Dieci et al., 2009). The C/D box snoRNAs are characterized by conserved 5′-end C box sequence (RUGAUGA) and 3′-end D box sequence (CUGA), and some contain less conserved internal copies of each box, known as C’ box and D’ box, respectively (Watkins et al., 2000; Kiss, 2002; Lui and Lowe, 2013). This subclass of snoRNAs guides methylation at the 2′-hydroxyl group (2′-O-ribose methylation) in rRNAs, transfer RNAs (tRNAs), and small nuclear RNAs (snRNAs) (Brown et al., 2003). Typically, H/ACA box snoRNAs have ACA motifs at the 3′-end and common hairpin-hinge-hairpin-tail secondary structures with H box (ANANNA) in the hinge region (Kiss-László et al., 1998; Lui and Lowe, 2013). H/ACA box snoRNAs guide pseudouridine modification in rRNAs and snRNAs, which convert uridine into pseudouridine (Ψ) (Tollervey and Kiss, 1997; Watkins and Bohnsack, 2012).

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.

Materials and Methods

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^®, 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₂, 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₂, 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^®, 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¹ 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 v1.2.8 was used according to protocols previously described (Cheng et al., 2018).

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 comparing them with rRNA sequences in H. brasiliensis with the following parameters: C/D box snoRNA > 14 and H/ACA box snoRNA > 40. If the snoRNA candidates were located in the exon, the threshold value would be set higher (C/D box snoRNA ≥ 20, H/ACA box snoRNA ≥ 45). To remove rRNAs and other structural gene sequences from the snoRNA candidates, the candidates were aligned to the eukaryotic rRNA sequences (download from ftp://ftp.ncbi.nlm.nih.gov/blast/db/), and the published rubber tree EST sequences (Ko et al., 2003; Cheng et al., 2016) using BLASTn v2.8.1 RSEM v1.2.8 (Li and Dewey, 2011) and EBSeq v1.32.0 (Leng et al., 2016) software were used to analyze the expression levels of snoRNA candidates in the leaves and latex with the false discovery rate (FDR) value set at 0.05.

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-hinge-hairpin-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.

FIGURE 1

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.

FIGURE 2

Figure 2. Flowchart of identification and prediction of small nucleolar RNAs in the rubber tree.

FIGURE 3

Figure 3. Identification of latex-abundant small nucleolar RNA (snoRNAs). (A) Heatmap of differential expressed (DE) snoRNA expression between latex and leaf by comparative transcriptomic analysis. (B) Reverse transcription polymerase chain reaction (RT-PCR) validation of 12 latex-abundant C/D subfamily snoRNAs. (C) RT-PCR validation of 19 latex-abundant H/ACA subfamily snoRNAs.

FIGURE 4

Figure 4. (A) Predicted targeting sites of C/D subfamily small nucleolar RNA (snoRNAs) on (A) 18S and (B) 28S rRNA. The nucleotides under 2′-O-methylation were indicated in red color. The C box was shown in blue, while the D box was in green. (C) Predicted targeting sites of H/ACA subfamily snoRNAs on 28S rRNA, the nucleotide under pseudouridylation was indicated in green color. H box is shown in yellow and the ACA box in purple.

FIGURE 5

Figure 5. Identification of snoRNAs positively correlated with latex yield in rubber trees. (A) Heatmap showing the qPCR expression of 26 validated snoRNAs in a Reyan66-2 × PR107 population with 50 F1 individual trees. (B) Heat map of the adjusted total solid content (TSC) (TSC/C) of F1 offspring. (C) Heat map of the Spearman correlation coefficient for each snoRNA between expression and TSC/C in F1 offspring. A high correlation coefficient was shown in red. (D–F) Plots showing the correlation coefficient of HbsnoR1, HbsnoR23, and HbsnoR46, respectively. At least three biological replicates were used for these studies.

FIGURE 6

Figure 6. Latex regeneration study identified HbsnoR28 as involved in latex regeneration regulation. (A,B) Demonstration of insulated tapping system (ITS). (C) Latex yield comparison between conventional S/2 d/3 tapping system and ITS with intervals of 6 to 72 h. No significant yield reduction was found for ITS 6 h. (D) Yield (TSC) drop trends of continuous tapping with 30-min intervals. A relatively stable TSC level was achieved after 3 h. The number denotes 14 of 50 F1 individual trees, respectively. (E) Profiles of HbsnoR28 expression and TSC/C under continuous tapping with the ITS. At least three biological replicates were used for these studies.

Results

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 according to the known sites in A. thaliana, O. sativa, and other model organisms. In total, 224 potential modification sites were identified on the rRNA sequences in the rubber tree, including 183 2′-O-methylation sites and 41 pseudouridylation sites (Table 1 and Supplementary 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.

TABLE 1

Table 1. List of putative modification sites on ribosomal RNA (rRNA) sequences in Hevea brasiliensis.

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).

TABLE 2

Table 2. Statistics of identification of latex-abundant small nucleolar RNAs (snoRNAs) in Hevea brasiliensis.

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).

TABLE 3

Table 3. List of 50 laticifer-abundant expressed small nucleolar RNA (snoRNA) candidates in Hevea brasiliensis.

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).

TABLE 4

Table 4. Targeted modification sites of 13 latex-abundant guide snoRNA on 18S and 28S rRNAs.

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² > 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² > 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² = 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).

TABLE 5

Table 5. The Spearman correlation coefficient between the expression of 13 latex-yield-related small nucleolar RNAs (snoRNAs) and latex regeneration ability.

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² = 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 high-throughput 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.

Funding

This research study was supported by grants from the National Natural Science Foundation of China (Grant Nos. 31301072 and 31870646) and Central Public-Interest Scientific Institution Basal Research Fund for Chinese Academy of Tropical Agricultural Sciences (1630022017004).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We are grateful to members of the Key Laboratory of Biology and Genetic Resources of Rubber Tree, Ministry of Agriculture and Rural Affairs, Rubber Research Institute, and Chinese Academy of Tropical Agricultural Science.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.731484/full#supplementary-material

Supplementary Figure 1 | Predicted modification sites in 5.8S, 18S, and 28S rRNA sequences in Hevea brasiliensis.

Supplementary Figure 2 | Spearman correlation coefficient between the expression of 10 laticifer-abundant small nucleolar ribonucleic acids (snoRNAs) and the total latex production and total solid content (TSC)/C (R² = 0.2–0.4).

Supplementary Table 1 | Sequences analysis of snoRNA libraries by high-throughput sequencing.

Supplementary Table 2 | Putative snoRNA candidates in rubber tree.

Supplementary Table 3 | Spearman correlation coefficient between the expression of 31 laticifer-abundant snoRNAs and the total latex production and TSC/C.

Supplementary Table 4 | Primers were used for quantitative reverse transcription polymerase chain reaction (qRT-PCR) in this study.

Footnotes

1. ^ https://gist.github.com/zachcp/c45433307a292513af41

References

Bachellerie, J.-P., Cavaillé, J., and Hüttenhofer, A. (2002). The expanding snoRNA world. Biochimie 84, 775–790. doi: 10.1016/S0300-9084(02)01402-5

CrossRef Full Text | Google Scholar

Bokszczanin, K. L., Krezdorn, N., Fragkostefanakis, S., Müller, S., Rycak, L., Chen, Y., et al. (2015). Identification of novel small ncRNAs in pollen of tomato. BMC Genomics 16:714. doi: 10.1186/s12864-015-1901-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, J. W. S., Echeverria, M., and Qu, L.-H. (2003). Plant snoRNAs: functional evolution and new modes of gene expression. Trends Plant Sci. 8, 42–49. doi: 10.1016/S1360-1385(02)00007-9

CrossRef Full Text | Google Scholar

Burroughs, A. M., Ando, Y., de Hoon, M. J. L., Tomaru, Y., Suzuki, H., Hayashizaki, Y., et al. (2011). Deep-sequencing of human Argonaute-associated small RNAs provides insight into miRNA sorting and reveals Argonaute association with RNA fragments of diverse origin. RNA Biol. 8, 158–177. doi: 10.4161/rna.8.1.14300

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C.-L., Liang, D., Zhou, H., Zhuo, M., Chen, Y.-Q., and Qu, L.-H. (2003). The high diversity of snoRNAs in plants: identification and comparative study of 120 snoRNA genes from Oryza sativa. Nucl. Acids Res. 31, 2601–2613. doi: 10.1093/nar/gkg373

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, H., Gao, J., Cai, H., Zhu, J., and Huang, H. (2016). Gain-of-function in Arabidopsis (GAINA) for identifying functional genes in Hevea brasiliensis. SpringerPlus 5:1853. doi: 10.1186/s40064-016-3523-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, H., Wang, Y., and Sun, M. (2018). “Comparison of gene expression profiles in nonmodel eukaryotic organisms with RNA-Seq,” in Transcriptome Data Analysis, eds Y. Wang and M. Sun (New York, NY: Springer New York), 3–16. doi: 10.1007/978-1-4939-7710-9_1

CrossRef Full Text | Google Scholar

Chow, J., and Heard, E. (2009). X inactivation and the complexities of silencing a sex chromosome. Curr. Opin. Cell Biol. 21, 359–366. doi: 10.1016/j.ceb.2009.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Chow, K.-S., Wan, K.-L., Isa, M. N. M., Bahari, A., Tan, S.-H., Harikrishna, K., et al. (2007). Insights into rubber biosynthesis from transcriptome analysis of Hevea brasiliensis latex. J. Exp. Bot. 58, 2429–2440. doi: 10.1093/jxb/erm093

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, W., Zhu, X., Skogerbø, G., Zhao, Y., Fu, Z., Wang, Y., et al. (2006). Organization of the Caenorhabditis elegans small non-coding transcriptome: genomic features, biogenesis, and expression. Genome Res. 16, 20–29. doi: 10.1101/gr.4139206

PubMed Abstract | CrossRef Full Text | Google Scholar

Dickenson, P. B. (1965). The Ultrastructure of the Latex Vessel of Hevea brasiliensis. Longdon: MacLaren and Sons Ltd.

Google Scholar

Dieci, G., Preti, M., and Montanini, B. (2009). Eukaryotic snoRNAs: a paradigm for gene expression flexibility. Genomics 94, 83–88. doi: 10.1016/j.ygeno.2009.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Doe, C. M., Relkovic, D., Garfield, A. S., Dalley, J. W., Theobald, D. E. H., Humby, T., et al. (2009). Loss of the imprinted snoRNA mbii-52 leads to increased 5htr2c pre-RNA editing and altered 5HT_2CR-mediated behaviour. Hum. Mol. Genet. 18, 2140–2148. doi: 10.1093/hmg/ddp137

PubMed Abstract | CrossRef Full Text | Google Scholar

Dupuis-Sandoval, F., Poirier, M., and Scott, M. S. (2015). The emerging landscape of small nucleolar RNAs in cell biology. Wiley Interdiscip. Rev. RNA 6, 381–397. doi: 10.1002/wrna.1284

PubMed Abstract | CrossRef Full Text | Google Scholar

Ender, C., Krek, A., Friedländer, M. R., Beitzinger, M., Weinmann, L., Chen, W., et al. (2008). A Human snoRNA with MicroRNA-Like functions. Mol. Cell. 32, 519–528. doi: 10.1016/j.molcel.2008.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Gébelin, V., Argout, X., Engchuan, W., Pitollat, B., Duan, C., Montoro, P., et al. (2012). Identification of novel microRNAs in Hevea brasiliensis and computational prediction of their targets. BMC Plant Biol. 12:18. doi: 10.1186/1471-2229-12-18

PubMed Abstract | CrossRef Full Text | Google Scholar

Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. doi: 10.1038/nbt.1883

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaspers, P., and Kangasjärvi, J. (2010). Reactive oxygen species in abiotic stress signaling. Physiol. Plant. 138, 405–413. doi: 10.1111/j.1399-3054.2009.01321.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawaji, H., Nakamura, M., Takahashi, Y., Sandelin, A., Katayama, S., Fukuda, S., et al. (2008). Hidden layers of human small RNAs. BMC Genomics 9:157. doi: 10.1186/1471-2164-9-157

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S. H., Spensley, M., Choi, S. K., Calixto, C. P. G., Pendle, A. F., Koroleva, O., et al. (2010). Plant U13 orthologues and orphan snoRNAs identified by RNomics of RNA from Arabidopsis nucleoli. Nucleic Acids Res. 38, 3054–3067. doi: 10.1093/nar/gkp1241

PubMed Abstract | CrossRef Full Text | Google Scholar

Kishore, S., Khanna, A., Zhang, Z., Hui, J., Balwierz, P. J., Stefan, M., et al. (2010). The snoRNA MBII-52 (SNORD 115) is processed into smaller RNAs and regulates alternative splicing. Hum. Mol. Genet. 19, 1153–1164. doi: 10.1093/hmg/ddp585

PubMed Abstract | CrossRef Full Text | Google Scholar

Kishore, S., and Stamm, S. (2006). The snoRNA HBII-52 regulates alternative splicing of the serotonin receptor 2C. Science 311, 230–232. doi: 10.1126/science.1118265

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiss, T. (2002). Small nucleolar RNAs: an abundant group of noncoding RNAs with diverse cellular functions. Cell 109, 145–148. doi: 10.1016/S0092-8674(02)00718-713

CrossRef Full Text | Google Scholar

Kiss-László, Z., Henry, Y., and Kiss, T. (1998). Sequence and structural elements of methylation guide snoRNAs essential for site-specific ribose methylation of pre-rRNA. EMBO J. 17, 797–807. doi: 10.1093/emboj/17.3.797

PubMed Abstract | CrossRef Full Text | Google Scholar

Ko, J.-H., Chow, K.-S., and Han, K.-H. (2003). Transcriptome analysis reveals novel features of the molecular events occurring in the laticifers of Hevea brasiliensis (para rubber tree). Plant Mol. Biol. 53, 479–492. doi: 10.1023/B:PLAN.0000019119.66643.5d

CrossRef Full Text | Google Scholar

Lau, N.-S., Makita, Y., Kawashima, M., Taylor, T. D., Kondo, S., Othman, A. S., et al. (2016). The rubber tree genome shows expansion of gene family associated with rubber biosynthesis. Sci. Rep. 6:28594. doi: 10.1038/srep28594

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, R. C., Feinbaum, R. L., and Ambros, V. (1993). The C. elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14. Cell 75, 843–854. doi: 10.1016/0092-8674(93)90529-y

CrossRef Full Text | Google Scholar

Leng, N., Dawson, J., and Kendziorski, C. (2016). EBSeq: An R Package for Differential Expression Analysis using RNA-seq Data. Available online at: http://www.biostat.wisc.edu/~ningleng/EBSeq_Package/EBSeq_Vignette.pdf (accessed December 19, 2016).

Google Scholar

Li, B., and Dewey, C. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12:323. doi: 10.1186/1471-2105-12-323

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, D., Qiao, H., Qiu, W., Xu, X., Liu, T., Jiang, Q., et al. (2018). Identification and functional characterization of intermediate-size non-coding RNAs in maize. BMC Genomics 19:730. doi: 10.1186/s12864-018-5103-5101

CrossRef Full Text | Google Scholar

Liang-Hu, Q., Qing, M., Hui, Z., and Yue-Qin, C. (2001). Identification of 10 novel snoRNA gene clusters from Arabidopsis thaliana. Nucl. Acids Res. 29, 1623–1630. doi: 10.1093/nar/29.7.1623

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Shi, C., Shi, C.-C., Li, W., Zhang, Q.-J., Zhang, Y., et al. (2019). The chromosome-based rubber tree genome provides new insights into spurge genome evolution and rubber biosynthesis. Mol. Plant 13, 336–350. doi: 10.1016/j.molp.2019.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, T.-T., Zhu, D., Chen, W., Deng, W., He, H., He, G., et al. (2012). A global identification and analysis of small nucleolar RNAs and possible intermediate-sized non-coding RNAs in Oryza sativa. Mol. Plant 6, 830–846. doi: 10.1093/mp/sss087

PubMed Abstract | CrossRef Full Text | Google Scholar

Lowe, T. M., and Eddy, S. R. (1999). A computational screen for methylation guide snoRNAs in yeast. Science 283, 1168–1171. doi: 10.1126/science.283.5405.1168

PubMed Abstract | CrossRef Full Text | Google Scholar

Lui, L., and Lowe, T. (2013). Small nucleolar RNAs and RNA-guided post-transcriptional modification. Essays Biochem. 54, 53–77. doi: 10.1042/bse0540053

PubMed Abstract | CrossRef Full Text | Google Scholar

McMullen, A. I. (1962). Particulate ribonucleoprotein components of Hevea brasiliensis latex. Biochem. J. 85, 491–495.

Google Scholar

Michel, C. I., Holley, C. L., Scruggs, B. S., Sidhu, R., Brookheart, R. T., Listenberger, L. L., et al. (2011). Small nucleolar RNAs U32a, U33, and U35a are critical mediators of metabolic stress. Cell Metab. 14, 33–44. doi: 10.1016/j.cmet.2011.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Nicoloso, M., Qu, L.-H., Michot, B., and Bachellerie, J.-P. (1996). Intron-encoded, antisense small nucleolar RNAs: the characterization of nine novel species points to their direct role as guides for the 2′-O-ribose methylation of rRNAs. J. Mol. Biol. 260, 178–195. doi: 10.1006/jmbi.1996.0391

PubMed Abstract | CrossRef Full Text | Google Scholar

Schubert, T., Pusch, M. C., Diermeier, S., Benes, V., Kremmer, E., Imhof, A., et al. (2012). Df31 protein and snoRNAs maintain accessible higher-order structures of chromatin. Mol. Cell. 48, 434–444. doi: 10.1016/j.molcel.2012.08.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Shearman, J. R., Sangsrakru, D., Ruang-areerate, P., Sonthirod, C., Uthaipaisanwong, P., Yoocha, T., et al. (2014). Assembly and analysis of a male sterile rubber tree mitochondrial genome reveals DNA rearrangement events and a novel transcript. BMC Plant Biol. 14:45. doi: 10.1186/1471-2229-14-45

PubMed Abstract | CrossRef Full Text | Google Scholar

Taft, R. J., Glazov, E. A., Lassmann, T., Hayashizaki, Y., Carninci, P., and Mattick, J. S. (2009). Small RNAs derived from snoRNAs. RNA 15, 1233–1240. doi: 10.1261/rna.1528909

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, C., Qi, J., Li, H., Zhang, C., and Wang, Y. (2007). A convenient and efficient protocol for isolating high-quality RNA from latex of Hevea brasiliensis (para rubber tree). J. Biochem. Biophys. Methods 70, 749–754. doi: 10.1016/j.jbbm.2007.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, C., Yang, M., Fang, Y., Luo, Y., Gao, S., Xiao, X., et al. (2016). The rubber tree genome reveals new insights into rubber production and species adaptation. Nat. Plants 2:16073. doi: 10.1038/nplants.2016.73

PubMed Abstract | CrossRef Full Text | Google Scholar

Tollervey, D., and Kiss, T. (1997). Function and synthesis of small nucleolar RNAs. Curr. Opin. Cell Biol. 9, 337–342. doi: 10.1016/S0955-0674(97)80005-80001

CrossRef Full Text | Google Scholar

Tupý, J. (1988). Ribosomal and polyadenylated RNA content of rubber tree latex, association with sucrose level and latex pH. Plant Sci. 55, 137–144. doi: 10.1016/0168-9452(88)90169-90160

CrossRef Full Text | Google Scholar

Watkins, N. J., and Bohnsack, M. T. (2012). The box C/D and H/ACA snoRNPs: key players in the modification, processing and the dynamic folding of ribosomal RNA. Wiley Interdiscip. Rev. RNA 3, 397–414. doi: 10.1002/wrna.117

PubMed Abstract | CrossRef Full Text | Google Scholar

Watkins, N. J., Ségault, V., Charpentier, B., Nottrott, S., Fabrizio, P., Bachi, A., et al. (2000). A common core rnp structure shared between the small nucleoar box C/D RNPs and the spliceosomal U4 snRNP. Cell 103, 457–466. doi: 10.1016/S0092-8674(00)00137-139

CrossRef Full Text | Google Scholar

Yamashita, S., Yamaguchi, H., Waki, T., Aoki, Y., Mizuno, M., Yanbe, F., et al. (2016). Identification and reconstitution of the rubber biosynthetic machinery on rubber particles from Hevea brasiliensis. eLife 5:e19022. doi: 10.7554/eLife.19022

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J.-H., Zhang, X.-C., Huang, Z.-P., Zhou, H., Huang, M.-B., Zhang, S., et al. (2006). snoSeeker: an advanced computational package for screening of guide and orphan snoRNA genes in the human genome. Nucleic Acids Res. 34, 5112–5123. doi: 10.1093/nar/gkl672

PubMed Abstract | CrossRef Full Text | Google Scholar