Magnesium Deficiency Induced Global Transcriptome Change in Citrus sinensis Leaves Revealed by RNA-Seq

Magnesium (Mg) deficiency is one of the major constraining factors that limit the yield and quality of agricultural products. Uniform seedlings of the Citrus sinensis were irrigated with Mg deficient (0 mM MgSO4) and Mg sufficient (1 mM MgSO4) nutrient solutions for 16 weeks. CO2 assimilation, starch, soluble carbohydrates, TBARS content and H2O2 production were measured. Transcriptomic analysis of C. sinensis leaves was performed by Illumina sequencing. Our results showed that Mg deficiency decreased CO2 assimilation, but increased starch, sucrose, TBARS content and H2O2 production in C. sinensis leaves. A total of 4864 genes showed differential expression in response to Mg deficiency revealed by RNA-Seq and the transcriptomic data were further validated by real-time quantitative PCR (RT-qPCR). Gene ontology (GO) enrichment analysis indicated that the mechanisms underlying Mg deficiency tolerance in C. sinensis may be attributed to the following aspects: (a) enhanced microtubule-based movement and cell cycle regulation; (b) elevated signal transduction in response to biotic and abiotic stimuli; (c) alteration of biological processes by tightly controlling phosphorylation especially protein phosphorylation; (d) down-regulation of light harvesting and photosynthesis due to the accumulation of carbohydrates; (e) up-regulation of cell wall remodeling and antioxidant system. Our results provide a comprehensive insight into the transcriptomic profile of key components involved in the Mg deficiency tolerance in C. sinensis and enrich our understanding of the molecular mechanisms by which plants adapted to a Mg deficient condition.


Introduction
Magnesium (Mg) is an essential macronutrient for plant growth and production. As the eighth most abundant element in the Earth's crust, soil Mg originates from weathering of minerals containing Mg, such as serpentine, bishopvillite, calcite, dolomite, magnesite, periclase, talcum and diopside. Due to high intensity agriculture, unawareness of the need for Mg fertilizer application and rainy climate, Mg deficiency is a frequently observed nutrient disorder in agricultural production, including the Citrus industry [1,2]. Understanding the physiological and molecular mechanisms underlying plant adaptation to Mg deficiency is of theoretical and practical value for high-quality and high-yield agricultural production.
Mg deficiency increased the content of ethylene and the expression levels of some genes encoding enzymes (1-aminocyclopropane-1-carboxylic acid synthase 2; 7; 8; 11) in the ethylene biosynthetic pathway in model plant Arabidopsis [29], implying that ethylene may play a key role in response to Mg status in plants [30].
Citrus is an evergreen fruit tree that is cultivated in tropical and sub-tropical areas, most of which have a heavy rainfall and high temperature climate. According to our previous survey of 319 Citrus orchards in Fujian province, China, soil acidification was a major problem in these orchards, with an average pH of 4. 34. In addition, 77.4% of soil samples were sub-optimum in exchangeable Mg and 35.6% of leaves were Mg deficient [2]. We have thus far investigated the metabolism of organic acid [31], photosynthesis and the antioxidant system [9], proteomic profile [1], cDNA-AFLP analysis [19] and the metabolism of key metabolites [18] in Citrus plants under Mg deficiency. In order to gain new insight into the molecular mechanism of plant adaption to low Mg condition and ensure green and sustainable development of agriculture, we used the RNA-Seq technique to investigate transcriptomic profile of C. sinensis under long-term Mg deficiency in this study. upregulated genes in leaves being ABA-responsive, even though no change in ABA content was observed [28]. Another study found that Mg deficiency increased the content of ethylene and the expression levels of some genes encoding enzymes (1-aminocyclopropane-1-carboxylic acid synthase 2; 7; 8; 11) in the ethylene biosynthetic pathway in model plant Arabidopsis [29], implying that ethylene may play a key role in response to Mg status in plants [30].

Effects of Mg
Citrus is an evergreen fruit tree that is cultivated in tropical and sub-tropical areas, most of which have a heavy rainfall and high temperature climate. According to our previous survey of 319 Citrus orchards in Fujian province, China, soil acidification was a major problem in these orchards, with an average pH of 4. 34. In addition, 77.4% of soil samples were sub-optimum in exchangeable Mg and 35.6% of leaves were Mg deficient [2]. We have thus far investigated the metabolism of organic acid [31], photosynthesis and the antioxidant system [9], proteomic profile [1], cDNA-AFLP analysis [19] and the metabolism of key metabolites [18] in Citrus plants under Mg deficiency. In order to gain new insight into the molecular mechanism of plant adaption to low Mg condition and ensure green and sustainable development of agriculture, we used the RNA-Seq technique to investigate transcriptomic profile of C. sinensis under long-term Mg deficiency in this study.

Effects of Mg Deficiency on Dry-Weight (DW), Mg Content, Leaf Gas Exchange Parameters and Soluble Sugar Content of C. sinensis
Mg deficiency (-Mg) leads to a symptom of chlorosis and protruding veins in C. sinensis leaf ( Figure 1). Mg deficiency dramatically decreased the DW and Mg content of root ( Figure 2A,D), stem ( Figure 2B,E) and leaf ( Figure 2C,F) in C. sinensis. Mg deficiency decreased CO2 assimilation and increased the intercellular CO2 concentration in C. sinensis leaves ( Figure 3A,B). Chlorophyll a fluorescence measurement indicated that Mg deficiency decreased the maximum quantum yield of primary photochemistry (Fv/Fm), electron transport flux per reaction center at t = 0 (ETo/RC) and increased the minimum fluorescence (Fo) ( Figure 3C-E). Due to the fundamental role of Mg in photosynthate transport, Mg deficiency impaired sucrose transport from source leaves to sink tissues, and thus both starch and sucrose were higher in -Mg leaves than in the control ones ( Figure  3F,G).

