Transcriptomic analyses reveal physiological changes in sweet orange roots affected by citrus blight

Background Citrus blight is a very important progressive decline disease of commercial citrus. The etiology is unknown, although the disease can be transmitted by root grafts, suggesting a viral etiology. Diagnosis is made by demonstrating physical blockage of xylem cells that prevents the movement of water. This test was used to identify symptomatic trees from four commercial groves in Florida. Total RNA extracts of phloem-enriched scaffold root tissues were prepared from seven trees that failed to take up water and from one healthy tree. These RNA extracts were used for transcriptomic analyses using paired end RNA-Seq from an Illumina 2500 system. The expression of transcripts annotated as polyprotein of citrus endogenous pararetrovirus were estimated by both RT-qPCR and RNA-Seq. Results Transcripts from seven RNA-Seq libraries from trees affected by citrus blight were compared to a control tree. 129–148 million RNA fragments (two paired-end reads/fragment) were generated per library and were mapped to the sweet orange reference genome. In response to citrus blight stress, genes encoding aquaporins, proteins with water channel activity and several cellulose synthase genes were down-regulated, whereas genes involved in lignin and glucosinolate biosynthesis were up-regulated. Transcripts encoding proteins in pathways of carbohydrate metabolism, nucleotide synthesis, signaling, hormone metabolism, secondary metabolism, transport, and biotic stress pathways were overwhelmingly down regulated in all libraries. Conclusion Reduced water intake and xylem plugging were observed in the trees tested and the changes in their transcriptome were analyzed. Plants adapted to reduced water flow by regulating primary and secondary metabolism, nuclear transport and hormone associated pathways. The patterns of energy generation, transcription, translation and protein degradation were consistent with irreversible decline. The down regulation of cellulose synthase transcripts and up regulation of transcripts related to lignin production likely lead to an imbalance in the pathways leading to wood formation, and may lead to the blockage of the xylem vessels seen as the cardinal symptom of citrus blight. Transcripts of a pararetrovirus were elevated in the transcriptome of roots used in this study.


Introduction
Citrus blight (CB) was described in Florida more than a century ago and 50 years of research on CB was reviewed in 1936 [41] and more recently [16]. CB is found in tropical or subtropical regions, notably Florida and Brazil but including the Caribbean, and portions of South Africa and Australia [9,29]. CB is noticeably absent from California, Texas and Asia. The disease affects only bearing trees, primarily grapefruit (Citrus paradisi Macf.) and sweet orange (Citrus sinensis L.) and has not been reported in greenhouse grown trees. Symptoms begin with a distinct loss of reflective sheen on the surface of leaves, a mild wilt, and zinc-deficiency symptoms in the foliage. Trees rapidly decline with extensive twig dieback, small fruit and off-season flowering [3,14] and levels of zinc become elevated in the wood [51]. CB has caused economic losses in excess of $60 million in Florida and was formerly the most significant plant health problem of Florida citrus with a loss of 650,000 trees/ year [48], but with the unprecedented crisis of huanglongbing (HLB) throughout the world, research efforts on CB have diminished.
Trees affected with blight will not recover from the disease if they are severely pruned, which can temporarily reduce the symptoms of HLB. Trees affected by both HLB and CB may decline and die more rapidly than when affected by either disease alone, since the physiology of both the phloem (HLB) and xylem (CB) are disrupted simultaneously. Transmission of blight has been demonstrated by root grafts [37,46,49], but not by bud or approach grafts or through soil [47]. Efforts to use epidemiological models to gain insight into whether CB is a disease or a physiological disorder have been inconclusive, as the models provide conflicting results. In one study CB spread within a grove following a linear model not typical of a pathogen and vectored disease [12], but in another study a logistic model typical of a vector transmitted disease was observed at high levels of disease [11]. Several causal agents for CB have been proposed, including Xylella fastidiosa, Fusarium solani, and an idaeovirus [17,23,27] as well as some abiotic factors, but none of these causes have been confirmed and the etiology of the disease remains unknown. Therefore, neither the isolation of a pathogen nor the detection of its presence in pre-symptomatic trees with PCR-based methods are possible. Confirmatory diagnosis of the disease is primarily based on demonstration of the blockage of water uptake by attempting to inject water into the trunks of affected trees by using a drill to make a hole in the trunk and a syringe to inject the water [33]. Water cannot be injected into the trunk of trees affected by CB. The presence of amorphous plugs [8,9] that physically fill the lumen of the xylem cells and measurement of elevated levels of zinc and soluble phenolics in the wood [51] are also used for diagnosis of CB. An immunoassay for a protein associated with citrus blight in bark tissues has also been used for presumptive diagnosis [6]. Recently, pararetrovirus transcripts have been found in both roots and leaves of blighted trees through Next Generation Sequencing (NGS) [43].
Because the etiology of CB remains unknown and there are no control methods or therapeutic treatments for CB, there is an urgent need for new data to provide both clues to a potential etiological agent and to provide insight that may lead to improved control and therapeutic measures. RNA-Seq has been widely used to provide insight into the pathophysiology of different plant diseases, including CB [55]. RNA-Seq has been particularly useful in the case of HLB and citrus sudden death disease, the causal agents of which are nonculturable [2,19,20,38]. Distinct plant defense responses can be revealed, with different patterns of regulation in response to different pathogens and abiotic stress. We wanted to know how alterations in the transcriptome correspond with the symptoms and physiological changes in advanced stages of CB.

