Tolerance evaluation of citrus rootstocks to Cu toxicity
To evaluate the tolerance of commonly used citrus rootstocks to Cu toxicity, trifoliate orange (TO), ‘Ziyang Xiangcheng’ (XC), red tangerine (RT), and ‘Shatian’ pummelo (ST) were subjected to excess Cu treatment. As shown in Fig. 1, after 25 d of treatment, top leaves of TO and ST showed a yellow color, while those of XC and RT were almost normal (similar to the phenotype of CK). Although the growth rate of plant height was significantly suppressed upon Cu toxicity, XC had a minimal impact among the four rootstocks. Chlorophyll contents were also significantly reduced upon Cu toxicity, but that in XC was obviously higher than that in TO, RT, and ST at 25 d and 40 d of treatment. In addition, excess Cu treatment resulted in a significant increase of MDA in TO, RT, and ST compared with CK, but increase of MDA in XC was inconspicuous. These results indicated that XC was the most tolerant to Cu toxicity among the tested rootstocks, followed by RT, while TO and ST were sensitive. Therefore, XC was selected to perform high-throughput RNA sequencing (RNAseq) for uncovering the underlying tolerant mechanisms.
mRNA expression profiles under Cu toxicity
Using pummelo as reference genome, we identified 30123 genes in leaves and roots of XC, and their log2FC values are presented as Volcano Plot pictures in Fig. 2a and 2b, which ranged from -8.2 to 7.9 in roots and from -5.3 to 4.8 in leaves. Among all of these genes, 5734 (2162 up-regulated and 3572 down-regulated) and 222 (132 up-regulated and 90 down-regulated) DEmRNAs were identified in Cu-treated root (CuR/CKR) and leaf (CuL/CKL) groups, respectively (Fig. 2c and Additional file 1: Table S1). Moreover, 137 DEmRNAs were common between CuR/CKR and CuL/CKL. A heat map of DEmRNAs showed the general expression profiles of DEmRNAs in each treatment and also showed that the three repeats of each treatment always clustered together while the Cu-treated group and the CK group were clustered separately (Fig. 2d).
To explore the functions of the DEmRNAs, GO annotation and GO enrichment analysis were performed. As shown in Additional file 10: Fig. S1a and b, 12 GO terms under biological process (BP), 10 GO terms under cellular component (CC), and 9 GO terms under molecular function (MF) were annotated for DEmRNAs of the leaf, while 16 GO terms under BP, 11 GO terms under CC, and 13 GO terms under MF were annotated for DEmRNAs of the root. In both the leaf and root, most DEmRNAs were annotated to the metabolic process, single-organism process, localization process and response to stimulus under BP, to membrane, cell, organelle, and extracellular region under CC, and to binding, catalytic activity, transporter activity, antioxidant activity, transcription factor activity and nutrient reservoir activity under MF. In addition, signaling process, detoxification process, molecular transducer activity, and metallochaerone activity were annotated in the root specifically. GO enrichment analysis showed that significantly enriched GO terms in the leaf included the lignin catabolic process, phenylpropanoid catabolic process, apoplast region, golgi subcompartment, extracellular region, oxygen oxidoreductase activity, and peptidase regulator activity, while those in the root included the photosynthesis process, lignin catabolic process, microtubule-based movement process, oxidation-reduction process, phosphorylation process, MCM complex, photosystem I reaction center, extracellular region, membrane, dehydrogenase activity, peroxidase activity, protein kinase activity, iron ion binding, and heme binding (Fig. 2e and Additional file 10: Fig. S1c).
LncRNA expression profiles under Cu toxicity
In total, we identified 243 known lncRNAs and 1033 novel lncRNAs in the root and leaf of XC by blasting to known lncRNAs of citrus in the GREENC database and performing CNCI, CPC, CPAT and PfamScan analysis (Fig. 3a). Comparison of the genomic characterizations of the lncRNAs with mRNAs showed that their transcripts were similar in length distribution, except lncRNA had relative higher numbers of long transcripts (> 4500bp) than mRNA; for exon number, a higher percentage of lncRNAs had 2 to 4 exons; in addition, lncRNAs had a shorter ORF length and lower FPKM value (Fig. 3b–e). At a cutoff with an absolute value of log2FC > 1 and FDR < 0.05, 164 (103 up-regulated and 61 down-regulated) and 5 (1 up-regulated and 4 down-regulated) DElncRNAs were identified in CuR/CKR and CuL/CKL groups, respectively (Fig. 3f and Additional file 2: Table S2). The log2FC values of DElncRNAs ranged from -10.0 to 13.2 in the root, and -11.1 to 8.8 in the leaf. The general expression profiles of DElncRNAs are shown in Fig. 3g. Similar to DEmRNAs and DElncRNAs in the Cu-treated group and CK group were clustered separately, while their three repeats always clustered together.
To explore the potential functions of these DElncRNAs, their cis- and trans-targeted mRNAs were predicted with bioinformatics methods (Additional file 3: Table S3). GO annotation of the targets of DElncRNAs in the root showed that they were annotated in 16 GO terms under BP (mainly involved in metabolic process, cellular process, single-organism process, localization and response to stimulus), 13 GO terms under CC (mainly involved in membrane, cell, organelle and macromolecular complex), and 10 GO terms under MF (mainly involved in binding, catalytic activity, transporter activity, electron carrier activity, transcription factor activity, and antioxidant activity) (Additional file 11: Fig. S2). GO enrichment analysis of these targets showed that significantly enriched GO terms were the photosynthesis process, phosphorylation process, oxidation-reduction process, lignin catabolic process, phenylpropanoid catabolic process, MCM complex, photosystem I reaction center, extracellular region, membrane, dehydrogenase activity, peroxidase activity, protein kinase activity, iron ion binding, and heme binding (Fig. 3h).
CircRNA expression profiles under Cu toxicity
In total, 2404 circRNAs were identified in the leaf and root of XC, and 60.48%, 28.62%, and 10.90% of them belong to the intergenic region type, exon type, and intron type, respectively (Fig. 4c). The sequence length distribution of circRNAs is shown in Fig. 4a, and most of them were 10000 bp to 50000 bp, or shorter than 1200 bp. Chromosome 2 (chr2) included maximum circRNAs, followed by chr3, chr5, and chr8 (Fig. 4b). Log2FC values of circRNAs are displayed in Fig. 4d and 4e, and 45 (28 up-regulated and 17 down-regulated) and 17 (11 up-regulated and 6 down-regulated) DEcircRNAs were identified in CuR/CKR and CuL/CKL groups, respectively (Fig. 4f and Additional file 4: Table S4), among which, only 1 DEcircRNAs were common between CuR/CKR and CuL/CKL. A heat map showed the general expression profiles of DEcircRNAs in each treatment, and most DEcircRNAs were more highly expressed in the CuR group (Fig. 4g).
miRNA expression profiles under Cu toxicity
In the present study, a total of 23333512, 24526156, 21822295, and 25871923 raw reads were generated from CKL, CKR, CuL and CuR libraries, respectively. Of these raw reads, we obtained 15050388, 15173260, 13474928 and 13748074 clean reads after removing adaptor sequences, low-quality sequences, and reads shorter than 18 nt and longer than 32 nt. The lengths of most clean reads were 20–24 nt; most were 21 nt, followed by 24 nt (Fig. 5a). Small RNA classification showed that 81% of clean reads were rRNA (42%) and unmatched (39%), and there were also 14% miRNA, 5% tRNA, and 1% of other types (Fig. 5b). From 14% of clean reads, we finally identified 149 known miRNAs and 336 novel miRNAs. The top 10 expressed miRNAs in each sample are shown in Fig. 5c and 5d, and miR166a-3p and nov-m0105-3p exhibited the highest expression abundance. From known miRNAs 12 (10 up-regulated and 2 down-regulated) and 3 (all up-regulated) DEmiRNAs, and from novel miRNAs 135 (26 up-regulated and 109 down-regulated) and 127 (42 up-regulated and 85 down-regulated) DEmiRNAs were identified in CuR/CKR and CuL/CKL groups, respectively (Fig. 5e and 5f and Additional file 5: Table S5). The general expression profiles of these DEmiRNAs are shown in Fig. 5g and 5h. Their expression levels exhibited obvious differences between CK and Cu-treated samples and between root and leaf samples.
Targeted mRNAs of these DEmiRNAs are listed in Additional file 6: Table S6. We found that 84.7% of DEmRNAs in the leaf (188/222) and 81.0% of DEmRNAs in the root (4642/5734) were targeted by one or multiple DEmiRNAs. GO annotation of the targets of DEmiRNAs in the root showed that most of them were annotated to the metabolic process, cellular process, single-organism process, localization process, and response to stimulus under BP, to membrane, cell, organelle and macromolecular complex under CC, and to binding, catalytic activity and transporter activity under MF (Additional file 12: Fig. S3a and b). In addition, most targets of known DEmiRNAs were down-regulated. GO enrichment analysis showed that significantly enriched GO terms of targets of known DEmiRNAs in root were photosynthesis, microtubule-based movement, carbohydrate metabolic process, DNA polymerase complex, microtubule, membrane, cellulose synthase activity, microtubule binding, and catalytic activity, while those of novel DEmiRNAs were microtubule-based movement process, phosphorylation process, protein modification process, oxidation-reduction process, DNA polymerase complex, kinesin complex, extracellular region, membrane, protein kinase activity, transferase activity, anion binding, and catalytic activity (Fig. 6a and 6b).
CeRNA regulatory network in response to Cu toxicity
To reveal the global regulatory network of protein-coding RNAs and ncRNAs under Cu toxicity, a ceRNA network was constructed using DEmiRNAs, DEmRNAs, DElncRNAs, and DEcircRNAs based on ceRNA theory. In total, 5739 DEmRNAs, 64 DElncRNAs, and 5 DEcircRNAs were predicted as targets of 251 miRNAs in the root and leaf. When their correlation was further filtered with PCC > 0.8, we obtained 3819 DEmiRNA-DEmRNA and 10 DEmiRNA-DElncRNA interactions in the root and 12 DEmiRNA-DEmRNA interactions in the leaf (Fig. 7). In this ceRNA network, Nov-m0238-3p, Nov-m0074-5p, Nov-m0183-3p, miR166c-5p, Nov-m0128-3p, Nov-m0328-5p, miR165a-5p, and miR535c involved in more than 100 nodes and were considered the core regulators. In addition, lncRNAs including TCONS_00012501, TCONS_00012960, TCONS_00025983, TCONS_00027274, TCONS_00034874, TCONS_00036810, TCONS_00042747, TCONS_00051884 and TCONS_00062474 participated in the network.
Considering that the network contains enormous information and each one cannot be displayed in the figure, we constructed a mini-ceRNA network by reducing the mRNA amount. We identified 284 known Cu-related mRNAs that were reported directly or indirectly in previous literatures from all DEmRNAs of the root (Additional file 7: Table S7). The mini-ceRNA network was constructed with these 284 DEmRNAs and all DEmiRNAs, DElncRNAs, and DEcircRNAs. Finally, only 261 DEmiRNA-DEmRNA and 10 DEmiRNA-DElncRNA interactions (PCC > 0.8) containing 18 DEmiRNAs, 129 DEmRNAs, and 9 DElncRNAs were obtained and are displayed in Fig. 8. DEmiRNAs including miR166a-5p, miR395c, miR535c, miR395k, miR166c-5p, miR165a-5p, and miR399a were involved in more than 20 nodes, and all of them were up-regulated in the CuR. A known Cu-related key miRNA, miR398b, was identified in the network, which down-regulated and interacted with 5 DEmRNAs and 2 DElncRNAs. In addition, many known Cu-related key mRNAs such as SPL (Cg5g011720, Cg6g012520, Cg7g016770), YSL (Cg7g013630), HMA (Cg5g002920, Cg5g002930, Cg4g021370), ABC transporter (Cg5g018290, Cg5g027620, Cg3g011050, Cg3g009290, and Cg5g021160 etc.), LAC (Cg6g004840, Cg7g002640, Cg6g004880) and ZIP (Cg8g019240, Cg9g029160, Cg2g037280) were down-regulated in the network.
qRT-PCR validated mRNA-miRNA-LncRNA interactions under Cu toxicity
To confirm the results of RNAseq and validate predicted interactions of miRNA and their targets preliminarily, 6 miRNAs (miR398b, miR8175, miR157c-3p, miR166a-5p, miR165a-5p, and Nov-m0284-5p) and their targeted mRNAs and lncRNAs were selected from the ceRNA network to determine their expression levels by qRT-PCR. As shown in Table 1, qRT-PCR results agreed well with RNAseq data except several low-abundance mRNAs were undetectable. In addition, miRNA and its targets showed quite correct up- or down-regulated relationships. For example, miR398 was down-regulated and all of its predicted targets were up-regulated; miR8175 was up-regulated and all of its targets were down-regulated. This result not only suggests reliability of RNAseq data in this study, but also validates the constructed ceRNA network in a certain way.