Mg Deficiency Increased Reactive Oxygen Species Generation, Membrane Lipid Peroxidation and Lignin Content
Mg deficiency significantly increased reactive oxygen species stress and membrane lipid peroxidation, revealed by the increased production of H2O2 and TBARS content, respectively ( Figure  3H,I). Moreover, in order to testify whether accumulated photosynthates could induce lignin biosynthesis, the lignin content was measured. Data showed that Mg deficiency could actually increase the lignin content in C. sinenesis leaves ( Figure 3J).    4). Difference among the treatments was analyzed by student's t-test. Different letters indicate a significant difference at p < 0.05.

Mg Deficiency Increased Reactive Oxygen Species Generation, Membrane Lipid Peroxidation and Lignin Content
Mg deficiency significantly increased reactive oxygen species stress and membrane lipid peroxidation, revealed by the increased production of H 2 O 2 and TBARS content, respectively ( Figure 3H,I). Moreover, in order to testify whether accumulated photosynthates could induce lignin biosynthesis, the lignin content was measured. Data showed that Mg deficiency could actually increase the lignin content in C. sinenesis leaves ( Figure 3J).

RNA-Seq, De Novo Assembly, Transcripts Annotation and Differentially Expressed Genes (DEGs) Identification
To explore the transcriptomic change of C. sinensis leaves in response to Mg deficiency, the control and Mg deficiency libraries were constructed by using NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Beverly, MA, USA) and cDNA fragments of preferentially 150~200 bp in length were purified with AMPure XP magnetic beads (Beckman Coulter Genomics, Danvers, MA, USA). Each treatment contained two biological replicates and each biological replicate comprised mixed samples from five different plants. We did not sequence three biological replicates considering that we used two treatments of C. sinensis and the data from the control and the Mg-deficient samples could mutually corroborate each other. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and 125 bp/150 bp paired-end reads were generated. By removing reads containing adapter, reads containing ploy-N and low quality reads from raw data, clean reads were eventually obtained. At the same time, Q20 and GC content in the clean data was calculated. The obtained libraries contained total reads ranging from 48,980,872 to 63,560,528, clean reads ranging from 46,926,606 to 60,296,780, clean bases ranging from 6.71 G to 9.04 G and Q20 values ranging from 96.27 to 97.23 (Table 1). Genome mapping revealed that total mapped reads, multiple mapped reads and uniquely mapped reads ranged from 35,701,067 to 44,649,336, 1,256,345 to 1,515,278, and 32,854,616 to 43,134,058, respectively. Here, 71.03% to 73.43% of the clean reads were uniquely mapped to the Citrus genome, which was consistent with our previous transcriptomic analysis (Table 1) [32]. The gene expression levels represented by FPKM showed that the proportion of 0~1 FPKM interval was nearly 50%, followed by 3~15, 15~60 and 1~3 FPKM intervals (Table 2). Interestingly, Mg deficiency increased the proportion of relatively high abundance genes (15~16, >60 FPKM interval) and decreased the proportion of relatively low abundance genes (0~1 FPKM interval). According to the criteria of log 2 (fold change) > 1 and FDR < 0.05, a total of 2176 downregulated and 2688 upregulated genes were identified as DEGs that responded to Mg deficiency in C. sinensis leaves (Tables S1 and S2). Gene ontology (GO) enrichment analysis showed that 28 GO terms with a corrected p < 0.05 were significantly enriched by using the method of Wallenius non-central hyper-geometric distribution (Table S3) [33]. In order to clarify which biological pathway was involved in the Mg deficiency tolerance in C. sinensis, 18 GO terms out of 28 enriched GO terms, which belong to the term type of biological process, would be the main topic of subsequent discussion (Figure 4).  Interestingly, Mg deficiency increased the proportion of relatively high abundance genes (15~16, >60 FPKM interval) and decreased the proportion of relatively low abundance genes (0~1 FPKM interval). According to the criteria of log2(fold change) > 1 and FDR < 0.05, a total of 2176 downregulated and 2688 upregulated genes were identified as DEGs that responded to Mg deficiency in C. sinensis leaves (Tables S1 and S2). Gene ontology (GO) enrichment analysis showed that 28 GO terms with a corrected p < 0.05 were significantly enriched by using the method of Wallenius non-central hyper-geometric distribution (Table S3) [33]. In order to clarify which biological pathway was involved in the Mg deficiency tolerance in C. sinensis, 18 GO terms out of 28 enriched GO terms, which belong to the term type of biological process, would be the main topic of subsequent discussion (Figure 4).

Mg Deficiency Increased Microtubule-Based Movement (GO: 0007018)
Mg plays a vital role in the transport of photosynthate through phloem from source tissue to sink tissue in plants [34]. Among the photosynthates measured in Mg-deficient plants, sucrose was the most pronounced one that accumulated in old matured leaves and resupplying Mg to Mg-deficient plants apparently restored phloem export of sucrose within 12 h [12]. Moreover, evidence suggested that Mg-ATP is the major complex of ATP in biological system and a fall in concentration of Mg-ATP at the phloem-loading site is the major reason for the inhibited sucrose loading in Mg-deficient leaves [14,35]. In the current study, Mg deficiency conspicuously influenced the gene expression level of kinesins with 22 kinesin-encoded genes upregulated and three downregulated (Table S4). According to the available literature, kinesin is a two-headed motor

