Evaluation of tolerance to Cu toxicity in the citrus rootstocks
To compare tolerance to Cu toxicity, four widely used citrus rootstocks, trifoliate orange (TO), ‘Ziyang Xiangcheng’ (XC), red tangerine (RT), and ‘Shatian’ pummelo (ST), were subjected to excessive Cu treatment, followed by evaluation of plant phenotypic and physiological parameters. 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 relative increase 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 by Cu treatment, but XC had a higher level than the other three rootstocks at 25 and 40 d of treatment. In addition, excessive Cu treatment resulted in a drastic increase of MDA in TO, RT, and ST relative to CK, but increase of MDA in XC was not conspicuous. These results indicated that XC was the most tolerant genotype to Cu toxicity among the tested rootstocks. Therefore, XC was selected for high-throughput RNAseq to reveal global transcriptome of mRNA, lncRNA, circRNA and miRNA in response to Cu toxicity.
Global responses of mRNA to Cu toxicity
From RNAseq data, we identified 30123 protein-coding genes in leaves and roots of XC by using pummelo as reference genome, 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, 3572 down-regulated) and 222 (132 up-regulated, 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). The number of DEmRNAs in the root was significantly higher than those in the leaf, suggesting that the root had more predominant responses to Cu toxicity. Moreover, 137 DEmRNAs were common between CuR/CKR and CuL/CKL, implying that they might participate in the basic response process under Cu toxicity. 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. Base on GO annotation, we found that most DEmRNAs in the leaf and root were annotated to the metabolic process, single-organism process, localization process and response to stimulus under biological process (BP), to membrane, cell, organelle, and extracellular region under cellular component (CC), and to binding, catalytic activity, transporter activity, antioxidant activity, transcription factor activity and nutrient reservoir activity under molecular function (MF) (Additional file 10: Fig. S1a and b). In particular, there were 317 DEmRNAs in response to stimulus, 596 in the membrane, 45 involved in antioxidant activity, 1 with metallochaperone activity, 32 with molecular transducer activity, 126 nucleic acid binding transcription factors, 18 with nutrient reservoir activity, 19 with receptor activity, and 188 with transporter activity in the root (Additional file 8: Table S8). GO enrichment analysis showed that the oxidation-reduction process, phosphorylation process, photosynthesis process, lignin catabolic process, phenylpropanoid catabolic process, membrane, dehydrogenase activity, peroxidase activity, protein kinase activity, iron ion binding, and heme binding were significantly enriched in the leaf and root (Fig. 2e and Additional file 10: Fig. S1c), and these GO terms of DEmRNAs possible have primarily participated in mitigation of Cu toxicity.
Global responses of lncRNA to Cu toxicity
Apart from mRNAs, we also identified 243 known lncRNAs and 1033 novel lncRNAs in XC from the RNAseq data 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 (> 4500 bp) 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, 61 down-regulated) and 5 (1 up-regulated, 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, the DElncRNAs in the Cu-treated group and CK group were clustered separately, while their three repeats were 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 terms under BP (mainly involved in metabolic process, cellular process, single-organism process, localization and response to stimulus), 13 terms under CC (mainly involved in membrane, cell, organelle and macromolecular complex), and 10 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).
Global responses of circRNA to 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, 17 down-regulated) and 17 (11 up-regulated, 6 down-regulated) DEcircRNAs were identified in CuR/CKR and CuL/CKL groups, respectively (Fig. 4f and Additional file 4: Table S4), among which, 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 highly expressed in the CuR group (Fig. 4g). These results suggest that the circRNAs are also involved in the responses to Cu toxicity.
Global responses of miRNA to Cu toxicity
In the present study, a total of 23333512, 24526156, 21822295, and 25871923 raw reads were generated from CKL, CKR, CuL and CuR small RNA 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 (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, 2 down-regulated) and 3 (all up-regulated) DEmiRNAs, and from novel miRNAs 135 (26 up-regulated, 109 down-regulated) and 127 (42 up-regulated, 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). GO enrichment analysis showed that the significantly enriched GO terms of the 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 SCC < -0.5, 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 are involved in more than 100 nodes, suggesting that they may act as 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; they included Cu-related regulators (SPL, YSL, CDPK, MAPK, and SUMO E3 Ligase, etc.), transporters (COPT, HMA, PPA, CDF, ZIP, OPT, and ABC transporter, etc.) and enzymes (LAC, CSD, CCS, PC, and COX, etc.) according to their functional description (Additional file 7: Table S7 and Additional file 8: Table S8). 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 (SCC < -0.5) 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 expression correlation between miRNAs and ceRNAs under Cu toxicity
To confirm the results of RNAseq and validate expression correlation of miRNA and their targets, six 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 for qRT-PCR. As shown in Table 1, the qRT-PCR results agreed well with RNAseq data except several low-abundance mRNAs were undetectable. In addition, the 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 negatively correlated expression between miRNAs and ceRNAs.
Expression patterns of candidate RNAs between XC and TO
To further understand the Cu tolerance mechanism of XC, comparative analysis of the expression of known Cu-related miRNAs and ceRNAs was performed between Cu-tolerant XC and Cu-sensitive TO. As shown in fig. 9, miR398b and its targets (Cg1g001620, Cg1g015400, Cg5g003300 and TCONS_0012960) were down-regulated or up-regulated in both roots of XC and TO at 1 d, 3 d and 5 d of Cu toxicity. However, the absolute values of log2FC in XC were significantly higher than those in TO at most time points. The similar changes were found for miR157c-3P and miR535c, as well as their targeted mRNAs and lncRNAs. These results suggest that the Cu-tolerant XC had occurred much more drastic expression of the Cu-related miRNAs, mRNAs and lncRNAs under Cu toxicity, which may be the important reason that leads to higher Cu tolerance in XC.