Adaptive Responses of Citrus grandis Leaves to Copper Toxicity Revealed by RNA-Seq and Physiology

Copper (Cu)-toxic effects on Citrus grandis growth and Cu uptake, as well as gene expression and physiological parameters in leaves were investigated. Using RNA-Seq, 715 upregulated and 573 downregulated genes were identified in leaves of C. grandis seedlings exposed to Cu-toxicity (LCGSEC). Cu-toxicity altered the expression of 52 genes related to cell wall metabolism, thus impairing cell wall metabolism and lowering leaf growth. Cu-toxicity downregulated the expression of photosynthetic electron transport-related genes, thus reducing CO2 assimilation. Some genes involved in thermal energy dissipation, photorespiration, reactive oxygen species scavenging and cell redox homeostasis and some antioxidants (reduced glutathione, phytochelatins, metallothioneins, l-tryptophan and total phenolics) were upregulated in LCGSEC, but they could not protect LCGSEC from oxidative damage. Several adaptive responses might occur in LCGSEC. LCGSEC displayed both enhanced capacities to maintain homeostasis of Cu via reducing Cu uptake by leaves and preventing release of vacuolar Cu into the cytoplasm, and to improve internal detoxification of Cu by accumulating Cu chelators (lignin, reduced glutathione, phytochelatins, metallothioneins, l-tryptophan and total phenolics). The capacities to maintain both energy homeostasis and Ca homeostasis might be upregulated in LCGSEC. Cu-toxicity increased abscisates (auxins) level, thus stimulating stomatal closure and lowering water loss (enhancing water use efficiency and photosynthesis).


Introduction
Copper (Cu) is required for the proper growth and development of plants, but it is extremely toxic in excess [1,2]. In order to control fruit and leaf fungal diseases, long-term and heavy application of Cu-based fungicides has caused Cu accumulation in soil of Citrus orchards. The content of available Cu in the soil increases with Citrus planting years. In the old Citrus orchards, excessive accumulation of Cu in soil is a common problem limiting Citrus production, especially in acidic soil, and is on the rise [3][4][5][6]. The physiological and molecular mechanisms of plant Cu-toxicity and Cu-tolerance have been investigated in some detail. Most studies, however, have focused on phenomena occurring at the roots, because Cu is mainly accumulated in roots exposed to Cu toxicity, and the reduction of root growth has been shown to be usually earlier than that of shoot growth [1,[7][8][9]. Less is known about how leaves deal with Cu-toxicity. Growing evidence has shown that Cutoxicity also influences biosynthesis of photosynthetic pigments, photosynthetic electron transport chain (PETC), CO 2 assimilation [2,10], production and detoxification of reactive oxygen species (ROS) [11,12], phenol metabolism [13], hormone biosynthesis [14], nitrogen (N) and carbohydrate metabolism [2,15], and cell wall formation [16]. In roots exposed to Cu-toxicity, most of the Cu is bound to the cell wall, thus preventing Cu from entering more sensitive root targets and sensitive shoots, and enhanc-

Seedling Growth and Cu Level in Roots, Stems and Leaves
'Shatian' pummelo [Citrus grandis) (L.) Osbeck] is one of the main rootstocks of pummelo. Recent work from our laboratory indicated that 'Shatian' pummelo had a higher tolerance to Cu-toxicity and was an ideal material to investigate the adaptive mechanism of Cu-toxicity [2,20]. Herein, 400 µM Cu was chosen as the Cu-toxicity treatment because it led to significant but not too severe alterations of biomass, nutrient uptake, photosynthesis and related parameters in C. grandis seedlings [2]. Additionally, we identified more differentially abundant proteins from leaves of 400 µM Cu-treated seedlings than from leaves of 200 or 300 Cu-treated seedlings [20]. As shown in Figure 1, root and shoot growth was inhibited greatly in 400 µM Cu-treated seedlings. Some fibrous roots became dark brown. Young leaf yellowing was observed in some plants. Compared to 0.5 µM Cu treatment (control), the Cu concentrations in leaves, stems and roots of 400 µM Cu-treated seedlings were increased by 494.7%, 408.0% and 929.5%, respectively. Cu concentrations were far higher in roots than in leaves and shoots of 400 µM Cu-treated seedlings. Therefore, most of Cu was accumulated preferentially in roots of 400 µM Cu-treated seedlings. mechanism underlying Cu-toxicity-induced reduction of photosynthesis at physiological and transcriptional levels.

Seedling Growth and Cu Level in Roots, Stems and Leaves
'Shatian' pummelo [Citrus grandis) (L.) Osbeck] is one of the main rootstocks of pummelo. Recent work from our laboratory indicated that 'Shatian' pummelo had a higher tolerance to Cu-toxicity and was an ideal material to investigate the adaptive mechanism of Cu-toxicity [2,20]. Herein, 400 μM Cu was chosen as the Cu-toxicity treatment because it led to significant but not too severe alterations of biomass, nutrient uptake, photosynthesis and related parameters in C. grandis seedlings [2]. Additionally, we identified more differentially abundant proteins from leaves of 400 μM Cu-treated seedlings than from leaves of 200 or 300 Cu-treated seedlings [20]. As shown in Figure 1, root and shoot growth was inhibited greatly in 400 μM Cu-treated seedlings. Some fibrous roots became dark brown. Young leaf yellowing was observed in some plants. Compared to 0.5 μM Cu treatment (control), the Cu concentrations in leaves, stems and roots of 400 μM Cu-treated seedlings were increased by 494.7%, 408.0% and 929.5%, respectively. Cu concentrations were far higher in roots than in leaves and shoots of 400 μM Cu-treated seedlings. Therefore, most of Cu was accumulated preferentially in roots of 400 μM Cu-treated seedlings.