Mg Deficiency Increased Microtubule-Based Movement (GO: 0007018)
Mg plays a vital role in the transport of photosynthate through phloem from source tissue to sink tissue in plants [34]. Among the photosynthates measured in Mg-deficient plants, sucrose was the most pronounced one that accumulated in old matured leaves and resupplying Mg to Mg-deficient plants apparently restored phloem export of sucrose within 12 h [12]. Moreover, evidence suggested that Mg-ATP is the major complex of ATP in biological system and a fall in concentration of Mg-ATP at the phloem-loading site is the major reason for the inhibited sucrose loading in Mg-deficient leaves [14,35]. In the current study, Mg deficiency conspicuously influenced the gene expression level of kinesins with 22 kinesin-encoded genes upregulated and three downregulated (Table S4). According to the available literature, kinesin is a two-headed motor protein that powers organelle transport along microtubules even in the absence of attachment to a cargo. Further study showed that two heads of kinesin can force the long distance movement along microtubule by hydrolyzing ATP molecules [36]. This gave rise to the possibility that Mg might perform its role by interacting with ATP to form Mg-ATP complex, which was subsequently utilized by kinesin to force cargo movement along microtubule. However, whether the lack of Mg-ATP caused by the depletion of Mg or the accumulation of cargo such as sucrose in C. sinensis leaves is the main reason for the up-regulation of the kinesin genes remains further investigation. Among the differentially expressed kinesin genes, some of them were associated with cell cycle regulation, such as genes encoding chromosome-associated kinesin (KIF4) and kinesin (centromere protein) like heavy chain-like protein (Table S4). This indicated that Mg deficiency may also disturb normal cell replication and proliferation in C. sinensis leaves.

