Eight-year-old H. brasiliensis cultivar CATAS7-33-97 SRJCs and CATAS7-33-97 DCs were cultivated in the experimental field of Guangdong NongKen Tropical Crop Science Institute. Three SRJC groups (10 trees in each group) and three DC groups (10 trees in each group). Latex (10 ml from each plant) was collected from SRJCs or DCs. The latex was allowed to drop directly into liquid nitrogen in an ice kettle. After the latex drops in liquid nitrogen, small spherical particles are formed immediately. About 2–3 min, 10 ml of latex can be collected. Then the frozen latex was transferred to an RNAase-free tube and stored in liquid nitrogen. After all, samples were collected, they were returned to the laboratory and stored at −70°C or used immediately. Before RNA extraction, the samples in each group were mixed in equal quantities, and then RNA extraction buffer was added to start RNA extraction.

The latex total RNA was isolated by the described method (Kush et al., 1990; Tang et al., 2007) with some modification. In brief, the latex powder was suspended in 10 ml RNA extraction buffer containing 1,000 mol.L⁻¹ Tris–HCl (pH 8.0), 100 mmol.L⁻¹ EDTA, 1.4 mol.L⁻¹ LiCl, 20% (v/v) sodium dodecyl sulfate and 5% (v/v) mercaptoethanol, and the suspension was incubated at room temperate for 20 min, and was centrifuged at 12,000 × g for 10 min at 4°C. The lower liquid was transferred to a fresh tube. Then an equal volume of phenol/chloroform/isoamyl alcohol (25:24:1, v/v) was added and mixed evenly. The mixture was centrifuged at 12,000 × g for 10 min at 4°C. The supernatant was transferred to a fresh tube, an equal volume of phenol/chloroform/isoamyl alcohol (25:24:1, v/v) was added again and mixed evenly. The mixture was centrifuged at 112,000 × g for 15 min at 4°C. The supernatant was transferred to a fresh tube, and 1/3 vol of 8 mol.L⁻¹ LiCl was added, and mixed well, and stored at 4°C overnight. The reaction mixture was centrifuged at 12,000 × g for 15 min at 4°C. The final pellet was washed with 70% (v/v) ethanol, air-dried, and dissolved in RNase-free water. Total RNA was monitored using agarose gel electrophoresis. The concentration and integrity of total latex RNA were measured using Qubit® RNA Assay Kit in Qubit® 2.0 Fluorometer (Life Technologies, CA, United States) and the RNA Nano 6,000 Assay Kit of the Bioanalyzer 2,100 system (Agilent Technologies, CA, United States).

For lncRNA sequencing, the ribosomal RNA was depleted from the latex total RNA according to the Epicenter Ribo-Zero Gold Kit (Illumina, San Diego, United States). Then, six ssRNA-seq libraries (3 replicates in each sample) were constructed by the protocol of the Illumina TruSeq™ RNA sample prep Kit (Illumina, CA, United States) and were sequenced on the Illumina Hiseq 2,500 platform (LC-Bio, China). For sRNA sequencing, six sequencing libraries (3 replicates in each sample) were constructed following the protocol of the NEB Next Ultra-small RNA Sample Library Prep Kit for Illumina (NEB, Ipswich, MA, United States), and sequencing was carried out on the Illumina Hiseq 2,500 platform. Library construction and sequencing were carried out by Novogene Co., LTD (Beijing, China).

To generate clean data from ssRNA-seq, raw reads of FastQ format were firstly processed through in-house Perl scripts using the quality control software FastQC.¹ In this step, clean data were obtained by removing reads containing adapter, reads containing more than 10% ploy-N, reads with average quality scores of each read of less than 20, and reads containing more than 40% of bases with quality scores less than 20. All obtained clean reads were mapped against the rubber tree genome (GenBank under the accession LVXX01000000, Tang et al., 2016) with STAR (v2.5.1b; Dobin et al., 2013) using the method of Maximal Mappable Prefix and assembled into transcripts by using the Cufflinks software (Trapnell et al., 2012). HTSeq v0.6.0 was used to count the reads numbers mapped to each gene (Anders et al., 2015). The raw count data were subsequently subjected to library-size normalization by edgeR software (Robinson et al., 2010). The expression level of transcripts was calculated using the Fragments Per Kilobase per Million mapped reads (FPKM).

According to a previously described method (Ding et al., 2019), lncRNAs were identified according to the following criteria: (i) transcripts that overlapped with known protein-coding genes on the same strand, transcripts with FPKM scores less than 0.5, and transcripts shorter than 200 nt were filtered; (ii) the transcripts with protein-coding potential were also discarded according to the evaluation of Coding-Potential Assessment Tool (Wang et al., 2013a), Coding Potential Calculator (Kong et al., 2007), and Coding-Non-Coding Index (Sun et al., 2013); and (iii) the transcripts with well-known protein domains were also removed based on the Pfam-hidden Markov models (Finn et al., 2008). The remaining transcripts were regarded as reliably expressed lncRNAs. DEGs, DELs between SRJCs and DCs were identified by DESeq2 by setting |log₂fold-change| > 1 and p < 0.05 (Love et al., 2014).