Demonstration of xylem plugging in blight-affected trees
Symptoms observed on the CB affected trees included severe twig dieback and defoliation, as seen on Hamlin sweet orange grafted on citrumelo rootstock (X Citroncirus spp.) (Fig. 1a). In one of the locations, very few trees in the citrus block remained for sampling, as the other trees had already been removed due to the severity of the citrus blight problem (Fig. 1b). The ability of the trees to take up water into their trunks by syringe injection was measured to select trees from which root samples were taken to prepare RNA extracts. Five of the seven trees sampled failed to take up any water in the assay, and the remaining two trees took up only 4 and 10 ml compared to more than 15 ml observed in healthy trees in a 30 s test (Table 1) [33].

Overview of the RNA-Seq data and transcriptome analysis
Each of the RNA-Seq technical replicate libraries was on average 7.25 Mb composed of 71.8 million reads, with 93.6% of reads with Q30 scores greater than 30 and mean quality of 35.55 (Additional file 1: Table S1). Two paired-end reads defined each fragment and four technical replicate libraries were pooled to become the final libraries for each of the eight trees tested ( Table 1). The Citrus sinensis transcripts were annotated using the sweet orange genome [52] and assigned to different Mapman functional bins by Mercator 3.6 (Additional file 2: Figure S1). Although precautions were taken to limit RNase activity during sampling and extraction of RNA from affected roots, the RNA integrity number (RIN) [44] for the seven libraries varied from 1 to 8. The percentage of RNA fragments that were successfully mapped to the citrus genome also varied from 42.6-95.1% (Table 1). When the seven libraries were integrated, 707 transcripts were co-regulated as compared to the control ( Table 2, Fig. 2). The general patterns of expression were highly consistent in all seven libraries.
Within the seven libraries, the co-differentially expressed transcripts (DETs) related to the biological categories of metabolic processes (GO:0008152) and cellular processes (GO:0009987) were most disrupted, followed by response to stimulus (GO:0050896) (Additional file 3: Figure S2A). Within the cellular component category, subcategories associated with cell parts, membranes and organelles were heavily affected (Additional file 3: Figure S2B). The metabolic pathways that were most heavily affected by citrus blight included pathways related to the cell wall, hormone signaling, proteolysis, transcription factors, and secondary metabolism. Pathogenesis related (PR-protein) genes and associated pathways were strongly affected and overwhelmingly repressed (Fig. 3).