Mg Deficiency Affected Carbohydrate Binding Activities (GO: 0030246)
The accumulation of carbohydrate in the mature leaves more or less impacted upstream and downstream biological process such as starch hydrolysis, amino acid metabolism, secondary compounds synthesis, plasma membrane transport activity [9,13,18,37]. As we expected, Mg deficiency influenced the expression level of an amount of genes related to carbohydrate binding activities compared to control ones (Table S5). Among these DEGs, ten genes were downregulated and 44 genes were upregulated by Mg deficiency. Interestingly, all the gene expression levels of five cold-inducible wall-associated kinases (WAKs) and 11 WAK-like genes were higher in the Mg-deficient leaves than in the control ones. Available evidence showed that WAKs belong to a protein family that was implicated primarily in signal transduction in host response to pathogen attack and particular in the perception of pectin-containing fragments [38][39][40]. Furthermore, Sivaguru et al. [41] found that WAK1 was involved in aluminum stress in Arabidopsis and thus may participate in heavy metal response. This result indicated that WAKs/WAK-like genes might participate in a global response to environmental stimuli including Mg deficiency. Similar to WAKs, Lectin-domain containing receptor protein kinase (LecLPK) is believed to play crucial roles in saccharide signaling as well as stress perception in plants [42]. Expression of the AtLecLPK1 was strongly induced by ABA, methyl jasmonate (MJA), salicylic acid (SA) and stress treatments including salt stress in Arabidopsis [43]. Here, three upregulated LecLPKs were isolated from -Mg C. sinensis leaves (Table S5 (Supplementary Materials)). This result corroborated that Mg deficiency, as a nutrient disorder, like other environment stimuli, could trigger a cascade of events along signal transduction and stress perception. Small glutamine-rich tetratricopeptide repeat (TPR)-containing protein (SGT) alpha (SGTA) belongs to a family of molecular co-chaperone proteins that harbored a TPR motif, acting as a protein-protein interaction module, which was capable of facilitating interaction with a diverse range of client proteins in animal [44]. In plants, except for a reference reporting SGT beta subunit that was associated with chromium (Cr) tolerance in spring wheat, limited data about SGT led to further studies on the function of this gene [45]. The upregulated SGTA in -Mg C. sinensis leaves means that SGTA may also respond to nutrient disorder in plants. Furthermore, among DEGs enriched in carbohydrate binding activities, genes encoding 2-acetamido-2-deoxy-d-galactose-binding seed lectin 2, 3-hexulose-6-phosphate isomerase, chloroplast post-illumination chlorophyll fluorescence increase protein isoform 1(fragment), glucoamylase and uncharacterized ribonuclease sll1290, which are implicated in plant defense against pathogen and herbivory [46], assimilating endogenous gaseous formaldehyde [47], chlororespiration in protecting plastoquinone (PQ) from over-reduction and thereby also photosystem (PS) I from photoinhibition [48], decomposition of starch [49], and posttranscriptional regulation of gene expression [50], respectively, were downregulated by Mg deficiency in C. sinensis leaves. These RNA-Seq data were consistent with our physiological and chlorophyll a fluorescence measurement, which showed that -Mg leaves accumulated a higher content of starch and suffered more severe photoinhibition ( Figure 2C-E). Brassinosteroid (BR) insensitive 1-associated receptor kinase, also known as leucine-rich repeat receptor like kinase, positively regulated plant growth, development and response to biotic and abiotic stresses by directly binding to the BR ligand, triggering a signal cascade in the cytoplasm that leads to the transcription of BR-responsive genes [51]. The upregulated gene expression levels of several BR insensitive 1-associated receptor kinases in C. sinensis leaves might imply that BRs participated in the response to Mg deficiency in Citrus plants (Table S5). Phosphorylation, especially protein phosphorylation, is a crucial biological regulation in plant growth and development. Our transcriptomic analysis revealed that Mg deficiency altered the expression of 439 genes involved in phosphorus metabolic process, of which, 337 and 309 genes were related to phosphorylation process and protein phosphorylation, respectively (Table S6).
Protein phosphorylation is the most abundant post-translational modification (PTM) of functional proteins and a key regulatory step in the eukaryotic signal transduction [52]. It was supposed to govern various processes such as CO 2 assimilation [53], chromatin stabilization [54], Ca signaling [55], hormone signaling [56], protein trafficking [57], protein turnover [58], photosynthate transport [59], circadian regulation [60] and pathogen-perception [61]. Here, among the 309 DEGs involved in protein phosphorylation (Table S6), 234 and 75 DEGs were upregulated and downregulated in C. sinensis leaves by Mg deficiency, respectively. Interestingly, more than one-third of the 309 DEGs belong to the serine-threonine protein kinase family, which is consistent with the concept that most phosphorylation events occur in serine or threonine residues in plants [52]. Meanwhile, as we expected, the conserved genes involved in protein phosphorylation, such as Protein kinase APK1B (PK), calcium-dependent protein kinase (CDPK), CBL-interacting protein kinase (CIPK), mitogen-activated protein kinase (MAPK)/kinase (MAPKK) and LRR protein kinase were differentially regulated by Mg deficiency. This result indicated that the sequential activation or inactivation of the downstream cascade may result in the changed processes of cell division, growth, differentiation, programmed cell death and environmental stress response [62]. Such expression alteration of the genes related to protein phosphorylation was also observed in Citrus and tobacco treated with boron deficiency [63,64], Poncirus trifoliata and soybean treated with phosphorus deficiency [65,66], rice treated with calcium deficiency [67], and cassava treated with low-temperature [68]. Evidence showed that G-type lectin S-receptor-like serine/threonine protein kinase (GsSRK) was a positive regulator of tolerance to salt stress in wild soybean, Gossypium barbadense and Arabidopsis [68,69]. Additionally, Vaid et al. [70] also found that GsSRKs are highly induced under abiotic stress conditions including heat, drought, salt, and cold in Arabidopsis and rice. Here, 18 differentially expressed GsSRKs were isolated from C. sinensis leaves and most of them were upregulated by Mg deficiency (Table S6). Besides protein phosphorylation, various genes involved in substrate phosphorylation were also altered in C. sinensis leaves by Mg deficiency (Table S6). Multidrug resistance-associated protein 2,6 (MRP2,6) can use the energy generated by the hydrolysis of Mg-ATP to facilitate the transmembrane movement of a variety of small molecules [71]. Here, we found that the expression level of MRP2,6 was downregulated by Mg deficiency in C. sinensis leaves (Table S6). This may further confirm that Mg deficiency hinders the transport of cargos including sugar due to the decreased Mg-ATP as we mentioned above. Interestingly, the gene encoding phototropin-1, which is an important mediator in the phototropic growth and stomata on-off in plants, was also downregulated by Mg deficiency (Table S6) [72], indicating that C. sinensis can self-adapt to relative high photo-radiation to ameliorate the photo-oxidative and photo-inhibition damage when the efficiency of photosynthetic electron transport and photosynthesis was impaired by Mg deficiency [9,19]. Furthermore, the genes encoding clathrin interactor 1, cyanidin-3-O-glucoside 2-O-glucuronosyltransferase, gamma-glutamyl hydrolase, lipase, L-lactate dehydrogenase, pyruvate kinase (cytosolic), ribonucleoside-diphosphate reductase small chain, etc., were upregulated in C. sinensis leaves by Mg deficiency (Table S6). The alteration of these genes implies that the endocytic pathway, glycolysis, anthocyanins modification, lipid metabolism, chromatin organization, etc., may be altered in response to Mg deficiency [31,73,74]. Genes encoding the proteins mentioned above should be further validated to understand their physiological role in tolerance to Mg deficiency in the future study. Isoprenoids are a fascinating family of compounds derived from the C5 precursor isopentenyl diphosphate (IPP) and dimethylallyl diphosphate (DMAPP), which contain primary and secondary metabolite [75]. Isoprenoids function in respiration (ubiquinone), photosynthesis (carotenoid, chlorophyll, tocopherol, plastoquinone), membrane architecture (sterol) and growth regulation (brassinosteroid, cytokinin, gibberellin, abscisic acid, strigolactone), whereas the secondary compounds mediate important interaction between plants and their environment [76]. In plants, there are two separate and biochemically different IPP biosynthesis pathways, named the mevalonic acid (MVA) and methylerythritol 4-phosphate (MEP) pathway. The MEP pathway simultaneously produces both IPP and DMAPP from pyruvate and glyceraldehyde 3-phosphate (GAP) in the plastids, whereas the MVA pathway synthesizes IPP from acetyl-CoA in the cytoplasm [77,78]. Here, we found that genes encoding geranylgeranyl pyrophosphate synthase (3 genes), probable 1-deoxy-d-xylulose-5-phosphate synthase (chloroplastic, 2 genes), putative geranyl diphosphate synthase (1 gene) and 2-C-methyl-d-erythritol 2,4-cyclodiphosphate synthase (chloroplastic, 1 gene), which involved in MEP pathway, were downregulated in C. sinensis leaves by Mg deficiency (Table S7). Furthermore, the generation of certain isoprenoid derivatives such as chlorophyll, ubiquinone, lycopene and zeaxanthin, may also be hindered by Mg deficiency due to the downregulated gene expression of protoporphyrinogen oxidase (chloroplastic), electron transfer flavoprotein-ubiquinone oxidoreductase, lycopene epsilon cyclase (chloroplastic), lycopene beta cyclase (chloroplastic/chromoplastic), zeaxanthin epoxidase (chloroplastic) (Table S7). Interestingly, the expression level of glycerol-3-phosphate dehydrogenase (SDP6, mitochondrial) was upregulated by Mg deficiency, which could result in the increasing conversion of dihydroxyacetone phosphate (DHAP) to glycerol-3-phosphate (G3P) and decrease DHAP to GAP shunt [79]. This result may further support the point that Mg deficiency decreased MEP-mediated isoprenoids biosynthesis in C. sinensis leaves. Such down-regulation of the MEP pathway was also observed in poplar leaves under UV-B radiation [80]. The reduction in the MEP pathway under Mg deficiency might result from elevated intercellular CO 2 ( Figure 3A), which may enhance higher consumption rate of cytosolic PEP (phosphoenolpyruvate) through increased PEP carboxylase (PEPC) activity [31], thus lowering the rate of PEP transport into the chloroplast, where PEP, in its dephosphorylated form (pyruvate), feeds into the MEP pathway [81]. Although understanding of the metabolic regulation of the MEP pathway has emerged in the past decade, knowledge about the dynamic change of the regulation steps of the MEP pathway's response to environmental stimuli, especially nutrient disorders, is limited. Future study should be aimed at a complete understanding of the molecular regulation of the MEP pathway under Mg deficiency as well as other nutrient disorder.