For sRNA analysis, clean reads were obtained by removing contaminated reads, low-quality reads, and short sequences with base lengths less than 18 nt from raw data. Then, the clean reads were aligned to the GenBank (Release 209.0) and Rfam (Release 11.0) databases, and the reads annotated as small nucleolar RNA, small nuclear RNA, small cytoplasmic RNA, ribosomal RNA, and translocation RNA were removed. The remaining clean reads were mapped to the rubber tree reference genome (Tang et al., 2016) by using the Bowtie software (Langmead et al., 2009). The known miRNAs were identified by searching against miRbase (Release 22; Griffiths-Jones et al., 2006) and PNRD (Yi et al., 2015) databases. Novel miRNAs were identified according to their genome positions and hairpin structures predicted by Mireap v0.2 software (Wen et al., 2010; Friedlander et al., 2011). The microRNAs expression levels were estimated by TPM (transcript per million) through the following criteria (Zhou et al., 2010): Normalization formula: Normalized expression = mapped read count/Total reads ×10⁶. Differential expression miRNA (DEMs) between SRJCs and DCs was determined using the DESeq R package (1.8.3). The p-values were adjusted using the Benjamini and Hochberg method (Benjamini and Hochberg, 1995). Corrected p-value of 0.05 was set as the threshold for significantly differential expression by default.

Identifying the targets of lncRNAs in cis- or trans regulation was implemented according to a previously described method (Wu et al., 2019; Fu et al., 2020). In brief, DELs and their nearby protein-encoding genes, which were located 100 kb downstream and upstream of DELs, were used to conduct co-expression analysis and identify the targets of lncRNAs in cis-regulation. DEL–DEG pairs that were nearly located and highly co-expressed were considered as cis-acting regulation. The target genes of trans-acting lncRNAs were identified by the correlation of expression levels. Pearson correlation coefficients (Pripp, 2018) between the expression of DELs and DEGs were calculated, and DEGs with an absolute value of the correlation coefficient of more than 0.95 were predicted as trans-target genes of DELs.

DELs as precursors of DEMs were identified by aligning DELs to miRbase (Release 22; Griffiths-Jones et al., 2006) by Blast. Those with identity higher than 90% were regarded as putative miRNA precursors. In addition, the software miRPara (Wu et al., 2011) was also used to predict DEM precursors from DELs.

To identify DELs that were targeted by DEMs, the DELs and DEMs data were predicted by psRNATarget (Dai et al., 2018) with expectation ≤3 based on perfect or near-perfect complementarity between the miRNAs and the target lncRNAs sequences. DELs acting as eTMs of DEMs were predicted by TAPIR (Bonnet et al., 2010) according to the following criteria (Wu et al., 2019): (i) perfect base pairing was required at the 2nd to 8th positions from the 5′ end of miRNA sequence; (ii) the total number of mismatches and G/U pairs within the miRNA and DEL pairing region (excluding the bulge) should be no more than three; (iii) a three-nucleotide bulge was only permitted at the 9th to 12th positions from the 5′ end of miRNA sequence.

Based on the above prediction results, the regulatory networks of DEL-DEG, DEM-DEG, and DEM-DEL were constructed visualized with Cytoscape (v3.8.0; Su et al., 2014). Subsequently, the interaction networks of DELs, DEMs, and DEGs were imported into Cytoscape for visualization (Su et al., 2014).

The enrichment of the Gene Ontology (GO) terms² and Kyoto Encyclopedia of Genes Genomes (KEGG) pathways³ in the target genes of DELs were performed by the cluster Profiler R package (Yu et al., 2012) setting p value <0.05 as significantly enriched.

cDNA was synthesized following the protocol of cDNA Synthesis Kit (TaKaRa, China). qRT-PCR was performed using SYBR Premix Ex Taq Kit (TaKaRa, China) on the Mx3005P Real-Time PCR System (Stratagene, United States). The HbACT7 was used as the internal reference gene (Li et al., 2016a). The primers are listed in Supplementary Table S1. The expression of lncRNAs and mRNAs were normalized using a 2^−∆∆CT method (Schmittgen and Livak, 2008).

The systematic identification of lncRNAs in rubber trees was performed in latex from SRJCs and DCs using ssRNA-seq. In total, 642 million clean reads were obtained by all ssRNA-seq data (Supplementary Table S2). A total of 48031 transcripts were obtained. Of these transcripts, 44,374 overlapping with 25,738 rubber tree protein-encoding genes were removed. A total of 3,657 lncRNAs were identified from the rubber tree, and 1,128 of these lncRNAs were novel lncRNAs (Supplementary Table S3). Approximately 50.4, 30.9, and 8.9% of lncRNAs had lengths of 200–600, 601–1,000, and 1,001–1,400 bp, respectively, and the remaining 10.8% of lncRNAs had a length of more than 1,400 bp (Figure 1A). Most (84.3%) lncRNAs contained two or three exons, and approximately 8.4, 2.8, and 1.7% of lncRNAs contained four, five, and six exons, respectively. The remaining 2.8% of lncRNAs contained more than six exons (Figure 1B). According to their genomic locations, 40.6% of lncRNAs were long intervening noncoding RNAs (lincRNAs), approximately 29.6 and 29.7% of lncRNAs were antisense and sense lncRNAs, and 0.1% of lncRNAs was overlapping lncRNAs, respectively (Figure 1C).