Primary and secondary metabolism
In all seven libraries, co-DETs involved in primary metabolism, including carbohydrate, lipid, amino acid and protein, as well as secondary metabolism were overwhelmingly down-regulated (Additional file 4: Table S2). Both sucrose and starch metabolism were especially downregulated (Additional file 4  [20]. Genes encoding PRP4 (proline-rich protein), XTR4 (xyloglucan endotransglycosylase-related protein) and pectin acetyl esterase protein, both involved in  cell wall modification were also induced ( Table 3). Concanavalin A-like lectin protein involved in post translational modification was also induced ( Table 3). Genes encoding a ubiquitin E3 RING zinc finger protein, an aspartyl protease and AATP1 (AAA-ATPase 1), associated with protein degradation in response to stress were also up-regulated (Table 3).
In contrast to the cellulose synthase genes that were downregulated, two genes involved in lignin biosynthesis and a suite of transcripts encoding aromatic monooxygenase CYP79A2 for sulfur containing aromatic glucosinolates were very sharply upregulated (Table 3

Energy and transport
Genes for mitochondrial electron transport and ATP synthesis via the NADH-DH complex were consistently down regulated (Additional file 4: Table S2 Bin 9.1) consistent with a reduced capacity for generation of energy. This may have contributed to a reduced capacity of the transportation system in the root phloem. Substances transported in the phloem include sugar, peptides, amino acids, potassium and metal ions. Notably, both major intrinsic proteins PIP (plasma membrane intrinsic proteins) and TIP (tonoplast intrinsic protein) with water channel activity were down-regulated ( Table 3). The repression of PIP1.1, PIP2.2 and TIP in the phloem is consistent with reduced water movement. Overexpression of ammonium transporters AMT2, AMT1;1 and

Hormones and plant defense responses
Metabolic pathways associated with abscisic acid, auxin, brassinosteroids, cytokinin, ethylene, gibberellin jasmonate and salicylate were all overwhelmingly repressed with induction of only ACCO1 (1-aminocyclopropane-1carboxylate oxidase 1), used to produce ethylene (Additional file 4:  (Table 3). Calciumdependent protein kinase 1 (CDPK1) and calcineurin Blike protein functioning in calcium signaling and phosphoinositide kinase in phospholipid signaling were repressed (Table 3). A few receptors, representing legumelectin, S-locus glycoprotein, wall associated kinase, ATPase E1-E2 and leucine-rich repeat, play roles in signaling were induced (Table 3). Notably, genes encoding Kunitz-type protease inhibitor and endochitinase A precursor were extremely down-regulated. There was also a set of up-regulated genes of unknown function that were not assigned into bins (Table 3).

Correlation between RNA-Seq data and RT-qPCR
The level of expression of genes encoding transcripts annotated as disease resistance/polyprotein/citrus endogenous pararetrovirus was estimated by both RNA-Seq and RT-PCR (Additional file 5 Table S3). The correlation of the two techniques was evaluated by the Spearman's rho value, which 0.70 for all seven libraries (Fig. 4).

Presence of genomic sequences from other organisms in the libraries
We examined the sequence libraries for evidence of other organisms that may be associated with citrus blight disease ( Table 4). As expected, the vast majority of the reads were from citrus, but all of the libraries also had reads from plant viruses, notably citrus tristeza virus. The libraries from the roots of trees with citrus blight also had significant numbers of reads mapped to pararetrovirus, CBaPRV, similar to pararetroviruses that we have described earlier from citrus blight infected trees [42]. We also observed large number of reads that mapped to Pseudomonas spp. Bacteria in the genus Pseudomonas have been observed previously in association with citrus blight [21,22]. Sequences of Xylella fastidiosa were not observed.  The quality of the RNA used to make the transcript libraries was variable when the libraries were prepared from blight affected trees, but not from the healthy control tree ( Table 1). The quality of the RNA library was inversely correlated with the number of reads that mapped to CBaPRV in each of the libraries (R 2 = 0.7114; Fig. 5a). There was no such correlation when the RIN was plotted against the number of reads that were mapped to Pseudomonas in the same libraries (R 2 = 0.0057; Fig. 5b).

Discussion
Many more genes were down-than up-regulated in response to citrus blight, accounting for 89% of the total number of genes whose expression was altered in the seven libraries, in general agreement with an earlier Ave. FC, average fold change of seven libraries NA, no gene symbol study [55]. In contrast with the earlier study, trees sampled in this study had fully developed citrus blight with symptoms of die back, foliar zinc deficiency and severely impaired water uptake. Xylem plugging leading to impaired transport of water is a physical change in the plants with blight. Over production of proline-rich proteins could help citrus plant cells to maintain water potential, thus enhancing tolerance to water stress as observed in transgenic tobacco [25]. Tobacco transformed with trehalose phosphate synthase (TPS1) and trehalose-6-phosphate synthase, exhibited drought tolerance [5] and improved photosynthetic efficiency, respectively under drought [26]. Presumably, the induction of trehalose assists the citrus plant adapt to dehydration caused by blight.
As a consequence of water loss, modifications of protein residues such as deamination or oxidation result in the disruption of native protein conformation or protein activity. The tree responds to this by activating the ubiquitination process and proteases to remove the damaged proteins and recycle the amino acids. Concanavalin A-like lectin protein (Con A-like) is involved in phosphorylation and alterations in protein phosphorylation suggest reversible protein phosphorylation is being used by the orange trees as a regulator in response to water deficit [15]. Thus, up-regulation of Con A-like could favor repair of damage in the plants. OMT1, CYP79A1, CYP79A2 and CYP79A9 proteins function in the biosynthesis of lignin and glucosinolates and their overexpression is likely to be associated with lignin and amorphous plugging in the xylem vessels [34,51].
Aquaporins, PIP (plasma membrane intrinsic proteins) and TIP (tonoplast intrinsic protein), are water channel proteins that facilitate water transport along transmembrane water potential gradients [39] and their downregulation is consistent with reduced water transport in the CB plants. Roots from trees suffering from citrus blight were previously shown to experience nitrogen  starvation, while nitrogen accumulated in the trunk wood, bark, and leaf tissue. This could be consistent with the up-regulation of transporters for ammonia. The repression of MDR1 possibly reflected the impairment of IAA gradient-driven transport in the roots, consistent with the decrease of IAA oxidase activities in declining trees. The repression of transcripts encoding zinc transporters is consistent with the accumulation of zinc in the roots as observed in the trunk phloem [3,4] and would contribute to the foliar symptoms of zinc deficiency since the ability to mobilize zinc is diminished. CB-affected trees also exhibit visible symptoms of wilting and general lack of flush [41,54]. CDPKs and calcineurin B-like proteins effect calcium signaling and gene expression in response to drought [31,50]. Their repression, along with phosphoinositide kinase indicated that calcium and phospholipid signaling pathways were blocked and thus amplification of the water stress signal after its perception was likely impaired. This is consistent with the negative regulation of nucleic acid, carbohydrate, lipid, protein cell organization and development associated pathways as well as a lack of water and demand for new growth, consistent with the earlier report [55]. Our data indicate that phloem and xylem function were heavily disturbed with accompanying unbalanced ion homeostasis under blight stress. It is especially notable that cellulose synthesis in the cell wall was sharply down regulated in roots of trees affected by CB. At the same time, biosynthetic pathways that lead to lignin production were sharply upregulated. In a healthy tree, cellulose fibers and lignin molecules are assembled into wood. It appears that in trees affected by CB the balance between cellulose and lignin production is perturbed leading to an accumulation of soluble lignin in xylem vessels [51], which could lead to the plugging of xylem vessels observed. Although water cannot be forced into the trunks by syringe injection of trees affected by CB, the trees do not die of wilt, but suffer dieback, lack of flush and fail to thrive. The uptake of water in blighted trees has been experimentally determined by gravity injection to be in the range of 5-43 ml/ day compared to 151-503 ml/day among the healthy control trees [12]. Thus, the transpiration stream is severely impaired rather than completely stopped.
Although the trees are very clearly under water stress, the expression of genes related to abscisic acid were not uniformly regulated either up or down, in contrast with results of trees at earlier stages of CB [55]. Expression of ACCO1 synthase was increased consistent with a generalized stress response, also in contrast with what was seen at earlier stages of the CB. The increased activity of the ubiquitin pathway for protein degradation suggests that the roots of trees with citrus blight process a great deal of denatured and damaged proteins which would otherwise support necessary physiological processes. Accelerated degradation of proteins would explain the general up regulation of transporters of oligopeptides and ammonium as the tree attempts to maintain protein homeostasis. Taken together, the data explain why trees with citrus blight are not productive but may linger for years in a 'zombie-like' state.
Transcript libraries can be mined to discover signature sequences to identify unknown organisms such as viruses that are associated with plant diseases of unknown etiology [30]. Sequences from non-citrus genomes were found in the libraries from trees affected by citrus blight. CTV is now ubiquitous in Florida and so is expected in the libraries but is not associated with citrus blight because citrus blight was present in Florida for nearly a century before CTV was introduced to the United States [7,41]. There are large numbers of retrotransposon and pararetrovirus-related sequences in the sweet orange genome [52] and some are located at a locus that confers resistance to CTV [53] and may be expressed in response to stressful conditions such as occur in citrus blight. The sequence read libraries were mapped to the reference genome of the CBaPRV before they were mapped to the citrus genome to account for the presence of these pararetrovirus sequences in the citrus genome. The number of reads of the CBaPRV were elevated in libraries prepared from roots of the blight affected trees (Table 4). Pseudomonads have been reported previously in association with citrus blight [21,22] and are likely to be present as opportunists. Xylella fastidosa was proposed as an etiologic agent for citrus blight [27,28] but was not present in our samples.
We have previously observed pararetrovirus sequences in small RNA libraries such as used for this study [43]. It is interesting that the number of CBaPRV reads present in the libraries was strongly and inversely correlated with the quality of the RNA in the libraries. Such a correlation was not observed for reads that mapped to Pseudomonas in similar numbers. Pseudomonads have been found in association with roots of trees with citrus blight [21,22] but their association with the roots is not likely to be related to citrus blight and they were not correlated with poor RNA quality as was CBaPRV.

Conclusion
Water deficit and xylem plugging are two core physical changes under blight stress and they are primarily responsible for the overwhelming repression of genes and disrupted pathways. Down regulation of cellulose production and up regulation of lignin production creates an imbalance for production of wood. The plants altered regulation of both primary and secondary metabolism, transport and hormone associated pathways but failed to reverse the citrus blight syndrome, which is the result of the co-ordination of physiological and biochemical alterations at the cellular and molecular level.

Plant materials, RNA extraction and sequencing
We used a direct water uptake assay in citrus trunks [33] to identify symptomatic trees from commercial groves located in four different growing regions (Collier, DeSoto, Lake and Polk counties) in Florida. Scaffold root samples (~1.5 cm diameter) from field trees that failed to take up water and healthy root samples from a 4-year old greenhouse-grown tree 'Navel orange' were collected. A greenhouse-grown tree was used as the healthy control because citrus blight does not occur in greenhouse grown trees and the incidence of blight in the groves sampled was such that an apparently healthy tree used as a control may in fact have been in early stage of citrus blight. Roots were placed on ice immediately after harvest and taken to the laboratory. Roots were washed and phloem-enriched tissues (bark) were scraped with a razor blade from roots and stored in RNA-Later (Life Technologies) at −20C. Aliquots of the phloem-enriched samples were ground in liquid nitrogen and RNA was extracted using the Trizol Reagent (Life Technologies). RNA quantity and quality was assessed with a Qubit 2 Fluorimeter (Life Technologies) and a Bioanalyzer 2100 [24]. The eight RNA extracts, seven blight samples and one healthy, were sent out for 2 X 100 bp HiSeq 2500 sequencing (Genewiz, Inc., South Plainfield, NJ).
In previous analyses of root samples from trees affected by citrus blight, 21 nt and 24 nt peaks following RNA-Seq were analyzed [42,43]. Several transcripts annotated as polyprotein/citrus endogenous pararetrovirus and disease resistance were found to have altered expression in response to citrus blight. Consistent expression changes of the tested genes were observed in the seven libraries from CB affected trees.
RT-qPCR assay of transcripts of polyprotein and disease related genes A total of 3 μg of RNA was reverse transcribed for firststrand cDNA syntheses with GoTaq® 2-step RT-qPCR system (Promega, Madison, WI) following the manufacturer's instructions. cDNA was diluted three-fold with 1 X TE and stored at -20°C for use. Gene-specific primers (Additional file 5: Table S3) were designed with an online tool (Integrated DNA Technologies, Coralville, IA) with melting temperatures of 60°C ± 5°C. qPCR reactions were performed with GoTaq® qPCR Master Mix (Promega) in a Bio-Rad CFX96 system. Twenty microliters of reaction mixture were added to each well with three replicated plant samples for each primer pair for each sample. The amplification program was set at 95°C for 3 min, and 40 cycles of 95°C for 10 s, 60°C for 30 s and 72°C for 30 s. Melting curves were analyzed to ensure that a single product was amplified. The GAPC2 gene has been verified to be a stable reference gene [36] and was used as the internal control in our study. The 2 -△△Ct method was applied for relative quantification [35].

Statistical analysis
One RNA sample from each of 8 trees was run in each of 4 lanes as technical replicates in an Illumina 2500 system to produce 32 libraries of paired-end reads (2 × 100 nucleotides). Raw sequence reads were trimmed to remove possible adapter sequences and nucleotides with poor quality (error rate < 0.05). After trimming, sequence reads shorter than 50 nucleotides were discarded (Genewiz). The trimmed reads were then mapped to the coding sequences of the reference genome Citrus sinensis 2.0 [52] with Bowtie2 [32] ( Table 1). The 4 technical replicates were pooled and the expression level of transcripts (log 2 fold-change) were calculated with GFOLD [18], which evaluates the statistical significance of foldchange values against the Poisson distribution rather than the normal distribution. Differentially expressed transcripts (DETs) were normalized as reads per kilobase per million reads (RPKM) with cut-off values, reads ≥200 and | log 2 FC | ≥ 2, P < 0.01. The functions of DETs were enriched by mapping to C. sinensis 2.0 with Mapman [45]. The mapping file was generated by Mercator using PlaBi as the reference database (http://www. plabipd.de/portal/web/guest/home1). Upset diagrams of co-regulated DETs of the seven libraries were created with TBtools [13]. DETs were also classified by alignment with Arabidopsis thaliana (e-value ≤ e-4) using Panther (http://pantherdb.org) [40].

Mapping of transcript reads to the sweet orange and other genomes
The raw sequence reads were pruned and cleaned using BBDuk version 37.66 with adapters.fa from the BBTools software suite [10]. The technical replicates were then pooled. The reads were mapped using Bowtie2 version 2.3.4 with default parameters [32] and 20 threads to reference genomes as follows: The reads were first mapped to All Plant Virus using the DPV web database [1], then to the CBaPCRV (NCBI XXXXXXX), then to the citrus genome (CGF_00317415) then to a large number of strains of Pseudomonas syringae and P. fluorescens such as (NC_004578.1, NC_004633.1, NC_004632.1, NC_ 007005.1, NC_007005.1) and finally to citrus mitochondria (NC_039463.1). After each mapping step, the remaining unmapped reads were used as input for Bow-tie2 (Table 4).