Mg Deficiency Affected the Genes Involved in Lipid Metabolism (GO:0006629)
Lipids are one of the major components of the plasma membrane, which is the interface between the cell and the environment [82]. Available literature showed that environmental stimuli could alter lipid content and lipid metabolism in higher plants. For instance, the content of polar lipid compounds in Arabidopsis and Virgilia divaricata was decreased by UV-B radiation [80], heat stress [83], P deficiency [84] and high light [14]. Moreover, RNA-Seq and proteomic studies showed that the genes or functional proteins involved in lipid metabolism could be altered by salinity stress in barley [85], low temperatures in cassava [86], P deficiency in Arabidopsis [87], boron deficiency in C. sinensis [63], and Mg deficiency in C. reticulata [19] and C. sinensis leaves [88]. Our previous studies showed that the expression levels of gene and protein abundance of phosphatidylinositol 3-kinase and 4-kinase were simultaneously upregulated by Mg deficiency in C. sinensis leaves, thus contributing to the Mg-deficiency-tolerance of plants [1,88]. Here, we successfully identified and quantified 194 DEGs related to lipid metabolism in C. sinensis leaves, of which, 86 DEGs were upregulated and 108 DEGs were downregulated by Mg deficiency (Table S8).
Phosphoinositide phospholipases C (PI-PLC) are the principle enzymes, which catalyze the initial step of phospholipid breakdown in cell membrane and generate multiple lipid derived second messengers such as diacylglycerol (DAG) and inositol 1,4,5-trisphosphate (IP3), which are involved in a wide range of cellular processes including plant development and stress response [82]. The alternation of PLC gene expression levels by environmental stimuli have been reported in Arabidopsis under heat [83], cold [89], dehydration and salt stress [90]. In the current study, Mg deficiency upregulated the expression levels of two PLC2 genes and downregulated two PLC isoform (PLC4 and PLC6) genes in C. sinensis leaves. Vossen et al. [91] also found that the different PLS genes in tomato leaves showed distinct expression patterns in response to pathogen infection. Previous research indicated that under phosphorus-limiting condition, membrane phospholipids provided a pool for inorganic phosphate by the activity of intracellular phospholipase C, which can be used for the synthesis of other essential phosphorus-containing biomolecule in Sinorhizobium meliloti [92]. The up-regulation of PLC2 genes may also be explained partially by this way, as Mg deficiency significantly decreased the P content in C. sinensis leaves [93]. Besides PLCs, other form of phospholipase such as patatin-like phospholipase family protein, phospholipase A1-I gamma and phospholipase A1-II gamma, were also differentially regulated by Mg deficiency (Table S8).
Fatty acids (FA), a major source of energy in plants, are the main components of plant membrane lipids. The increasing of membrane unsaturated FA catalyzed by fatty acid desaturase, plays a critical role in ion permeability, appropriate fluidity and stability of membrane [94]. The expression of a fatty acid desaturase gene (stearoyl-ACP desaturase SSI2) was upregulated almost 10 fold under cold stress in cassava apical shoots, implying that this gene may play a role in enhancing membrane fluidity [86]. Here, we report that genes encoding Omega-3 fatty acid desaturase (endoplasmic reticulum) and acyl-acyl-carrier-protein desaturase (chloroplast) and genes encoding the sphingolipid long chain base delta 8 desaturase, temperature-sensitive omega-3 fatty acid desaturase (chloroplast), omega-6 fatty acid desaturase (chloroplast) and omega-6 fatty acid desaturase isozyme 2 (endoplasmic reticulum) were upregulated and downregulated in C. sinensis leaves by Mg deficiency, respectively (Table S8). The differentially regulated FA desaturase genes means that fine-tuning of the saturation of fatty acids in different subcellular organelles may confer resistance to Mg deficiency and subsequent side effects downstream in C. sinensis, although such resistance could not compensate for the adverse injury as the TBARS content was increased in C. sinensis leaves by Mg deficiency ( Figure 3I).
Very long chain fatty acids (VLCFAs) are involved in diverse biological functions such as membrane constituent, surface barrier and seed storage compound. The first step of VLCFAs biosynthesis is the condensation of two carbons to an acyl-CoA, which is catalyzed by 3-keto acyl-CoA synthase (KCS) [95]. Here, we found three KCS genes (3-ketoacyl-CoA synthase 11, 17 and 21) were downregulated by Mg deficiency (Table S8). This is consistent with the results obtained from C. sinensis under B deficiency [63], rice under Mg deficiency [96], Arabidopsis under drought stress [97] and flax (Linum usitatissimum L.) under saline-alkaline stress [98]. The down-regulation of KCS genes may indicate that the biosynthesis of VLCFAs was slowed down by adverse conditions including Mg deficiency. Besides KCS genes, the differential regulation of genes encoding acyl-acyl-carrier-protein thioesterase type B, acyl-coenzyme A oxidase, acyl-protein thioesterase 2, lecithin-cholesterol acyltransferase and protein WAX2 (fatty acid hydroxylase domain containing CER1-like protein) by Mg deficiency implied that the balance of fatty acid metabolism might be impaired in C. sinensis leaves (Table S8). Interestingly, a gene encoding the probable mannitol dehydrogenase, which catalyzes the conversion of mannitol to mannose, was upregulated by Mg deficiency in C. sinensis leaves. Previous research indicated that mannose could induce stomatal closure through methods mediated by ROS production [99]. This is accordant with the higher level of ROS production and stomatal closure represented as the H 2 O 2 production and stomatal conductance in C. sinensis leaves ( Figure 3H) [18].