In total, 3,013 and 2,982 lncRNAs were identified in SRJCs and DCs, respectively. A total of 196 DELs, including 57 downregulated and 139 upregulated lncRNAs, were identified between SRJCs and DCs (Figure 2A). Among DELs, 38 lncRNAs were unique in DCS, and 99 lncRNAs were unique in SRJCs (Figure 2B).

GO enrichment and KEGG pathway analyses were performed on the 1796 target genes of DELs. Under three GO categories, the term with the largest proportion in biological process was cellular nitrogen compound metabolic process (Figure 3). The analysis results of KEGG pathway enrichment showed that the target genes of DELs were associated with photosynthesis, plant hormone signal transduction, glutathione metabolism, oxidative phosphorylation, protein export, oxidative phosphorylation, and amino acids biosynthesis (Figure 4).

A total of 122 lncRNA–mRNA pairs related to cis-acting regulation were obtained (Supplementary Table S4). Among these lncRNA–mRNA pairs, XR_002492989.1 was spaced 4.26 kb upstream of XM_021798588.1, which encoded a 3-hydroxy-3-methylglutaryl-coenzyme A reductase (HMGR) that catalyzes the conversion of HMG-CoA to mevalonic acid, which is a central step for the biosynthesis of a number of natural products (Figure 5A). XR_002492331.1 was located 6.69 kb upstream of XM_021793515.1, which encoded a polyphenol oxidase (Figure 5B), and XR_002492403.1 was located 6.67 kb upstream of XM_021794198.1, which encoded a peroxidase 64-like (Figure 5C), suggesting that XR_002492331.1 and XR_002492403.1 were related to reactive oxygen species metabolism. XR_002490996.1 was located 9.93 kb upstream of XM_021783538.1, which encoded a MADS-box transcription protein (Figure 5D), and TCONS_00023828 was located 0.45 kb upstream of XM_021814952.1, which encoded a WRKY transcription factor (Figure 5E). TCONS_00003846 was spaced 4.19 kb upstream of XM_021783906.1, which encoded a DNA (cytosine-5)-methyltransferase CMT3 that participated in the DNA methylation (Figure 5F). These results suggested that these cis-acting DELs might involve in the regulation of their adjacent genes connected with natural rubber biosynthesis, reactive oxygen species metabolism, transcription regulation, and DNA methylation.

A total of 8 DELs and 4 genes were randomly selected to examine the expression by qRT-PCR. qPCR results showed the expression of selected DELs and genes were highly consistent with ssRNA-seq data (Figure 6).

A total of 15,778,815 and 11,396,302 raw reads were obtained in SRJCs and DCs, respectively. After removing reads less than 18 nt, adapter sequences, and low-complexity sequences, averages of 3,354,892 and 3,867,160 unique reads were obtained (Supplementary Table S5). More than 80% of the sRNAs had lengths of 18–24 nt.

In total, 136 miRNAs, including 108 novel and 28 known miRNAs, were found in DCs and SRJCs. About 41.96, 29.23, and 14.61% of these miRNAs had lengths of 21, 24, and 22 nt, respectively (Figure 7A). Among 136 miRNAs, 133 were common to DCs and SRJCs, and two and one miRNAs were specific to DCs and SRJCs, respectively (Figure 7B). Eleven differential expression miRNAs, including four upregulated and seven downregulated miRNAs, were found between DCs and SRJCs (Figure 7C).

The construction of the lncRNA-miRNA-mRNA network is helpful to understand the regulatory role of DELs in rubber trees. We constructed the lncRNA–miRNA–mRNA network based on the co-expression of mRNA and DELs, and DELs–miRNA interaction pairs. As shown in Figure 8, the lncRNA–miRNA–mRNA network contained 5 miRNAs (i.e., hbr-miR2118, hbr-miR47, hbr-miR6169, hbr-miR6171, and hbr-miR6485), 7 lncRNAs (e.g., TCONS_00093868, XR_002494646.1, and XR_002493724.1), and 23 mRNAs (e.g., XM_021806494.1 encoding small rubber particle protein (SRPP), XM_021836414.1 encoding superoxide dismutase, XM_021788744.1 encoding ethylene-responsive transcription factor ERF091, XM_021788744.1 and XM_021836818.1 encoding DNA (cytosine-5)-methyltransferase DRM2, and XM_021803144.1 encoding histone deacetylase 9 (Supplementary Table S7).