Functional Annotation and Cu-Toxicity-Responsive Genes
Here, we identified 573 downregulated and 715 upregulated genes, including 60 upregulated and 45 downregulated transcription factors (TFs; Figure 3A -C and Tables S5 and  S6). Cluster analysis showed that the general expression profiles of DEGs were clustered separately in leaves of seedlings treated with 0.5 and 400 μM Cu, but were clustered together in three biological replicates per treatment ( Figure 3D). Cu-toxic effects on mean (± SE, n = 7 except for 5 for Chl and Car and 4 for RWC, PCs and MTs) CO 2 assimilation (A), stomatal conductance (g s , (B)), transpiration rate (Tr, (C)), ratio of intercellular to ambient CO 2 concentration (C i /C a , (D)), water use efficiency (WUE, (E)), chlorophyll (Chl) a + b (F), Chl a/b (G), carotenoids (Car, (H)), Car/Chl (I), relative water content (RWC, (J)), phytochelatins (PCs, (K)) and metallothioneins (MTs, (L)) in leaves. Different letters above the bars indicate a significant difference at p < 0.05.

Functional Annotation and Cu-Toxicity-Responsive Genes
Here, we identified 573 downregulated and 715 upregulated genes, including 60 upregulated and 45 downregulated transcription factors (TFs; Figure 3A-C and Tables S5 and S6). Cluster analysis showed that the general expression profiles of DEGs were clustered separately in leaves of seedlings treated with 0.5 and 400 µM Cu, but were clustered together in three biological replicates per treatment ( Figure 3D).  (TFs, C), and cluster analysis of DEGs (D) in leaves of C. grandis seedlings exposed to Cu-toxicity (LCGSEC).
All the assembled high-quality unigenes were first blasted against the National Centre for Biotechnology Information (NCBI) non-redundant protein sequences (NR) database using BLASTX with a cut-off E-value of 10 −5 . Majority of these genes displayed a All the assembled high-quality unigenes were first blasted against the National Centre for Biotechnology Information (NCBI) non-redundant protein sequences (NR) database using BLASTX with a cut-off E-value of 10 −5 . Majority of these genes displayed a significant sequence identity to Citrus sinensis, Citrus clementina and Citrus unshiu, which contributed 52.21%, 24.67% and 16.45% of the total assembled genes, respectively ( Figure 4A).
All unigenes and DEGs were submitted to euKaryotic Orthologous Groups (KOG) classification for functional prediction ( Figure 4B,C). There were 15,219 annotated genes (671 DEGs) assigned to 25 (23) KOG classifications. For all annotated genes, general functional prediction only (3007) contained the most genes, followed by posttranslational modification, protein turnover, chaperones (1553) and signal transduction mechanisms (1399). For DEGs, KOG classification involving the highest number of DEGs was general functional prediction only (97), followed by signal transduction mechanisms (69) and posttranslational modification (68).
A total of 882 DEGs were assigned to 199 Gene Ontology (GO) terms in cellular component, which of three (viz., photosystem (GO:0009521), photosystem I (PSI; GO:0009522), and PSII (GO:0009523) were significantly enriched with an adjusted p < 0.05. A total of 866 DEGs were mapped to 528 GO terms in molecular function, including 17 significantly enriched GO terms. Tetrapyrrole binding (GO:0046906) was the most significantly enriched GO term in molecular function, followed by pigment binding (GO:0031409) and oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, NAD(P)H as one donor, and incorporation of one atom of oxygen (GO:0016709). A total of 763 DEGs were assigned to 1573 GO terms in biological process, including 30 significantly enriched GO terms. The most significantly enriched GO term in biological process was monocarboxylic acid metabolic process (GO:0032787), followed by monocarboxylic acid biosynthetic process (GO:0072330), and photosynthesis, light harvesting in photosystem I (GO:0009768) (Table S8).

qRT-PCR Analysis
Except for Cg1g015970 and Cg2g038560, Cu-toxicity-induced expression alterations of the other 18 DEGs from RNA-Seq matched well with those from qRT-PCR. There was a significant positive linear correlation between Cu-toxicity-induced alterations of expression levels for the 20 DEGs obtained by qRT-PCR and those obtained by RNA-Seq ( Figure S1 and Table S5). Thus, the RNA-Seq results were reliable.

qRT-PCR Analysis
Except for Cg1g015970 and Cg2g038560, Cu-toxicity-induced expression alterati of the other 18 DEGs from RNA-Seq matched well with those from qRT-PCR. There a significant positive linear correlation between Cu-toxicity-induced alterations of exp sion levels for the 20 DEGs obtained by qRT-PCR and those obtained by RNA-Seq ( Fig  S1 and Table S5). Thus, the RNA-Seq results were reliable.

Hormones in Leaves
We detected 34 hormones in leaves, including two ABA and its metabolic prod ABA-GE and IAA were detected only in LCGSEC. The concentrations of the other hormones, total CKs, total GAs, total JAs, total SAs and total SLs in leaves were not significantly altered by Cu-toxicity.

Increased Immobilization of Cu in Roots, and Cu Homeostasis and Detoxification in Leaves
Our results demonstrated that Cu-toxicity increased the accumulation of Cu in C. grandis roots (Figure 1), thus limiting Cu to more sensitive shoots and enhancing C. grandis Cu-tolerance [17].
To cope with Cu-toxicity, plants have evolved a conserved and complex network of proteins to maintain Cu homeostasis, including Cu transporters, Cu chaperones and Cubinding proteins [21,40]. Here, we obtained 12 upregulated and 15 downregulated genes related to Cu homeostasis in LCGSEC (Table 1). Transporters responsible for the transport of Cu into the cytoplasm are the high affinity Cu transporter (COPT) family. COPT1 is involved in Cu acquisition and transport into leaves [41]. In Arabidopsis, tonoplast COPT5 is important for the export of Cu from the vacuole [42]. Here, we identified one downregulated COPT5 gene, three downregulated and one upregulated COPT1 genes in LCGSEC. It is worth mentioning that a total of six COPT1 genes (Cg8g023340, Cg8g023350, Cg8g023360, Cg8g023370, Cg8g023380 and Cg6g005770) were identified in C. grandis leaves (Table S3). Wan et al. [18] reported that Cu-toxicity-induced downregulation of COPT5 was greater in higher Cu-tolerant HF/Mp than in less Cu-tolerant HF/Mb ['Hanfu' (Malus domestica) scions grafted on M. prunifolia (Mp) and M. baccata (Mb), respectively] leaves, and that the expression level of COPT1 was lower in leaves of HF/Mp than in leaves of HF/Mb when exposed to Cu-toxicity. This might be an adaptive strategy to Cu-toxicity by reducing leaf Cu concentration and preventing the release of vacuolar Cu into the cytoplasm. This agrees with the increased accumulation of Cu in LCGSEC ( Figure 1). Cu is transported by the COPT/Ctr-like proteins in its reduced form Cu(I), but most of the bioavailable Cu form in soil is Cu(II). The reduction of Cu (II) to Cu (I) may facilitate the uptake of Cu in roots [40]. Ferric reduction oxidase 2 (FRO2) plays a role in Fe uptake and homeostasis [43]. Additionally, FRO2 can act as a Cu-chelate reductase and facilitate the uptake of Cu [40,44]. FRO7 is involved in Fe uptake by plastids (chloroplasts) [45]. Here, we obtained one upregulated FRO2 (Cg5g041700) gene, and four upregulated FRO7 genes (Cg1g023140, Cg2g021360, Cg5g010410 and Cg6g025130) and one downregulated (Cg8g013730) FRO7 genes in LCGSEC. The Cu-toxicity-induced upregulation of FRO2 and FRO7 might be an adaptive strategy by increasing leaf and chloroplast Fe uptake, because excess Cu reduced Fe concentration in C. grandis leaves [2]. In addition to the COPT family, zinc (Zn)-regulated transporter (ZRT)-and iron-regulated transporter (IRT)-like proteins (ZIPs) may play a role in Cu uptake [40]. Here, we obtained two downregulated ZIP genes (Zn transporter 4, chloroplastic and Zn transporter 1) in LCGSEC. Wan et al. [18] observed that Cu-toxicity downregulated the expression of ZIP2 and ZIP4 in leaves, with a greater degree in HF/Mp than in HF/Mb leaves. The transport of the nicotianamine-metal complexes across plant cell membranes is carried out by the members of the Yellow Stripe-Like (YSL) family [21]. Besides maintaining Fe homeostasis, YSL transporters are involved in distribution and redistribution of Cu [40]. Wan et al. [18] found that YSL3 was upregulated and downregulated in HF/Mb and HF/Mp leaves, respectively. Cu-toxicity-induced downregulation of YSL3 and YSL5 in C. grandis leaves might reduce the transport of Cu-nicotianamine from older to younger leaves, thus protecting younger leaves against Cu-toxicity [46].
Cu chaperones can assist Cu intracellular homeostasis by their Cu-chelating ability [47]. Here, we identified one upregulated Cu chaperone for superoxide dismutase (CCS, Cg5g009340) in LCGSEC. del Pozo et al. [48] reported similar results in roots and shoots of Arabidopsis seedlings exposed to Cu-toxicity. Wan et al. [18] reported that Cu-toxicity induced the expression of CCS in HF/Mp leaves, but not in HF/Mb leaves. Herein, we obtained three upregulated genes involved Cu ion binding (Cg2g001710, Cg5g007370 and Cg5g009340) and three upregulated Cu protein genes (Cg8g018870, Cg2g018560 and Cg3g024840) in LCGSEC, implying that these genes played a role in C. grandis Cu-tolerance by binding (cytoplasmic) free Cu ions [40]. However, Cu-toxicity induced the expression of plastocyanin (Cg3g024680) and two cytochrome c oxidase (Cg9g013180 and CgUng010240) genes in C. grandis leaves. This agrees with the report that excess Cu reduced the mRNA transcript levels of plastocyanin in Arabidopsis leaves [49]. Cu chelators play a role in internal accumulation mechanisms, in which the complexation of Cu can increase Cu immobilization in organelles such as vacuole or cell wall. The present findings and our recent work demonstrated that Cu-toxicity increased the accumulation of PCs, MTs (Figure 2), TRP ( Figure 5), GSH, total phenolics and lignin [50] in C. grandis leaves, implying that internal accumulation mechanisms played a role in Cu-tolerance of C. grandis leaves.

Cu-Toxic Effects on Cell Wall Metabolism in Leaves
As shown in Table S9, 32 downregulated and 20 upregulated genes related to cell wall metabolism were identified in LCGSEC. Wall-associated receptor kinases (WAKs) play a role in cell expansion and in defense against abiotic stress in plants. Hou et al. [51] demonstrated that WAKL4 expression was induced by excess Cu in Arabidopsis, and that an Arabidopsis WAKL4 T-DNA insertional mutant was hypersensitive to excess Cu. Xia et al. [16] observed that RNAi-mediated WAK11 knockdown lowered rice Cu-tolerance through enhancing Cu level in cytoplasm of roots and shoots and lowering Cu concentration in the cell wall (pectin and hemicellulose) of roots and shoots due to increased degree of pectin methylesterification, possibly because of decreased activity of pectin methylesterase in roots and shoots. Here, we identified five downregulated and one upregulated pectinesterase genes involved in pectin de-esterification, one downregulated omega-hydroxypalmitate O-feruloyl transferase involved in the cell wall pectin biosynthetic process, and two downregulated and two upregulated WAKs in LCGSEC, implying that pectin biosynthesis and the degree of pectin methylesterification were downregulated and upregulated in these leaves, respectively, thus decreasing and increasing Cu concentration in the cell wall and cytoplasm, respectively, and hence lowering Cu-tolerance. Xyloglucan is a major hemicellulose component in the cell wall of dicotyledonous plants [52]. Xyloglucan endotransglucosylase/hydrolases (XTHs) were shown to catalyze either the hydrolysis of xyloglucan through xyloglucan endohydrolase (XEH) activity and/or the endotransglycosylation of xyloglucan through xyloglucan endotransglucosylase (XET) activity, thus loosening the cell wall. Zhu et al. [53,54] observed that xth31 and xth17 Arabidopsis mutants had decreased xyloglucan content, slower root elongation, and less aluminum (Al) level in the root tips and cell wall, but higher Al-tolerance than wild-type plants. Here, we obtained two downregulated XTHs, indicating that Cu-toxicity might reduce xyloglucan level in LCGSEC, thus decreasing Cu accumulation. This agrees with a report showing that the Cu level was reduced in the cell wall hemicellulose and pectin in WAK11-RNAi transgenic rice roots and leaves [16]. To conclude, excess Cu might impair leaf cell wall metabolism, thus inhibiting leaf growth and lowering Cu-tolerance.

Cu-Toxic Effects on Pigment Metabolism, Photosynthesis, and Carbon, Carbohydrate and Energy Metabolisms in Leaves
Cu-toxicity-induced decreases in g s and Tr ( Figure 2) agrees with our transcriptome data that Cu-toxicity inhibited the expression of E3 ubiquitin-protein ligase RZFP34 involved in stomatal opening, but induced the expression of heat shock 70 kDa protein 1/2/6/8 involved in stomatal closure. However, Cu-toxicity-induced reduction of CO 2 assimilation was not only explained by reduced g s , because C i /C a ratio displayed an increasing trend in LCGSEC ( Figure 2). The reduction in Chl level ( Figure 2) might be due to reduced biosynthesis, as indicated by reduced expression of genes involved in Chl biosynthesis, and increased catabolization, as indicated by increased expression of genes involved in Chl catabolism, while the reduction in Car level ( Figure 2) might be due to decreased biosynthesis, as indicated by decreased expression of geranylgeranyl pyrophosphate synthase involved in Car biosynthesis. Cu-toxicity-induced inhibition of photosynthesis could not be explained by reduced photosynthetic pigment levels alone, because Cu-toxicity affected CO 2 assimilation much more than photosynthetic pigments ( Figure 2). Here, we isolated 35 downregulated and 12 upregulated genes involved in photosynthesis (ko00195; 12 downregulated and one upregulated genes), photosynthesis-antenna proteins (ko00196; eight downregulated genes), PSI (GO:0009522; 12 downregulated and one upregulated genes), PSII (GO:0009523; 14 downregulated and one upregulated genes), PETC (GO:0009767; three downregulated and one upregulated genes), PSII oxygen evolving complex (OEC, GO:0009654; four downregulated genes), photosynthesis, light reaction (GO:0019684; 12 downregulated and one upregulated genes) and thylakoid (GO:0009579; 29 downregulated and nine upregulated genes). By contrast, we only identified six upregulated genes involved in photosynthesis, dark reaction (GO:0019685) and carbon fixation in photosynthetic organisms (ko00710) (Table S10). Thus, it is reasonable to assume that Cu-toxicity-induced reduction in leaf CO 2 assimilation was mainly caused by an impaired light reaction, including (a) whole PETC (viz., photosynthesis-antenna proteins, light reaction, PSI, PSII, and PSII OEC) and (b) thylakoid [2,9,34], rather than by an impaired dark reaction. This is also supported by the report that quite a few of the genes related to PETC were downregulated in leaves of rice plants exposed to Cu-toxicity [36].
We obtained 43 upregulated and 16 downregulated genes involved in carbon, carbohydrate and energy metabolisms, including carbon metabolism (ko01200; 24 upregulated and six downregulated genes), starch and sucrose metabolism (ko00500; 15 upregulated and five downregulated genes), glycolysis/gluconeogenesis (ko00010; 16 upregulated and six downregulated genes), glycolytic process (GO:0006096; eight upregulated genes); Pyr metabolism (ko00620; 17 upregulated and two downregulated genes), citrate cycle [tricarboxylic acid (TCA) cycle] (ko00020; eight upregulated genes), ATP biosynthetic process (GO:0006754; eight upregulated and two downregulated genes), and ATP generation from ADP (GO:0006757; eight upregulated genes) in LCGSEC (Table S11). Thus, these physiological processes might be upregulated in these leaves. Cu-toxicity-induced upregulation of carbon metabolism and starch and sucrose metabolism agrees with increased accumulation of glucose, fructose, sucrose and starch in LCGSEC [2]. The increased accumulation of starch might be due to increased biosynthesis, as indicated by the upregulated expression of two genes encoding 1,4-alpha-glucan branching enzyme and one gene encoding isoamylase involved in starch biosynthesis, and reduced catabolization, as indicated by downregulated expression of two genes (α-amylase and β-amylase) correlated to starch catabolic process. The higher level of sucrose might be mainly due to elevated biosynthesis, as indicated by enhanced expression of two genes encoding sucrose synthase and one gene encoding sucrose-phosphate synthase, rather than by reduced catabolization, as indicated by elevated expression of one gene encoding β-fructofuranosidase. The elevated levels of glucose and fructose might be associated with increased formation, as indicated by enhanced expression of one β-fructofuranosidase gene involved in sucrose degradation, and two β-glucosidase and two glucan endo-1,3-β-glucosidase 1/2/3 genes involved in glucose formation.
Stressed plants often suffer from an energy deficit. Stress tolerance is highly correlated to energy availability in plants. The upregulation of glycolysis in roots of rice and roots and leaves of Citrus seedlings exposed to Al has been suggested to be an adaptive mechanism through maintaining basic respiration and meeting an elevated requirement for energy [37,38,55,56]. Mitochondria, the site of respiration and ATP synthesis, have a key role in energy metabolism. Thus, the Cu-toxicity-induced upregulation of genes involved in glycolysis/gluconeogenesis, glycolysis, Pyr metabolism, TCA cycle, ATP biosynthetic process and ATP generation from ADP might be an adaptive strategy for maintaining energy homeostasis and preventing energy shortage.

Cu-Toxic Effects on Thermal Dissipation, ROS Scavenging and Cell Redox Homeostasis in Leaves
More excess excitation energy might exist in LCGSEC due to reduced CO 2 assimilation ( Figure 2). If not rapidly dissipated, the excess excitation energy can potentially stimulate the generation of ROS. Excess excitation energy can be safely removed by xanthophyll cycle-dependent thermal dissipation before it reaches PSII reaction centers (RCs). Photorespiration can also dissipate excess light energy by consuming NADPH and ATP [57]. Cu-toxicity-induced upregulation of violaxanthin de-epoxidase involved in xanthophyll cycle and three genes involved in photorespiration (Cg4g009780, Cg6g010200 and Cg4g021720) agrees with the elevated demand for excess light energy dissipation. However, Cu-toxicity inhibited the expression of catalase (CAT) involved in scavenging H 2 O 2 produced by photorespiration (Table S12). Because Cu-toxicity stimulated the generation of ROS in C. grandis leaves [2], both the expression of some genes correlated to ROS scavenging and the levels of some antioxidants should be altered in LCGSEC. Superoxide dismutase (SOD) can rapidly dismutase superoxide anion into H 2 O 2 and O 2 . Our findings that Fe/Mn-SOD and Cu/Zn-SOD (CSD) were downregulated and upregulated in LCGSEC, respectively is in agreement with the reports that Fe-SOD abundance was reduced and CSD abundance was elevated in Cu-sufficient Arabidopsis leaves [58] and that Cu-toxicity increased CSD abundance, but had no influence on Fe-SOD abundance in C. grandis leaves [20]. Plants that suppress Fe-SOD and induce CSD under Cu-toxicity can keep superoxide anion scavenging and prevent the Cu-toxic effect on photosynthesis by buffering Cu concentration [11,58,59]. Similarly, Cu-toxicity induced the expression of CCS in leaves (Table S12), as reported in Arabidopsis [59]. CCS has been shown to bind Cu ions and deliver them specifically to CSD, thereby activating CSD [60]. Additionally, both CCS and CSD have been suggested to play a role in Cu homeostasis [18,44]. Ascorbate (ASC) peroxidase (APX) catalyzes the reduction of H 2 O 2 into H 2 O using ASC as the reducing agent. Davletova et al. [61] reported that cytosolic APX1 played a key role in protecting the chloroplast against photooxidative damage. Here, Cu-toxicity induced the expression of cytosolic APX2 (Cg6g002810) in leaves (Table S12). Besides protecting plant cells from oxidative damage by quenching reactive molecules by the addition of GSH, glutathione S-transferases (GSTs) play a role in the detoxification of HMs including Cu [62]. Lim et al. [62] found that transgenic Dianthus superbus plants overexpressing a tobacco Tau class GST (Nt107) had higher biomass, CO 2 assimilation and Cu accumulation than wild type plants when exposed to excess Cu, concluding that transgenic plants enhanced PCs biosynthesis, thereby sequestering and detoxifying excess Cu. Here, we identified one upregulated and seven downregulated GSTs in LCGSEC (Table S12). This agrees with our report that LCGSEC had increased accumulation of GSH [50]. However, Leng et al. [34] identified 27 upregulated and three downregulated GSTs in grape leaves exposed to Cu-toxicity. Peroxidases (PODs) can reduce H 2 O 2 to H 2 O using a wide variety of substrates as electron donor. Here, we obtained two upregulated and five downregulated PODs in LCGSEC. However, only seven upregulated PODs were observed in grape leaves exposed to Cu-toxicity [34]. Additionally, we obtained three downregulated and five upregulated genes related to cell redox homeostasis (GO:0045454) in LCGSEC. In conclusion, some genes correlated to thermal dissipation, photorespiration, ROS scavenging and cell redox homeostasis (Table S12) and the levels of some antioxidants such as PCs, MTs, TRP (Figures 2 and 5), GSH and total phenolics [50] were upregulated in LCGSEC, but they could not sufficiently protect these leaves from oxidative damage, because MDA level [50] and electrolyte leakage [2] were increased in LCGSEC.

Cu-Toxic Effects on Calcium Signaling and MAPK Signaling in Leaves
Calcium (Ca) signaling and MAPK signaling are the major signaling networks involved in the toxicity of HMs including Cu [24]. We obtained 29 downregulated and 11 upregulated genes involved in Ca homeostasis, including Ca ion transmembrane transporter activity (GO:0015085; five downregulated and four upregulated genes), Ca-mediated signaling (GO:0019722; one upregulated and one downregulated genes), CaM binding (GO:0005516; 12 downregulated and five upregulated genes), CaM-dependent protein kinase activity (GO:0004683; one upregulated gene), Ca ion binding (GO:0005509; 14 downregulated and three upregulated genes), and Ca channel activity (GO:0005262; two upregulated genes) in LCGSEC (Table S13). In plants, various Ca 2+ -binding proteins can act as Ca 2+ -sensors to monitor the alterations of cytosolic Ca 2+ concentration ([Ca 2+ ] cyt ), including calmodulins (CaMs), calcineurin B-like proteins (CBLs), CaM like proteins (CMLs), and Ca 2+ -dependent protein kinases (CDPKs) [24]. CaMs work through binding to and regulating the activities of diverse downstream target proteins called "CaM-binding proteins" (CaMBPs), which provide another level of specificity for Ca signaling [63]. Ca 2+ transport molecules, Ca 2+ buffers and Ca 2+ sensors are involved in the maintenance of Ca 2+ homeostasis [64]. Here, we identified three downregulated [Ca-transporting ATPase 9, plasma membrane-type (ACA9, Cg5g012700) and ACA12 (Cg2g020060 and Cg5g034830)] and one upregulated [ACA12 (Cg3g010620)] AUTOINHIBITED Ca 2+ -ATPase (ACA) genes involved in the translocation of Ca 2+ from the cytosol out of the cell or into the organelle, one upregulated Ca-transporting ATPase 4, endoplasmic reticulum (ER)-type (ECA4, Cg6g003500) involved in the translocation of Ca 2+ from the cytosol into an endomembrane compartment, two downregulated vacuolar cation/proton exchanger 3 (CAX3) genes (Cg6g006900 and Cg8g024060) involved in the translocation of Ca 2+ into the vacuole; one upregulated Ca uniporter protein 2, mitochondrial (Cg3g022310) involved in the uptake of mitochondrial Ca, and one upregulated glutamate receptor 2.8 (GLR2.8, Cg4g022490), a non-selective cation channel involved in the regulation of Ca 2+ influx into cell in LCGSEC (Table S13), indicating that Ca 2+ influx into the cytosol and Ca 2+ influx into endomembrane compartment (ER and mitochondrium) might be enhanced, but Ca 2+ efflux out of cytosol into the cell exterior and sequestration into vacuole might be downregulated in these leaves, thereby maintaining Ca 2+ homeostasis in the cytosol and the endomembrane compartment, because excess Cu reduced Ca level in C. grandis leaves [2] and that Cg3g010620 (ACA12) might be involved in the translocation of Ca 2+ from the cytosol into the organelle. Here, we obtained 14 downregulated and three upregulated genes involved in Ca 2+ binding in LCGSEC, indicating that the amount of Ca 2+ bound to Ca 2+ buffering proteins might be reduced in these leaves, thus contributing to Ca 2+ homeostasis in the cytosol. Additionally, most of DEGs [12 downregulated (Cg2g020060, Cg3g013020, Cg3g012090, Cg8g018400, Cg2g014880, Cg8g018410, Cg5g018340, Cg8g018470, Cg5g006140, Cg5g034830, Cg4g002250 and Cg5g012700) and five upregulated (Cg9g029620, Cg9g027140, Cg2g043490, Cg3g010620 and Cg8g010170) genes] involved in CaM binding were downregulated in LCGSEC.
We obtained one upregulated GLR2.8 (Cg4g022490) and one downregulated protein RALF-like 32 (RALFL32; Cg9g014450) related to Ca-mediated signaling, and one upregulated [Ca/CaM-dependent protein kinase (CaM kinase) II (CCaMKII, Cg8g010170)] and six downregulated [four CMLs (Cg1g003570, Cg2g011280, Cg3g017200 and Cg6g011100), CaM (Cg3g021480) and CBL (Cg5g018200)] genes encoding Ca sensors in LCGSEC (Table S13). A deletion of AtGLR3.5 caused premature senescence and a decrease in Chl level in Arabidopsis leaves [65]. A deletion of AtGLR3.4 led to decreased photosynthetic yield of PSII and non-photochemical quenching [66]. ZmCCaMK is considered to play a key role in ABA-and BR-induced antioxidant protection in maize leaves [67,68]. Rapid alkalinization factor (RALF) led to rapid alkalization of the cell wall by mediating a transient elevation of [Ca 2+ ] cyt , thus inhibiting cell growth in tissue culture [69]. Transgenic Arabidopsis plants overexpressing RALF22 or RALF23 displayed decreased growth and less tolerance to salt stress [70]. Thus, Cu-toxicity-induced upregulation of GLR2.8 and CCaMKII and downregulation of RALFL32 might contribute to C. grandis Cu-tolerance. Except for CCaMKII, all the other six Ca sensor genes were downregulated in LCGSEC, probably contributing to Ca 2+ homeostasis, because all the six Ca sensors were Ca 2+ binding proteins.
We identified 29 downregulated and 12 upregulated genes involved in the MAPK signaling pathway-plant (ko04016) in LCGSEC (Table S14). MAPK-signaling cascade includes three protein kinases [MAPKs, MAPK kinases (MAPKKs) and MAPKK kinases (MAPKKKs)] acting in a sequential manner to activate different downstream targets. MAPKs can be stimulated by specific metal ligands or indirectly by ROS produced due to metal stress [24]. Yeh et al. [71] reported that Cu 2+ induced MAPK activation by distinct ROS generating systems in rice roots. Additionally, ABA induced a rapid and transient MAPK activation in pea leaves [72]. Thus, Cu-toxicity might increase ABA accumulation ( Figure 5) and ROS generation [2], thus inducing MAPKK4/5 (MKK4/5) expression in C. grandis leaves. This agrees with a report that two novel rice MAPK genes (OsMSRMK3 and OsWJUMK1) in rice leaves were induced by Cu-toxicity, ABA and H 2 O 2 [73]. These results demonstrated the involvement of MAPKs in mediating Cu-toxicity. Protein phosphatases 2C (PP2Cs) are the negative regulators of stress-induced receptor kinase signaling, MAPK pathways and ABA signaling [74]. Here, we identified five upregulated genes encoding protein phosphatase 2C (PP2C) in LCGSEC (Table S14). Overexpression of an ABA, salt and drought inducible rice PP2C gene, OsPP108 conferred Arabidopsis tolerance to drought, mannitol, salt stress and ABA insensitivity [75]. Liu et al. [76] indicated that A. thaliana protein phosphatase 2C G Group 1 positively regulated salt stress in an ABA-dependent manner. Therefore, the upregulation of PP2Cs might be an adaptive response to Cu-toxicity. To conclude, our findings demonstrated the involvement of the MAPK signaling pathway in Cu-tolerance and toxicity.

Cu-Toxic Effects on Biosynthesis and Signaling of Phytohormones in Leaves
As shown in Tables S15-S17, we identified 24 upregulated and 12 downregulated genes involved in plant hormone signal transduction (ko04075), 37 upregulated and 28 downreg-ulated genes involved in auxin (AUX)-, GA-, ETH-, BR-, ABA-, SA-, CK-and JA-mediated or activated signaling pathways, 22 upregulated and 15 downregulated genes involved in ABA, CK, AUX, BR, SA, JA and GA metabolism in LCGSEC. Cu-toxicity increased the concentrations of ABA, ILA, TRP, and cZ9G, decreased the concentrations of BAP7G and 5DS, but did not affect the concentrations of the other 28 hormones in leaves ( Figure 5). These results indicated the involvement of hormone signaling in Cu-tolerance of C. grandis leaves.
Ouzounidou and Ilias [77] demonstrated that IAA alleviated Cu-toxic effects on sunflower seedlings by improving photosynthesis and WUE. Leng et al. [34] reported that all AUX response factors (ARFs) were inhibited and most genes encoding IAA synthase and AUX/IAA proteins were induced in grape leaves exposed to Cu-toxicity, suggesting that the IAA regulated genes might positively regulate grape development and Cu-tolerance. Here, we obtained four upregulated and two downregulated genes correlated to AUX signal transduction, 11 upregulated and eight downregulated genes correlated to AUX-activated signaling pathway, four upregulated and one downregulated genes correlated to AUX biosynthetic process, and seven upregulated and three downregulated genes correlated to AUX transport in LCGSEC. Additionally, Cu-toxicity increased the levels of IAA, TRP and total AUXs in leaves ( Figure 5). These results suggested that AUX biosynthesis, levels and signal signaling might be upregulated in LCGSEC, which might enhance Cu-tolerance of C. grandis leaves via promoting CO 2 assimilation and WUE.
We detected six upregulated [five PP2Cs and one ABA responsive element binding factor (ABF)] and two downregulated (two ABA receptor PYL4) and 11 upregulated and 8 downregulated genes involved in ABA signal transduction and ABA-mediated signaling pathway, respectively in LCGSEC, implying that ABA signaling might be upregulated in these leaves. However, one upregulated and four downregulated genes correlated to the ABA biosynthetic process were identified in LCGSEC. The downregulation of ABA biosynthesis-related genes caused by Cu-toxicity might be due to feedback inhibition caused by large increase of ABA level in these leaves [78]. This agrees with the reports that both the level of ABA and the expression of genes [CsPP2C5, CsABI1, Cucumis sativus SNF1related kinase 2.3 (CsSnRK2.3), CsSnRK2.4 and most CsPYLs] involved in ABA signaling were upregulated in cucumber seeds exposed to Cu-toxicity [79], and that the level of ABA was increased in leaves of sunflower seedlings exposed to Cu-toxicity [80]. ABA is well-known for its ability to stimulate stomatal closure. Cu-toxicity-induced increase in ABA level might stimulate stomatal closure, thus reducing transpiration water loss, because g s , Tr and RWC were reduced in LCGSEC ( Figure 2). ABA has been suggested to have a positive and synergistic relationship with the biosynthesis of GSH and PCs [81]. Cu-toxicity-induced increase of ABA concentration might improve the biosynthesis of GSH and PCs, as indicated by increased accumulation of GSH [50] and PCs ( Figure 5) in LCGSEC. Zehra et al. [26] reported that ABA alleviated Cu-toxicity-induced oxidative damage in A. annua leaves by reducing ROS production. The increased accumulation of ABA is also in agreement with the increased requirement for ROS scavenging [2].
One upregulated histidine-containing phosphotransfer peotein 4 (AHP4) and one downregulated two-component response regulator ARR-B family (ARR-B) involved in CK signal transduction, two upregulated genes (Cg1g013900 and Cg5g001740) involved in CK-activated signaling pathway, and one upregulated CK riboside 5 -monophosphate phosphoribohydrolase LOG7 involved in the conversion from inactive CK nucleotides to the biologically active freebase forms [82] were identified in LCGSEC (Table S17). Additionally, Cu-toxicity increased cZ9G levels, and decreased BAP7G levels, but did not significantly alter the levels of total CKs and other CKs in leaves ( Figure 5). Thus, both the levels of the active free-base form CKs and CK signaling might be upregulated in these leaves. Thomas et al. [83] reported that isopentenyltransferase (IPT)-induced CKs in transgenic tobacco enhanced Cu-tolerance and decreased Cu-toxicity-induced lipid peroxidation in leaves, which was explained by an upregulated expression of a MT-like gene (MT-L2) in leaves. Cu-toxicity might increase the levels of active CKs, thus enhancing the level of MTs and Cu-tolerance of C. grandis leaves.

Plant Culture and Cu Treatments
Plant culture and Cu treatments were carried out as described by Li et al. [2] and Huang et al. [50]. Briefly,~13-week-old 'Shatian' pummelo (Citrus grandis (L.) Osbeck) seedlings grown in 6 L pots (two plants pot −1 ) containing sand in a greenhouse with natural photoperiod at Fujian Agriculture and Forestry University, Fuzhou (26 • 5 N, 119 • 14 E) with annual average sunshine hours, temperature and relative humidity of~1600 h, 20 • C and 76%, respectively [84], were irrigated 6 times weekly for 6 months with freshly prepared nutrition solution at a Cu level of 0.5 or 400 µM from CuCl 2 until there was nutrient solution leaking out of the small hole at the bottom of the pot (~500 mL). The pH of the nutrient solutions was adjusted to 4.8 with HCl to prevent Cu precipitation. Twenty pots (a total of 40 plants) in each treatment were arranged at random. Thereafter, the most recent, fully expanded mature leaves at~7 weeks of age were selected for all analyses. After leaf gas exchange measurements, leaves without petioles, winged leaves and midribs and leaf discs of 6-mm-diameter were harvested from the same plants used for gas exchange measurements at sunny noon, frozen in liquid N 2 , and then stored at -80 • C until RNA and metabolite extraction. The unsampled seedlings were selected for the measurements of Cu in leaves, stems and roots, and RWC in leaves.

Cu Concentration in Leaves, Stems and Roots
The fully expanded mature leaves without petioles, midribs and winged leaves, the middle sections of stems and the fibrous roots were taken for subsequent analysis [85]. Leaf, stem and root Cu concentration was determined with a NexION 300X Inductively Coupled Plasma Mass Spectrometer (ICP-MS, PerkinElmer, Shelton, CT, USA) after 0.2 g of samples were digested in 5:1 (v:v) of HNO 3 :HClO 4 [2].
Leaf Chl a, Chl b and carotenoids (Car) concentrations were measured after being extracted with 80% (v:v) acetone [86].
Leaf RWC were assayed using weighing method [87]. Leaf GSH and total non-protein thiols (TNP-SH) were extracted and assayed according to Garg and Kaur [88]. Four leaf discs (0.2826 cm 2 in size) were extracted with 2 mL of ice-cold 5% (w:v) sulphosalicylic acid. For the GSH assay, 0.5 mL of supernatant was mixed with 0.6 mL of 100 mM (pH 7.0) phosphate buffer and 40 µL of 1 mM 5,5-dithiobis-2nitrobenzoic acid (DTNB). After 2 min, the absorbance was read at 412 nm. For the TNP-SH assay, 100 µL of supernatant was mixed with 0.5 mL of 0.1 M phosphate buffer (pH 7.0) containing 0.5 mM EDTA and 0.5 mL of 1 mM DTNB. After 10 min, the absorbance was read at 412 nm. The difference between TNP-SH and GSH was considered to represent PCs.
Leaf MTs was extracted and measured according to Malik et al. [89]. Four leaf discs (0.2826 cm 2 in size) were homogenized in 1.6 mL of ice-cold solution containing 0.5 M sucrose, 20 mM Tris-HCl buffer, pH 8.6 and 0.01% β-mercaptoethanol. After being centrifuged at 30,000× g for 20 min. 1 mL of supernatant was mixed with 1 mL of cold ethanol and 80 µL of chloroform; the mixtures were then centrifuged at 6000× g for 10 min at 4 ºC. The collected supernatant was mixed with 1 mg RNA and 40 µL of 37% HCl and subsequently with 3 mL of cold ethanol. The sample was maintained at -20 • C for 1 h, then centrifuged at 6000× g for 10 min. The MTs containing pellet was washed with 87% ethanol and was re-suspended in 150 µL of 0.25 M NaCl and 150 µL of 1 M HCl containing 4 mM EDTA. A volume of 4.2 mL NaCl (2 M) containing 0.43 mM DTNB buffered with 0.2 M Na-phosphate (pH 8) was added to the sample at room temperature. The sample was finally centrifuged at 3000× g for 5 min, and the supernatant was measured as absorbance at 412 nm.

Leaf RNA Extraction and RNA-Seq
Total RNA was extracted from~200 mg of frozen leaves mixed equally from four plants (one plant pot −1 ) using Recalcitrant Plant Total RNA Extraction Kit (Bioteke Corporation, Beijing, China). There were three biological replicates per treatment. A total of six sequencing libraries were constructed according to Guo et al. [37,38] and sequenced on Illumina HiSeq platform (Illumina Inc., San Diego, CA, USA) at Wuhan MetWare Biotechnology Co., Ltd. (www.metware.cn, accessed on 1 March 2021). The raw transcriptome data have been deposited in NCBI SRA database (accession number: PRJNA702620).

RNA-Seq Analysis
The raw sequencing reads were filtered using fastp v0.7.0 [90]. High-quality clean reads obtained after filtering were mapped to the genome of pummelo downloaded from the genome website (http://citrus.hzau.edu.cn/orange/download/index.php, accessed on 1 January 2021) directly using HISAT2. The mapped reads of each sample were assembled using StringTie v1.3.4b [91]. The expression level of each gene was given as fragment per kilobase of transcript per million mapped reads (FPKM). Differential expression analysis between two samples was carried out by DESeq2 software using the counts from featureCounts [92]. Then, we used Benjamini and Hochberg's method to adjust the p-valve and obtain the false discovery rate (FDR). Screening criteria for DEGs were |log 2 (fold change)| ≥ 1 and FDR < 0.05. Gene functions were annotated according to Swiss-Prot, TrEMBL, KEGG, NR, GO, KOG and Protein family (Pfam) [93].

qRT-PCR Validation
Twenty DEGs were randomly selected for qRT-PCR validation. Total RNA was extracted from leaves according to the method described above. There were 3 replicates per treatment. The sequences of the Forward and Reverse primers designed using Primier version 5.0 (Premier Biosoft International, CA, USA) were listed in Table S18. qRT-PCR was run in 2 technical replicates. U4/U6 small nuclear ribonucleoprotein PRP31 (PRPF31; Cg7g019550) and actin (Cg1g026080) were selected as internal standards [37].

Conclusions
Thirty-two upregulated and 20 downregulated genes related to cell wall metabolism were identified in LCGSEC, implying that Cu-toxicity might impair cell wall metabolism, thus reducing leaf growth and lowering Cu-tolerance. Cu-toxicity decreased and increased the expression of genes involved in Chl biosynthesis and degradation, respectively, thus reducing Chl concentration, while Cu-toxicity-induced reduction in Car concentration might be mainly due to decreased biosynthesis. The reduction in leaf CO 2 assimilation caused by Cu-toxicity might be due to impaired PETC, as indicated by the downregulated expression of related genes. Although some genes related to thermal dissipation, photorespiration, ROS scavenging and cell redox homeostasis and some antioxidants were upregulated in LCGSEC, but antioxidant systems as a whole could not protect LCGSEC from oxidative damage. In addition to reducing Cu transport from roots to shoots, several adaptive responses might occur in LCGSEC. A model for the adaptive responses of C. grandis leaves to Cu-toxicity was proposed through the integration of the present findings and available data in the previous literatures ( Figure 6). LCGSEC displayed enhanced capacities to maintain homeostasis of Cu via reducing Cu uptake by leaves and preventing vacuolar Cu into the cytoplasm, as indicated by altered expression of genes encoding Cu transporters, Cu chaperones and Cu-binding proteins and to improve internal detoxification of Cu, as indicated by increased accumulation of selected Cu chelators. LCGSEC displayed increased capacities to maintain both energy homeostasis by upregulating the expression of genes involved in energy (ATP) production and Ca homeostasis by altering the expression of related genes. Cu-toxicity increased the concentrations of ABAs (AUXs), thus stimulating stomatal closure and reducing water loss (improving WUE and photosynthesis).