Mg Deficiency Affected the Genes Involved in Metabolic Process (GO: 0008152) and Biological Process (GO: 0008150)
A total of 2113 DEGs were clustered into metabolic process with 1142 DEGs and 971 DEGs upregulated and downregulated by Mg deficiency, respectively, whereas a total of 609 more DEGs were clustered into biological process than into metabolic process. Because some of these DEGs overlapped with those clustered into the GO terms discussed above, hereafter we mainly focus on the DEGs involved in the conspicuously altered pathway that differed with those discussed above.
Plant hormones simultaneously coordinate physiological responses to biotic and abiotic stresses [100].
Here, we found that the expression levels of key genes such as S-adenosylmethionine synthase (SAMS), 1-aminocyclopropane 1-carboxylate synthase (ACCS) and 1-aminocyclopropane-1-carboxylate oxidase (ACCO), which are involved in ethylene biosynthesis pathway, and several ethylene-responsive transcription factors (ERFs) were all significantly upregulated by Mg deficiency (Table S9). Available studies show that ethylene interacts with nutrient uptake (e.g., P and N) and controls plant response under growth-limiting condition or stress. Environmental stress could speed up the ethylene production rate in plants [101]. The up-regulation of genes involved in ethylene biosynthesis might imply that Mg deficiency enhances ethylene production of C. sinensis leaves, which is consistent with the results obtained in Arabidopsis under Mg deficiency [29]. Besides the DEGs related to ethylene, genes encoding indole-3-acetic acid inducible 3, indole-3-acetic acid inducible 9, auxin efflux carrier component, auxin-induced proteins and auxin-responsive proteins, which are involved in auxin response, were also induced by Mg deficiency in C. sinensis leaves (Table S9). The up-regulation of transport proteins/genes and adjustment of metabolic processes of auxin may help root hair elongation and enhance nutrient absorption under adverse conditions [66,101]. According to the up-regulation of auxin related genes, we speculate that Mg deficiency might also enhance auxin production in C. sinensis leaves.
As we described above, Mg deficiency significantly decreased CO 2 assimilation and F v /F m , ET o /RC and increased the intercellular CO 2 concentration and F o in C. sinensis leaves ( Figure 3A-E). This result is highly concordant with RNA-Seq data, which show that genes encoding chloroplast apparatus such as chloroplast oxygen-evolving complex (OEC)/thylakoid lumenal 25.6 Da protein, chloroplast photosynthetic water oxidation complex 33 kDa subunit, chloroplast ribose-5-phosphate isomerase and photosystem I (PS I) reaction center subunits, thioredoxins (THX), protoporphyrinogen oxidase (PPO) and photosystem II (PS II) core complex proteins were downregulated by Mg deficiency in C. sinensis leaves (Table S9). One of the apparent symptoms of Mg deficiency is the thickening, suberification and even splitting of leaf veins due to the accumulation of callose or lignin (Figure 1; Figure 3J). Here, we found that several DEGs encoding cellulose synthase A, callose synthase, laccase and cytochrome P450, which are responsible for the biosynthesis of cellulose, callose or lignin, were upregulated by Mg deficiency in C. sinensis leaves (Table S9).
Transcription factors tightly control plant development, nutrient acquisition, biotic and abiotic stresses [102]. Here, we found that the gene expression levels of two AP2, two bHLH and seven MYB containing transcription factors were downregulated, meanwhile, the gene expression levels of five MYB transcription factors, 19 NAC domain-containing proteins and 11 WRKY transcription factors were upregulated in C. sinensis leaves ( Figure 5; Table S9). This result was in accordance with the research of Hermans et al. [28], who also found that Mg deficiency altered the expression levels of several transcription factors and upregulated the gene expression level of AP2 in Arabidopsis. Recent evidence showed that MYB transcription factor could act as a sensory switch regulating lignin biosynthesis in woody plant cells [103]. The up-regulation of MYB transcription factor may also involve in the accumulation of lignin content in C. sinensis leaves under Mg deficiency ( Figure 3J). Mg deficiency decreased light energy utilization and phloem loading of photosynthate, thus increased the oxidative stress and the accumulation of photosynthate in source leaf cells [9,22]. As a strategic response, DEGs involved in cell wall remodeling such as xyloglucan endotransglucosylase/hydrolases (XTH), and the antioxidant system such as peroxidase (POD), S-adenosylmethionine-dependent methyltransferase (SAM-MT), glycoside hydrolase (GH), cytochrome P450 (CYP450) were upregulated by Mg deficiency in C. sinensis leaves, indicating their probable roles in cell wall loosening and antioxidation under Mg deficiency ( Figure 5; Tables S9 and S10) [19]. A previous study showed that Mg deficiency led to leaf cell senescence and necrosis in rice and bean, respectively [104]. Local necrosis may facilitate the infection of pathogen and bacteria, therefore, pathogen and disease resistant related genes such as pathogenesis-related protein 5, pathogenesis-related protein 10, TIR-NBS-LRR-TIR type disease resistance protein, chitinase-like protein 2, late blight resistance protein R1-A and TMV resistance protein N, etc., were upregulated in C. sinensis leaves under Mg deficiency (Tables S9 and S10).

Plant Culture and Mg Treatment
Uniform seeds of C. sinensis were sown in a plastic tray containing clean river sand and kept moist under natural light and temperature in the greenhouse of Fujian Agriculture and Forestry University. After the seeds germinated, tender seedlings were irrigated with 1/4 strength nutrient solution every other day. Seven weeks later, uniform seedlings with two leaves and one sprout were transplanted to 6 L pottery pots containing clean river sand and irrigated with 500 mL 1/2 strength nutrient solution every two days. There were two seedlings per pot. The full strength nutrient solution contained the following macronutrients (in mM): KNO 3

Leaf Gas Exchange and Chlorophyll a Fluorescence
Leaf gas exchange was measured by using a CIRAS-2 portable photosynthesis system (PP systems, Amesbury, MA, USA) at an artificial CO 2 concentration (380 µmol mol -1 ) supplied by a CO 2 cylinder and a controlled light intensity of 1000 µmol m −2 s −1 between 10:00 am and 11:30 am on a clear day. During measurements, leaf temperature and vapor pressure deficit were 30.1 ± 0.2 • C and 1.37 ± 0.1 kPa, respectively. Chlorophyll a fluorescence parameters were measured by using Handy PEA parable instrument (Hansatech, Norfolk, UK) on the dark-adapted leaves.

Plant Dry-Weight (DW), Mg and Leaf Soluble Carbohydrates Content
At the end of experiment, ten plants per treatment from different pots were separated into roots, stems and leaves, loaded into kraft bags and oven dried until at a constant weight at 70 • C. After DW was measured by an electronic balance, plant tissues were ground into fine powder. Mg concentration in roots, stems and leaves was measured by atomic absorption spectroscopy after dry ashing and digestion with 1 M HCl [107]. Sucrose, fructose, glucose and starch in leaves and roots were extracted and measured according to the methods described by Yang et al. [9] and Jones et al. [108]. Briefly, about 100 mg samples were extracted three times with 80% (v/v) ethanol at 80 • C. After the extracts were evaporated, the resulting pellets were thoroughly dissolved with 3 mL ddH 2 O. Glucose was measured in 1 mL reaction mixture of 100 mM inidazole-HCl (pH 7.9), 5 mM MgCl 2 , 0.5 mM NAD, 1 mM ATP, 2 units G6PDH and 2 units hexokinase. Fructose and sucrose was measured at the same reaction solution mentioned above by adding 2 units phosphoglucose isomerase and 20 units invertase, respectively. Starch was determined by enzyme kinetic method as glucose equivalents with G6PDH and hexokinase.
4.4. Measurements of H 2 O 2 Production, TBARS and Lignin Content in C. sinensis leaves H 2 O 2 production was measured according to Yang et al. [109]. The content of H 2 O 2 was calculated with an extinction coefficient of ε = 26.6 cm −1 mM −1 . TBARS was extracted and measured according to the methods described by Hodges et al. [110]. For lignin measurement, leaf samples were homogenized in 95% ethanol. After centrifugation at 1500 g for 5 min, the pellet was washed three times with 95% ethanol and twice with ethanol-hexane (1:2, v/v). The washed pellet was allowed to air-dry. The lignin content was determined according to the method described by Morrison with some modifications [111]. There were four replicates for H 2 O 2 , TBARS and lignin contents.

Total RNA Extraction and RNA-Seq
Samples from five different plants were pooled into one biological replicates that used to RNA extraction and further sequence analysis. There were two biological replicates per treatments. Total RNA were extracted from leaves and roots samples by using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Degradation and contamination of total RNA were monitored by 1% agarose electrophoresis. The purity and integrity of total RNA were measured by Nanodrop 2000 and Qubit 4 Fluorometer (ThermoFisher, New York, NY, USA). Only the intact RNA samples with the ratio of A260 (absorbance at 260 nm) to A280 between 1.8 and 2.0 were submitted for further analysis. Message RNA (mRNA) was enriched from 3 µg total RNA by using magnetic beads with oligo dT (Qiagen, Valencia, CA, USA). The resulting mRNA was then fragmented into short fragments (200 nt) and used for first-strand cDNA synthesis with random hexamer-primers. Subsequently, second strand cDNA was synthesized by using DNA polymerase I and purified by AMPure XP beads (Beckman Coulter, Beverly, MA, USA). Purified second-strand cDNA were then put through end repair, 3' end adenosine tailing, sequencing adaptor connection, fragment selection by AMPure XP beads and finally PCR amplification. The library quality was assessed using the Agilent Bioanalyzer 2100 system (Agilent, Santa Clara, CA, USA). After cluster generation, the library amplicons were sequenced on the Illumina HiSeq 4000 platform (Illumina Inc., San Diego, CA, USA) using the paired-end technology in a single run. Illumina GA Pipeline (Version 1.6) was used to perform the original image process to sequences, base calling and quality value calculation, in which 125 bp/150 bp paired-end reads were generated.

RNA Reads Mapping and Analysis of DEGs
Clean reads were obtained by removing adapter reads, reads with more than 10% N (uncertain base) and low quality reads (the proportion of Q20 more than 50%). After getting rid of the reads mapping to ribosomal RNA, the remaining reads were mapped to the C. sinensis genome published by Huazhong Agricultural University (http://citrus.hzau.edu.cn/orange/) by TopHat2 software [112]. The gene expression level was calculated by FPKM (Fragments Per Kilobase of transcript per Million mapped reads) method [113]. Genes with an adjusted p-value < 0.05 and fold change (based on log 2 ) >1 or <−1 identified by DESeq software were considered DEGs. Functional annotation of DEGs was carried out by the method described by Young et al. [33] and Kanehisa et al. [114].

RT-qPCR Analysis of DEGs
Total RNA extraction and first-strand cDNA synthesis were performed by the methods described above. There were three biological replicates and each replicate contained pooled samples from five different plants (one plant per pot). Gene special primer pairs were designed by Premier Primer 5.0 software (Premier Biosoft International, Palo Alto, CA, USA) according to the sequences published in the citrus genome mentioned above (Table S11). RT-qPCR was carried out in a mixture containing 10 µL Bestar TM qPCR MasterMix SYBR Green (DBI Bioscience, Shanghai, China), 0.4 µL of 10 µM forward primer, 0.4 µL of 10 µM reverse primer, 2 µL cDNA diluted template and 7.2 µL ddH 2 O. RT-qPCR was performed by using a CFX96 Touch TM Deep Well Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). The cycling conditions were 30 s at 95 • C, followed by 40 cycles of 95 • C for 10 s, 60 • C for 20 s, 72 • C Cr 20 s. Relative gene expression was calculated by using the ddCt algorithm. For the normalization of different samples and replicates, the polyubiqutin gene (Cs2g20650) was used as an internal standard and the control samples were used as references whose expression levels of each gene were set to 1.

Experimental Design and Statistical Analysis
There were 20 pots for each treatment in a completely randomized design. Experiments were conducted with 3-10 replicates (one plant per replicate). There were ten replicates for plant biomass, five replicates for Mg contents, three replicates for RT-qPCR analysis, four replicates for TBARS and H 2 O 2 content measurement, respectively. Results were displayed as means ±SD for n = 3-10. Means were separated by the student's t-test at p < 0.05.

Conclusions
Mg deficiency decreased CO 2 assimilation, but increased starch, sucrose, TBARS contents and H 2 O 2 production in C. sinensis leaves. A total of 4864 genes showed differential expression in response to Mg deficiency as revealed by RNA-Seq and the transcriptomic data were further validated by real-time quantitative RT-qPCR. GO enrichment analysis indicated that the mechanisms underlying Mg deficiency tolerance in C. sinensis may be attributed to the following aspects: (a) enhancing microtubule-based movement and cell cycle regulation (kinesin, Sugar carrier protein C); (b) elevating signal transduction in response to biotic and abiotic stimuli (WAK, LecLPK, SGTA, etc.); (c) altering biological processes by tightly controlling phosphorylation especially protein phosphorylation (PK, GsSRK, LRRPK, MAPKKK, etc.); (d) down-regulating light harvest and photosynthesis process due to accumulation of carbohydrates (OEC, PS II reaction center W, PPO, THX, etc.); (e) up-regulating cell wall remodeling and antioxidant system (XTH, POD, SAM-MT, GH, CYP450, etc.) (The full name of acronyms could be found in Table S12). Our results provide a comprehensive insight into the transcriptomic profile of key components involved in the Mg deficiency tolerance in C. sinensis and enrich our understanding of the molecular mechanisms by which plants adapted to Mg deficient condition.  Table S6. DEGs enriched in phosphorus metabolic process (GO:0006793, GO:0016310 and GO:0006468). Table S7. DEGs enriched in isoprenoid, terpenoid and carotenoid metabolic process (GO:0006720, GO:0006721 and GO:0016116). Table S8. DEGs enriched in lipid metabolism (GO:0006629). Table S9. DEGs enriched in metabolic process (GO: 0008152). Table S10. DEGs enriched in biological process (GO: 0008150). Table S11. Primer pairs used for RT-qPCR analysis. Table S12. List of the acronyms of some nouns.