Transcriptome comparison and gene coexpression network analysis provide a systems view of citrus response to ‘CandidatusLiberibacter asiaticus’ infection
© Zheng and Zhao; licensee BioMed Central Ltd. 2013
Received: 11 April 2012
Accepted: 9 January 2013
Published: 16 January 2013
Huanglongbing (HLB) is arguably the most destructive disease for the citrus industry. HLB is caused by infection of the bacterium, Candidatus Liberibacter spp. Several citrus GeneChip studies have revealed thousands of genes that are up- or down-regulated by infection with Ca. Liberibacter asiaticus. However, whether and how these host genes act to protect against HLB remains poorly understood.
As a first step towards a mechanistic view of citrus in response to the HLB bacterial infection, we performed a comparative transcriptome analysis and found that a total of 21 Probesets are commonly up-regulated by the HLB bacterial infection. In addition, a number of genes are likely regulated specifically at early, late or very late stages of the infection. Furthermore, using Pearson correlation coefficient-based gene coexpression analysis, we constructed a citrus HLB response network consisting of 3,507 Probesets and 56,287 interactions. Genes involved in carbohydrate and nitrogen metabolic processes, transport, defense, signaling and hormone response were overrepresented in the HLB response network and the subnetworks for these processes were constructed. Analysis of the defense and hormone response subnetworks indicates that hormone response is interconnected with defense response. In addition, mapping the commonly up-regulated HLB responsive genes into the HLB response network resulted in a core subnetwork where transport plays a key role in the citrus response to the HLB bacterial infection. Moreover, analysis of a phloem protein subnetwork indicates a role for this protein and zinc transporters or zinc-binding proteins in the citrus HLB defense response.
Through integrating transcriptome comparison and gene coexpression network analysis, we have provided for the first time a systems view of citrus in response to the Ca. Liberibacter spp. infection causing HLB.
Orange juice is among the largest beverage industries in the world. Sweet orange is mainly produced in the subtropical areas in the countries of China, US and Brazil and the Mediterranean basin regions. Sweet orange (Citrus sinensis) belongs to the Citrus genus that includes several other species such as tangerine, mandarin and grapefruit. In horticultural practice, citrus is asexually propagated through grafting the scion onto the stock which is grown through the seeds. The scion has been bred for the desired traits of fruit quality while the stock is mostly selected for supporting the optimal scion growth and increased resistance to biotic and abiotic stresses.
Among the major biotic factors which frequently challenge tree growth and fruit development, Huanglongbing (HLB) or called citrus greening is one of the most destructive diseases. HLB was first reported in 1919 in southern China, and very recently it has been reported in almost all major citrus production areas [1–3]. For example, in Florida alone, HLB has caused the loss of several billion dollars since 2005 when HLB was first reported, ranging from 30–100% of loss in fruit production in citrus groves . HLB is caused by infection of the bacterium, Candidatus Liberibacter spp., which is spread to plants via the vector Asian citrus psyllid (Diaphorina citri) or through grafting of a diseased shoot. The HLB bacterium has three species, Ca. Liberibacter asiaticus (Las), Ca. Liberibacter africanus (Laf) and Ca. Liberibacter americanus (Lam). The genome of the Las species was recently published, with a size of approximately 1.23 Mb . It has been generally accepted that, after infection or inoculation, the HLB bacteria migrate through phloem and, by accumulating there, causes the formation of sieve plug [2, 3, 5–7]. Consequently, the transport of nutrients (such as carbohydrates and amino acids) from the source leaves to various sinks (such as young leaves, fruits and roots) are compromised or even blocked in severely infected plants, leading to the alterations in carbohydrate metabolism for metabolic flow and exhibiting such phenotypes as yellow and blotchy mottles on leaves, variegated fruits and poor root growth [6–8]. Because of the huge impact of HLB in the citrus industry, plant pathologists and horticulturists have long sought after the HLB resistance mechanism in citrus.
A recent survey suggests the existence of genetic variations among different citrus species, varieties and stocks . In general, mandarin, sweet orange and grapefruit are relatively more susceptible to the HLB bacterial infection, while sour orange, lemon, lime, and citrange (a hybrid of sweet orange and the stock C. trifoliata) are less susceptible. This raises the possibility that HLB resistance can be achieved through genetic means. Nevertheless, breeding for the HLB resistance through crossing will be a daunting task, given the complex genetic backgrounds, the nature of asexual propagation and the relatively long juvenile period for citrus. Therefore, many researchers have turned their attentions to finding the target genes that are required or critical for the citrus host response to the HLB bacteria. Transcriptome analysis has been used as a straight forward approach to identify the genes whose expression is altered in citrus leaves in response to the HLB inoculation [5, 10–12]. These studies led to the identification of several hundred or thousand genes that are up- or down-regulated by the HLB bacterial infection. The majority of these genes can be grouped into metabolism, transport and response to stimulus. However, these studies varied significantly in terms of study design (including experimental materials and infection stages) and data analysis (for example, different statistical methods and various fold-change cutoffs). Furthermore, there is a lack of comparison of the results from these different experiments. In addition, how these HLB bacterium-regulated genes are connected in a system remains unknown.
To provide a systems view of citrus response to the HLB bacterial infection, we first performed a comparative study of the previously reported transcriptome datasets. Our results show that there are 21 probe sets (representing up to 19 genes) are commonly up-regulated and a number of genes that are specific to early, late or very late stages of inoculation. Furthermore, using the Pearson correlation coefficient (Pcc)-based unweighted gene coexpression analysis, we constructed an HLB response network. This citrus gene coexpression network consists of 3,507 Probesets and 56,857 interactions. We then mapped certain categories of the HLB responsive genes to the HLB response network, resulting in the formation of several important subnetworks including metabolism, transport, signaling, defense response and hormone response. Taken together, through comparative transcriptome analysis and construction of a citrus gene coexpression analysis, we have provided a systems view of citrus response to the Ca. Liberibacter infection causing HLB.
An overview of comparative analysis of HLB transcriptomes
A list of citrus GeneChip studies used in this analysis
Weeks after inoculation
No. of significantly regulated Probesets reported in this study
The citrus GeneChip contains a total of 30,173 Probesets. Because the Affymetrix company has not provided a comprehensive annotation analysis for those Probesets, it is not known how many unique citrus genes are actually represented in the chip. Therefore, we decided to analyze the number of Probesets that are significantly regulated in response to HLB. The data pre-processing was described in Methods. In brief, those Probesets with the calls of present (P) or marginal (M) in at least two chips in each of the four reports were included in our analysis. For the identification of significantly regulated genes, the adjusted LPE approach was used because of its power in analyzing small samples . In our analysis, a two-fold cutoff was used, resulting in various numbers of genes that were either up- or down-regulated in each of the six studies (see Table 1). The HLB regulated genes for each study were listed in Additional file 1. If the genes significantly regulated in at least one study were added together, we found that a total of 3,345 and 3,230 Probesets were up-regulated and down-regulated, respectively. These Probesets are called “HLB responsive genes” in this study (Additional file 1). The percentage of “HLB responsive genes” identified in this comparative analysis (22%) is similar to that of the bacterial pathogen responsive genes in Arabidopsis . This indicates that either the disease response mechanism could be somehow conserved or these four reports have probably identified most of the HLB responsive genes in the citrus genome. Surprisingly, the study-wise comparison showed that the proportion of the significantly regulated genes overlapped in two of the six studies varied dramatically (from 0.05% to 75%) (see Additional file 2).
Comparative studies reveal commonly regulated and stage-specifically regulated genes by HLB
A list of commonly up-regulated HLB responsive genes in various studies
Arabidopsis gene annotation
GA-responsive GAST1 homolog; BR, GA and ABA responsive expression; located in cell wall; GASA1 (GAST1 PROTEIN HOMOLOG 1)
Response to stimulus
NAC096 (Arabidopsis NAC domain containing protein 96); transcription factor
APL3 (sugar-inducible); glucose-1- phosphate adenylyltransferase
Carbohydrate metabolic process
Metal ion binding
LOX2 (LIPOXYGENASE 2); lipoxygenase; response to JA, bacterium, fungus and wounding; JA biosynthetic process
Response to stimulus
WBC11 (WHITE-BROWN COMPLEX HOMOLOG PROTEIN 11); ATPase, coupled to transmembrane movement of substances / fatty acid transporter; response to stress and ABA
Transport; Response to stimulus
Unknown protein (auxin-ethylene interaction, gemivirus induced)
CCT motif family protein (ethylene, gemivirus-induced)
Phospholipase/carboxylesterase family protein
Other metabolic process
Transcript Assignment: DN622378; Human protein tyrosine phosphatase type IVA
Grape RAP2-4-like ethylene transcription factor
UCRPT02_65A04_b Poncirus trifoliata Roots with Iron Deficiency cDNA clone
Transcript Assignment: CX308038; Hypothetical protein
ZIP5; cation/metal ion transmembrane transporter
ZIP5; cation/metal ion transmembrane transporter
O-methyltransferase family 2 protein
Other metabolic process
Glutaredoxin family protein
BGLU11 (BETA GLUCOSIDASE 11); catalytic/ cation binding / hydrolase, hydrolyzing O-glycosyl compounds
Carbohydrate metabolic process
BGLU11 (BETA GLUCOSIDASE 11); catalytic/ cation binding / hydrolase, hydrolyzing O-glycosyl compounds
Carbohydrate metabolic process
Proton-dependent oligopeptide transport (POT) family protein
A list of early stage-specific HLB responsive genes
Arabidopsis gene annotation
ATHSP90.1 (HEAT SHOCK PROTEIN 90.1); ATP binding / unfolded protein binding
PLA2A (PHOSPHOLIPASE A 2A); lipase/ nutrient reservoir
potassium channel tetramerisation domain-containing protein
GCN5-related N-acetyltransferase (GT) family protein
XTR4 (XYLOGLUCAN ENDOTRANSGLYCOSYLASE 4); hydrolase, acting on glycosyl bonds
AGP16 (ARABINOGALACTAN PROTEIN 16)
Immediate-early fungal elicitor family protein
TPS14 (TERPENE SYNTHASE 14); S-lilool synthase
CWLP (CELL WALL-PLASMA MEMBRANE LINKER PROTEIN); lipid binding
FAD-binding domain-containing protein
ERF5 (ETHYLENE RESPONSIVE ELEMENT BINDING FACTOR 5); D binding / transcription activator/ transcription factor
CWLP (CELL WALL-PLASMA MEMBRANE LINKER PROTEIN); lipid binding
GRAM domain-containing protein / ABA-responsive protein-related
ATAF1; transcription activator/ transcription factor
ATAF1; transcription activator/ transcription factor
ATAF1; transcription activator/ transcription factor
GATL10 (Galacturonosyltransferase-like 10); polygalacturote 4-alpha-galacturonosyltransferase/ transferase, transferring hexo
Response to oxidative stress; LOCATED IN: plasma membrane
TT7 (TRANSPARENT TESTA 7); flavonoid 3'-monooxygese/ oxygen binding
ELI3-2 (ELICITOR-ACTIVATED GENE 3–2); aryl-alcohol dehydrogese/ mannitol dehydrogese
ARR9 (RESPONSE REGULATOR 9); transcription regulator/two-component response regulator
ZIP1 (ZINC TRANSPORTER 1 PRECURSOR); zinc ion transmembrane transporter
AthCHIB (Arabidopsis BASIC CHITINASE); chitinase
AtGDU3 (Arabidopsis GLUTAMINE DUMPER 3)
Trypsin and protease inhibitor family protein/Kunitz family protein
In addition, 103 up-regulated and 74 down-regulated Probesets are specific to the “late” stage of Las infection (Additional file 3). Interestingly, these Probesets represent some genes that belong to the categories of metabolism of carbohydrate, nitrogen and lipids, hormone IAA metabolism, response to chemical stimulus, endomembrane systems and extracellular regions. In addition, while several genes involved in cell wall property regulation (such as expansin) are up-regulated, some genes encoding transcription factors (such as MYB15/52/103) and protein kinases are down-regulated. The most striking feature is that only seven Probesets represent the very late stage-specific genes (Additional file 3). These include the genes that are most closely related to Arabidopsis C domain containing protein 71 (a transcription factor), a copper-binding family protein, a trypsin and protease inhibitor family protein/Kunitz family protein, a myosin heavy chain-related protein, two basic chitinase (CHIB) and one unknown protein encoded by At1g42430. The small number of genes belonging to this very late stage-specific category is likely due to the various experimental conditions because only 26 Probesets are commonly up- or down-regulated even in the four studies within the same very late stage of Las infection (Additional file 1). Nevertheless, as this group of genes were identified from four studies specifically at the very late stage compared to only one study for early and late stages respectively, they could be more reliable than groups of early- and late-stage-specific genes.
Construction and characterization of gene coexpression network for citrus response to HLB
We next determined the robustness of our network across each dataset using the cross-validation approach. We randomly left out one dataset and reconstructed the gene co-expression networks using the remaining three datasets. The resulting four networks were then compared to the network based on all four datasets in terms of network connectivity rank of each gene according to the suggestion described elsewhere . There were strong, highly significant connectivity correlations (R = 0.75, 0.83, 0.83, and 0.80 respectively, each associated with a p-value of less than 2.2e-16) between the network based on all four datasets and the ones reconstructed from any combination of the three datasets. This suggests a high degree of preservation of gene co-expression patterns across the networks based on different datasets.
Second, we performed a GO enrichment analysis for the Probesets in the HLB response network. Among 30,173 Probesets, 22,775 (or 75%) have the Arabidopsis gene ID (AtGID) as their closest orthologs or homologs (Additional file 6). Therefore, these Probesets were assigned GO terms based on the most recent Arabidopsis GO assignment. The remaining Probesets were given three general GO terms: “biological process” (GO:0008150), “molecular function” (GO:0003674), and “cellular component” (GO:0005575). GO enrichment analysis using the hypergeometric statistical method with the Hochberg false discovery rate (FDR) adjustment showed that many GO terms were overrespresented in the HLB response network. Among the overrepresented GO terms (Additional file 5), the nodes belonging to the following six categories were color coded in the HLB response network (see Figure 1): carbohydrate metabolic process (213 Probesets), nitrogen and amino acid metabolic process (207), transport (311), defense response (175), hormone response (238) and signaling (218). The nodes for each of these six categories, together with the nodes belonging to some highly overrepresented GO terms such as response to stress, lipid metabolic process, cell wall and membrane part, were listed in Additional file 7. The p-values of the overrepresented GO terms were listed in Additional file 5.
We also performed a GO enrichment analysis for the hub genes. We arbitrarily divided the 2,247 hubs into two categories: minor hubs (3–99 interactions) and major hubs (> = 100 interactions) and their overrepresented GO terms were summarized in Additional file 8. The major hubs have 13 overrepresented GO terms: carbohydrate metabolic process, primary metabolic process, metabolic process, secondary metabolic process, lipid metabolic process, cellular amino acid and derivative metabolic process, cellular process, localization, transport, establishment of localization, regulation of anatomical structure size, regulation of cell size, and regulation of cellular component size. In addition to these 13 GO terms, the minor hubs have 16 additional overrepresented GO terms, such as response to stimulus (including biotic and abiotic), response to stress, regulation of biological quality and signal transduction.
Analysis of the defense and hormone response subnetworks
Given the importance of carbon and nitrogen metabolism, transport, signaling, defense response and hormone response in the citrus response to the HLB bacterial infection and in general plant defense response, the subnetworks for these six categories were constructed by mapping the Probesets belonging to these categories into the HLB response network. The resulting edges (interactions) were listed in Additional file 7.
Analysis of the early stage HLB response subnetwork
Cit.12214.1.S1_s_at represents a transcription factor most closely related to Arabidopsis NAC096. Mapping this Probeset as the seed node to the HLB response network with the second degree neighbors resulted in an NAC096 subnetwork (Figure 8B). Two medium-size hubs were identified in this subnetwork: Cit.10032.1.S1_x_at and Cit.15606.1.S1_at, both of which were up-regulated transcriptionally by the Las infection (Additional file 1). Cit.10032.1.S1_x_at represents a GA-responsive GAST1 homolog and has been reported to be responsive to other hormones such as BR and ABA [37, 38]. Cit.15606.1.S1_at has interactions with 15 Probesets and is closely related to At1g80130 which encodes Arabidopsis tetratricopeptide repeat (TPR)-like superfamily protein and is responsive to oxidative stress . Given that both GA response and oxidative stress response have been implicated an important role in a relatively resistant variety US-897 in response to the Las infection at the very late stage , our preliminary analysis of the NAC096 subnetwork supports that transcriptional control involving hormone response and oxidative stress response might also be important even at the early stage of the HLB bacterial attack.
Subnetwork analysis reveals transport process as a key component in the HLB response core subnetwork
There are 13 Probesets grouped into the category of carbohydrate metabolic process and 11 Probesets that belong to the hormone response category. For the category of carbohydrate metabolic process, Cit.13437.1.S1_s_at forms a larger hub with 11 interactions, and Cit.17155.1.S1_at forms a smaller hub with seven interactions. Cit.13437.1.S1_s_at represents a citrus gene similar to Arabidopsis APL3 encoding a glucose-1-phosphate adenylyltransferase. Cit.17155.1.S1_at represents a gene closely related to BGLU11 (Beta-glucosidase 11) hydrolysis of O-glycosyl compounds. For the hormone response category, Cit.19674.1.S1_s_at forms a larger hub with 15 interactions, and Cit.10032.1.S1_x_at and Cit.25840.1.S1_s_at form smaller hubs with seven and six interactions respectively. As described above, Cit.19674.1.S1_s_at represents a gene closely related to LOX2 encoding a lipoxygenase and exhibiting response to JA. In Arabidopsis, LOX2 has also been shown to be involved in JA biosynthesis in response to wounding [40, 41] and recently in disease development . As described previously, Cit.10032.1.S1_x_at represents a GA-responsive GAST1 homolog and is connected to the NAC096 transcription factor subnetwork in the HLB early response subnetwork (Figures 7 and 8B). Interestingly, Cit.25840.1.S1_s_at represents a gene very similar to Arabidopsis WBC11 (White-brown complex homologous protein 11) which encodes an ATPase coupled to transmembrane movement of substances or fatty acid transporter . This small hub is responsive to ABA and salt stress but is also involved in fatty acid transport, implying a potential role for hormone signaling in the control of transport process.
The remaining two large hubs in the HLB response core subnetwork are formed by Cit.12172.1.S1_s_at and Cit.15630.1.S1_at. Cit.12172.1.S1_s_at represents a putative O-methyltransferase family 2 protein most closely related to the protein encoded by At4g35160. At4g35160 is only annotated as a general GO term “methylation”, and predicted to contain a winged helix-turn-helix transcription repressor DNA-binding domain without any functional implication. This hub includes 31 interactions, and most of the interactions are with the Probesets related to transport process. Cit.15630.1.S1_at represents a gene closest to At4g33040 which encodes a glutaredoxin family protein. It connects to a transportor hub through Cit. 17265.1.S1_at (a glycine-rich protein without any specific functional annotation) and the two hormone response hubs through Cit.17398.1.S1_at (which represents e gene encoding an unknown protein). In Arabidopsis, At4g33040-encoded glutaredoxin family protein is annotated as a protein involved in cell redox homeostasis, which could have a potential role in response to oxidation stress.
Analysis of a phloem protein subnetwork implicates a potential role for zinc transport in the citrus HLB defense response
The transcriptomes in citrus in response to the HLB bacterial infection have been well documented in four previous reports [5, 10–12], but the information regarding the interactions between the differentially expressed genes is lacking. Through the combination of transcriptome comparative study and gene coexpression network analysis, we have provided for the first time a systems view of how the citrus host plant exerts a genome-wide response to the HLB bacterial infection.
First, we have constructed an HLB response network involving 3,507 Probesets with 56,287 interactions. Using the transcriptome datasets and orthology-based or experimentally verified protein-protein interaction datasets, gene-gene interactions or interactomes have been constructed in the model plants including Arabidopsis and rice (for example, [16, 45–52]) and occasionally in non-model plants such as soybean  and barley . However, there was no report on gene-gene interaction networks in citrus prior to our work. We used the Pcc method to construct a gene coexpression network in citrus, with a focus on the HLB response mechanism. The citrus gene coexpression network will be very useful for the citrus researchers to visualize the subnetworks specific for certain biological processes (such as carbohydrate metabolic process or signaling not analyzed in detail in this report), or to search some potential gene-gene interactions for certain genes or a group of genes in the future. The Citrus Gene Interaction Networks (CitGIN) database has been constructed and made available to the research community to query through the Internet (http://xt.cric.cn/cgr/CitGIN.php).
Second, our analysis of the defense subnetwork has shown that many defense hubs and hormone hubs are intertwined or overlapped. Although the roles of hormone and defense response genes have been discussed in the four previous reports [5, 10–12], our network analysis further indicates that hormone response is interconnected to defense response in citrus when challenged by the HLB bacteria. This may lead to the development of integrating hormone and disease response pathways as a potentially more effective genetic means to improve the citrus resistance to HLB.
Third, our comparative studies of transcriptomes have led to the identification of subsets of commonly up-regulated and stage-specific HLB responsive genes. In contrast to those four GeneChip reports where various statistical methods and fold-change cutoffs were used, we used the same procedure for the analysis of all of the transcriptome datasets. Furthermore, by mapping the subset of commonly up-regulated genes into the HLB response network, we have found that the genes belonging to the categories of carbohydrate metabolic process, transport and hormone response are positioned as the large hubs in the HLB response core subnetwork. This indicates that these three processes constitute a core subnetwork for the citrus host response to the HLB bacterial infection. In addition, we propose that transport is a key component in this HLB response core subnetwork.
Fourth, using PP2 gene as an example of applying the HLB response network, our subnetwork analysis provides an intriguing possibility for the zinc transporter or zinc binding proteins to act with PP2 protein in response to the HLB bacterial infection. PP2 proteins belong to a large gene family in higher plants. However, they have not been assigned a specific biological process, and thus their biological function remains unknown. They are predicted to bind carbohydrates and have been implicated a role in the formation of sieve plug or replacement phloem [5, 10, 12, 44]. Some of the PP2 genes from other organisms such as melon, cucumber and Arabidopsis are specifically or preferentially expressed in companion cells but their protein products are translocated in sieve elements . This indicates a role for PP2 proteins not only in intracellular signaling but also in long distance intercellular communication . Recent evidence show that PP2-type proteins might be involved in aphid-mediated virus transmission in cucumber  and overexpression of PP2-A1 in Arabidopsis increased resistance to a phloem-feeding insect , indicating a possibly more active role for PP2 family proteins in defense response in plants. In citrus, the Probeset Cit.35955.1.S1_at, which is closely related to Arabidopsis PP2-B8, was dramatically up-regulated at late stage and very late stages. The most surprising feature of the PP2-B8 subnetwork is that the 20 Probesets, which are the first degree neighbors of Cit.35955.1.S1_at, are interconnected frequently between each other. This indicates that these genes might be regulated by the precise coordination of various signaling pathways through transcription factors, chromatin modification or remodeling proteins or other factors. Furthermore, seven of the 20 interacting Probesets encode proteins involved in transport, consistent with our proposal that transport is a key component in the HLB response core subnetwork. In addition, three of the seven transporters are predicted to transport zinc, and the PP2 subnetwork also contains four Probesets which represent the genes encoding zinc-binding proteins. Intriguingly, HLB disease symptom was initially thought to be related to zinc deficiency and the zinc transport system is required for virulence in other organisms [1, 2], and therefore the PP2 subnetwork analysis indicates that zinc transporters or zinc binding proteins may have a potentially important role for citrus to respond to the HLB bacterial infection. Taken together, our analysis using the HLB response network can lead to an intriguing but testable hypothesis (which cannot be done without the network analysis) regarding the role of PP2 proteins and zinc transport system or zinc binding proteins in citrus HLB defense response.
It should be noted that there are some potential limitations in our network study. The first one is GO enrichment analysis. The agriGO web tool, which is based on the hypergeometric method and used in this work, does not take into account the local dependency of GO terms. Using the four algorithms provided in the topGO R package which are proposed to eliminate local dependencies, we have found that four of the six hormone GO terms determined to be overrepresented by agriGO are also overrepresented, while the two other hormones (ABA and ethylene) have their child GO terms being truly overrepresented. Therefore, different algorithms or statistical methods in GO enrichment analysis will probably lead to some differences in terms of the overrepresented GO terms for the nodes in the HLB response network.
The second limitation is due to the small sample size. Computational prediction of gene-gene interactions usually requires large sample size; however relatively small number of samples were recently used to construct gene coexpression networks specific to certain aspects of biology (such as [50, 52, 57]). In our analysis, we used the transcriptome datasets described in four previous reports [5, 10–12]. Among these, only one study is available for early stage and late stage respectively, while there are four studies for very late stage (Table 1). We have found that only a small number (222) of significantly regulated Probesets can be identified for early stage, while almost 600 and 2,000–4,000 differentially expressed Probesets can be found for late and very late stages respectively. The variation in the number of differentially expressed genes at different stages could be caused by the difference in experimental conditions given that different ages and varieties of trees and different sources of inoculants were used in different years in those four reports. However, this variation might lead to some sort of bias towards the very late stage genes. To minimize the possibility that the interactions we have detected were the result of random events due to the small sample size, we have selected a high Pcc cutoff value which has led us to believe that the interactions are more likely statistically significant rather than by random and that the topology of the HLB response network is quite similar to most biological networks. Furthermore, the cross-validation result shows a high degree of preservation of gene coexpression patterns, suggesting that the HLB response network is at least moderately robust and biologically relevant. Therefore, despite some limitations due to the small sample size and the experimental variations, the network reported here should be quite useful for the citrus research community and have provided some novel insights into the citrus HLB defense mechanisms. When larger scale transcriptome datasets become available in the future, similar network analysis will provide a comprehensive picture of the gene networks in citrus.
The most daunting challenge in the citrus postgenomic era remains how to identify the best candidate genes for functional dissection of the HLB response mechanism and for genetic modification with an ultimate goal of improving the HLB resistance in citrus. Genetic variations of HLB susceptibility  clearly shows the potential towards dissection of genetic mechanisms of HLB resistance, but understanding the inheritance patterns and subsequently cloning the disease genes requires a long term effort because of long juvenile phase and complex reproductive biology for citrus. Recent developments have shed some lights into the identification of key hub genes as candidate regulatory genes. For example, a seed germination study found that 22–50% of the Arabidopsis hub genes identified from the seed germination network actually have physiological functions in the control of seed germination . Therefore, the hub genes identified in this report may potentially be the first batch of candidates for the functional test in HLB resistance in citrus.
Through integration of transcriptome comparison and gene coexpression network analysis, we have provided novel insights into the mechanism by which citrus host plants respond to the HLB bacterial infection. Specifically, several biological processes are important in the citrus HLB response network, including carbohydrate metabolic process, nitrogen and amino acid metabolic process, transport, defense response, signaling and hormone response. Furthermore, our results have led us to propose that transport is a key component in the HLB response core subnetwork. This systems view of citrus response to the Ca. Liberibacter spp. infection will be a critical first step towards dissecting the genetic mechanisms of HLB response and ultimately improving HLB resistance in citrus.
Data collection and preprocessing
Raw data for citrus Affymetrix GeneChip analysis published by Fan et al.  (GSE29633) and Albrecht and Bowman  (GSE30502) were downloaded from NCBI. Raw data published in  (GSE33459) and  were kindly provided by Drs. Bowman and Wang, respectively. These .cel files were read into R and preprocessed using “rma” function  and normalized using the “normalize.quantiles.robust” function. After quantile normalization, Probesets with an absent (A) call were removed using the “pma” function. Probesets with the calls of present (P) or marginal (M) in at least two samples in each of the four reports above were included in the analysis. All of the statistical analysis and gene expression network construction were performed in the R environment.
Analysis of significantly regulated genes
The adjusted local pooled error (LPE) method  was used to identify differentially expressed transcripts, as this method has been shown to provide high power in analyzing microarray data with small sample size. A gene was called statistically significant if its permutation-based false discovery rate (FDR) p-value was smaller than 0.05 and at least a two-fold change was observed.
Network construction and visualization
For computational reasons, up to 10,000 of the Probesets with highest expression levels were selected from each of the datasets described in the four reports [5, 10–12]. The HLB responsive genes identified in this study (described above) were then added to this list and duplicated ones were removed, resulting in a total of 10,668 common Probesets for each of the four datasets. Gene coexpression network was constructed from the preprocessed files using R package “weighted correlation network analysis” . Following the protocol for constructing gene co-expression network using multiple datasets , we first calculated Pearson correlation matrix for each dataset. We then obtained an overall weighted correlation matrix based on the number of samples used in that dataset. The weight for each correlation matrix was defined as , where n i was the number of samples for ith dataset, n max was the maximum number of samples in all datasets, and s was the number of datasets used. Two nodes were determined to be connected if the absolute value of the Pearson correlation coefficient (Pcc) exceeded 0.93. The threshold of 0.93 was selected such that it gave the best overall fit to each dataset based on the criteria such as the scale-free topology model fitting index, mean network connectivity, and network density. Cytoscape v2.8.2  was used to visualize the networks and Photoshop was used to edit the images.
GO analysis and Arabidopsis orthology prediction
Because of the lack of citrus genome annotation for the Probesets in the Affymetrix chip, the Probesets were used for all analysis. They were annotated using Arabidopsis orthologs or homologs. The Probesets were annotated by searching against the Arabidopsis genome using the tool provided in HarvEST database (http://harvest-web.org). GO terms were assigned to the citrus Probesets based on their corresponding Arabidopsis gene ID (AtGID). For those without AtGID, general GO terms were assigned: “biological process” (GO:0008150), “molecular function” (GO:0003674), and “cellular component” (GO:0005575). GO enrichment analysis was performed using the hypergeometric statistical method with Hochberg FDR adjustment in the AgriCO website (http://bioinfo.cau.edu.cn/agriGO/analysis.php) as described elsewhere .
Candidatus Liberibacter asiaticus
Pearson correlation coefficient
This work was supported in part by a subproject from the Chinese Ministry of Agriculture “Program 948” (Grant No. 2011-G21(2)) and a project from Chongqing Science and Technology Commission (Grant No. cstc2012gg-yyjsB80004). We are grateful to Dr. Kim Bowman (USDA Horticultural Research Center) and Dr. Nian Wang (University of Florida) for providing their published citrus GeneChip data. We greatly appreciate Dr. Matt Geisler (Southern Illinois University) for his help and advice in the network analysis and visualization. We also thank two reviewers for their constructive comments and suggestions.
- Zhao XY: Citrus yellow shoot disease (huanglongbing) in China-a. Review. Proc Int Citriculture. 1981, 1: 466-469.
- de Graça JV: Citrus greening disease. Annu Rev Phytopathol. 2005, 29: 109-136.View Article
- Wang N, Li W, Irey M, Albrigo G, Bo K, Kim JS: Citrus huanglongbing. Tree Forestry Science Biotechnology. 2009, 3: 67-72.
- Duan Y, Zhou L, Hall DG, Li W, Doddapaneni H, Lin H, Liu L, Vahling CM, Gabriel DW, Williams KP, et al: Complete genome sequence of citrus huanglongbing bacterium, ‘Candidatus Liberibacter asiaticus’ obtained through metagenomics. Mol Plant Microbe Interact. 2009, 22 (8): 1011-1020. 10.1094/MPMI-22-8-1011.View ArticlePubMed
- Kim JS, Sagaram US, Burns JK, Li JL, Wang N: Response of sweet orange (Citrus sinensis) to ‘Candidatus Liberibacter asiaticus’ infection: microscopy and microarray analyses. Phytopathology. 2009, 99 (1): 50-57. 10.1094/PHYTO-99-1-0050.View ArticlePubMed
- Achor DS, Etxeberria E, Wang N, Folimonova SY, Chung KR, Albrigo LG: Sequence of anatomical symptom observations in Citrus affected with Huanglongbing disease. Plant Pathology Journal. 2010, 9: 56-64.View Article
- Koh EJ, Zhou L, Williams DS, Park J, Ding N, Duan YP, Kang BH: Callose deposition in the phloem plasmodesmata and inhibition of phloem transport in citrus leaves infected with “Candidatus Liberibacter asiaticus”. Protoplasma. 2012, 249: 687-697. 10.1007/s00709-011-0312-3.View ArticlePubMed
- Fan J, Chen C, Brlansky RH, Gmitter FG, Li ZG: Changes in carbohydrate metabolism in Citrus sinensis infected with ‘Candidatus Liberibacter asiaticus’. Plant Pathol. 2010, 59: 1037-1043. 10.1111/j.1365-3059.2010.02328.x.View Article
- Folimonova SY, Robertson CJ, Garnsey SM, Gowda S, Dawson WO: Examination of the responses of different genotypes of citrus to huanglongbing (citrus greening) under different conditions. Phytopathology. 2009, 99 (12): 1346-1354. 10.1094/PHYTO-99-12-1346.View ArticlePubMed
- Albrecht U, Bowman KD: Gene expression in Citrus sinensis (L.) Osbeck following infection with the bacterial pathogen Candidatus Liberibacter asiaticus causing Huanglongbing in Florida. Plant Sci. 2008, 175: 291-306. 10.1016/j.plantsci.2008.05.001.View Article
- Fan J, Chen C, Yu Q, Brlansky RH, Li ZG, Gmitter FG: Comparative iTRAQ proteome and transcriptome analyses of sweet orange infected by “Candidatus Liberibacter asiaticus”. Physiol Plant. 2011, 143 (3): 235-245. 10.1111/j.1399-3054.2011.01502.x.View ArticlePubMed
- Albrecht U, Bowman KD: Transcriptional response of susceptible and tolerant citrus to infection with Candidatus Liberibacter asiaticus. Plant Sci. 2011, 185–186: 118-130.PubMed
- Murie C, Nadon R: A correction for estimating error when using the Local Pooled Error Statistical Test. Bioinformatics. 2008, 24 (15): 1735-1736. 10.1093/bioinformatics/btn211.View ArticlePubMed
- Tao Y, Xie Z, Chen W, Glazebrook J, Chang HS, Han B, Zhu T, Zou G, Katagiri F: Quantitative nature of Arabidopsis responses during compatible and incompatible interactions with the bacterial pathogen Pseudomonas syringae. Plant Cell. 2003, 15 (2): 317-330. 10.1105/tpc.007591.PubMed CentralView ArticlePubMed
- Miller JA, Horvath S, Geschwind DH: Divergence of human and mouse brain transcriptome highlights Alzheimer disease pathways. Proc Natl Acad Sci USA. 2010, 107 (28): 12698-12703. 10.1073/pnas.0914257107.PubMed CentralView ArticlePubMed
- Geisler-Lee J, O’Toole N, Ammar R, Provart NJ, Millar AH, Geisler M: A predicted interactome for Arabidopsis. Plant Physiol. 2007, 145 (2): 317-329. 10.1104/pp.107.103465.PubMed CentralView ArticlePubMed
- Yanhui C, Xiaoyuan Y, Kun H, Meihua L, Jigang L, Zhaofeng G, Zhiqiang L, Yunfei Z, Xiaoxiao W, Xiaoming Q, et al: The MYB transcription factor superfamily of Arabidopsis: expression analysis and phylogenetic comparison with the rice MYB family. Plant Mol Biol. 2006, 60 (1): 107-124. 10.1007/s11103-005-2910-y.View ArticlePubMed
- Nurmberg PL, Knox KA, Yun BW, Morris PC, Shafiei R, Hudson A, Loake GJ: The developmental selector AS1 is an evolutionarily conserved regulator of the plant immune response. Proc Natl Acad Sci USA. 2007, 104 (47): 18795-18800. 10.1073/pnas.0705586104.PubMed CentralView ArticlePubMed
- Yang JY, Iwasaki M, Machida C, Machida Y, Zhou X, Chua NH: betaC1, the pathogenicity factor of TYLCCNV, interacts with AS1 to alter leaf development and suppress selective jasmonic acid responses. Genes Dev. 2008, 22 (18): 2564-2577. 10.1101/gad.1682208.PubMed CentralView ArticlePubMed
- Dohmann EM, Kuhnle C, Schwechheimer C: Loss of the CONSTITUTIVE PHOTOMORPHOGENIC9 signalosome subunit 5 is sufficient to cause the cop/det/fus mutant phenotype in Arabidopsis. Plant Cell. 2005, 17 (7): 1967-1978. 10.1105/tpc.105.032870.PubMed CentralView ArticlePubMed
- Mukhtar MS, Carvunis AR, Dreze M, Epple P, Steinbrenner J, Moore J, Tasan M, Galli M, Hao T, Nishimura MT, et al: Independently evolved virulence effectors converge onto hubs in a plant immune system network. Science. 2011, 333 (6042): 596-601. 10.1126/science.1203659.PubMed CentralView ArticlePubMed
- Du Z, Zhou X, Ling Y, Zhang Z, Su Z: agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 2010, 38: W64-W70. 10.1093/nar/gkq310.PubMed CentralView ArticlePubMed
- Alexa A, Rahnenfuhrer J, Lengauer T: Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 2006, 22 (13): 1600-1607. 10.1093/bioinformatics/btl140.View ArticlePubMed
- Alexa A, Rahnenführer J: Gene set enrichment analysis with topGO. http://bioconductororg/packages/28/bioc/html/topGOhtml 2012
- Vlot AC, Dempsey DA, Klessig DF: Salicylic Acid, a multifaceted hormone to combat disease. Annu Rev Phytopathol. 2009, 47: 177-206. 10.1146/annurev.phyto.050908.135202.View ArticlePubMed
- Zhang X, Francis MI, Dawson WO, Graham JH, Orbović V, Triplett EW, Mou Z: Over-expression of the Arabidopsis NPR1 gene in citrus increases resistance to citrus canker. Eur J Plant Pathol. 2010, 128: 91-100. 10.1007/s10658-010-9633-x.View Article
- Vailleau F, Daniel X, Tronchet M, Montillet JL, Triantaphylides C, Roby D: A R2R3-MYB gene, AtMYB30, acts as a positive regulator of the hypersensitive cell death program in plants in response to pathogen attack. Proc Natl Acad Sci USA. 2002, 99 (15): 10179-10184. 10.1073/pnas.152047199.PubMed CentralView ArticlePubMed
- Pandey SP, Roccaro M, Schon M, Logemann E, Somssich IE: Transcriptional reprogramming regulated by WRKY18 and WRKY40 facilitates powdery mildew infection of Arabidopsis. Plant J. 2010, 64 (6): 912-923. 10.1111/j.1365-313X.2010.04387.x.View ArticlePubMed
- Zeng LR, Vega-Sanchez ME, Zhu T, Wang GL: Ubiquitination-mediated protein degradation and modification: an emerging theme in plant-microbe interactions. Cell Res. 2006, 16 (5): 413-426. 10.1038/sj.cr.7310053.View ArticlePubMed
- Hua Z, Vierstra RD: The cullin-RING ubiquitin-protein ligases. Annu Rev Plant Biol. 2011, 62: 299-334. 10.1146/annurev-arplant-042809-112256.View ArticlePubMed
- Lee DH, Choi HW, Hwang BK: The pepper E3 ubiquitin ligase RING1 gene, CaRING1, is required for cell death and the salicylic acid-dependent defense response. Plant Physiol. 2011, 156 (4): 2011-2025. 10.1104/pp.111.177568.PubMed CentralView ArticlePubMed
- van Verk MC, Bol JF, Linthorst HJ: WRKY transcription factors involved in activation of SA biosynthesis genes. BMC Plant Biol. 2011, 11: 89-10.1186/1471-2229-11-89.PubMed CentralView ArticlePubMed
- Falk A, Feys BJ, Frost LN, Jones JD, Daniels MJ, Parker JE: EDS1, an essential component of R gene-mediated disease resistance in Arabidopsis has homology to eukaryotic lipases. Proc Natl Acad Sci USA. 1999, 96 (6): 3292-3297. 10.1073/pnas.96.6.3292.PubMed CentralView ArticlePubMed
- Zhu S, Jeong RD, Venugopal SC, Lapchyk L, Navarre D, Kachroo A, Kachroo P: SAG101 forms a ternary complex with EDS1 and PAD4 and is required for resistance signaling against turnip crinkle virus. PLoS Pathog. 2011, 7 (11): e1002318-10.1371/journal.ppat.1002318.PubMed CentralView ArticlePubMed
- Peart JR, Cook G, Feys BJ, Parker JE, Baulcombe DC: An EDS1 orthologue is required for N-mediated resistance against tobacco mosaic virus. Plant J. 2002, 29 (5): 569-579. 10.1046/j.1365-313X.2002.029005569.x.View ArticlePubMed
- Hu G, DeHart AK, Li Y, Ustach C, Handley V, Navarre R, Hwang CF, Aegerter BJ, Williamson VM, Baker B: EDS1 in tomato is required for resistance mediated by TIR-class R genes and the receptor-like R gene Ve. Plant J. 2005, 42 (3): 376-391. 10.1111/j.1365-313X.2005.02380.x.View ArticlePubMed
- Raventos D, Meier C, Mattsson O, Jensen AB, Mundy J: Fusion genetic analysis of gibberellin signaling mutants. Plant J. 2000, 22 (5): 427-438. 10.1046/j.1365-313X.2000.00759.x.View ArticlePubMed
- Bouquin T, Meier C, Foster R, Nielsen ME, Mundy J: Control of specific gene expression by gibberellin and brassinosteroid. Plant Physiol. 2001, 127 (2): 450-458. 10.1104/pp.010173.PubMed CentralView ArticlePubMed
- Luhua S, Ciftci-Yilmaz S, Harper J, Cushman J, Mittler R: Enhanced tolerance to oxidative stress in transgenic Arabidopsis plants expressing proteins of unknown function. Plant Physiol. 2008, 148 (1): 280-292. 10.1104/pp.108.124875.PubMed CentralView ArticlePubMed
- Bell E, Mullet JE: Characterization of an Arabidopsis lipoxygenase gene responsive to methyl jasmonate and wounding. Plant Physiol. 1993, 103 (4): 1133-1137. 10.1104/pp.103.4.1133.PubMed CentralView ArticlePubMed
- Bell E, Creelman RA, Mullet JE: A chloroplast lipoxygenase is required for wound-induced jasmonic acid accumulation in Arabidopsis. Proc Natl Acad Sci USA. 1995, 92 (19): 8675-8679. 10.1073/pnas.92.19.8675.PubMed CentralView ArticlePubMed
- Thatcher LF, Manners JM, Kazan K: Fusarium oxysporum hijacks COI1-mediated jasmonate signaling to promote disease development in Arabidopsis. Plant J. 2009, 58 (6): 927-939. 10.1111/j.1365-313X.2009.03831.x.View ArticlePubMed
- Panikashvili D, Savaldi-Goldstein S, Mandel T, Yifhar T, Franke RB, Hofer R, Schreiber L, Chory J, Aharoni A: The Arabidopsis DESPERADO/AtWBC11 transporter is required for cutin and wax secretion. Plant Physiol. 2007, 145 (4): 1345-1360. 10.1104/pp.107.105676.PubMed CentralView ArticlePubMed
- Dinant S, Clark AM, Zhu Y, Vilaine F, Palauqui JC, Kusiak C, Thompson GA: Diversity of the superfamily of phloem lectins (phloem protein 2) in angiosperms. Plant Physiol. 2003, 131 (1): 114-128. 10.1104/pp.013086.PubMed CentralView ArticlePubMed
- Ma S, Gong Q, Bohnert HJ: An Arabidopsis gene network based on the graphical Gaussian model. Genome Res. 2007, 17 (11): 1614-1625. 10.1101/gr.6911207.PubMed CentralView ArticlePubMed
- De Bodt S, Proost S, Vandepoele K, Rouze P, Van de Peer Y: Predicting protein-protein interactions in Arabidopsis thaliana through integration of orthology, gene ontology and co-expression. BMC Genomics. 2009, 10: 288-10.1186/1471-2164-10-288.PubMed CentralView ArticlePubMed
- Mao L, Van Hemert JL, Dash S, Dickerson JA: Arabidopsis gene co-expression network and its functional modules. BMC Bioinforma. 2009, 10: 346-10.1186/1471-2105-10-346.View Article
- Lin M, Zhou X, Shen X, Mao C, Chen X: The predicted Arabidopsis interactome resource and network topology-based systems biology analyses. Plant Cell. 2011, 23 (3): 911-922. 10.1105/tpc.110.082529.PubMed CentralView ArticlePubMed
- Bassel GW, Lan H, Glaab E, Gibbs DJ, Gerjets T, Krasnogor N, Bonner AJ, Holdsworth MJ, Provart NJ: Genome-wide network model capturing seed germination reveals coordinated regulation of plant cellular phase transitions. Proc Natl Acad Sci USA. 2011, 108 (23): 9709-9714. 10.1073/pnas.1100958108.PubMed CentralView ArticlePubMed
- Peng FY, Weselake RJ: Gene coexpression clusters and putative regulatory elements underlying seed storage reserve accumulation in Arabidopsis. BMC Genomics. 2011, 12: 286-10.1186/1471-2164-12-286.PubMed CentralView ArticlePubMed
- Hamada K, Hongo K, Suwabe K, Shimizu A, Nagayama T, Abe R, Kikuchi S, Yamamoto N, Fujii T, Yokoyama K, et al: OryzaExpress: an integrated database of gene expression networks and omics annotations in rice. Plant Cell Physiol. 2011, 52 (2): 220-229. 10.1093/pcp/pcq195.PubMed CentralView ArticlePubMed
- Aya K, Suzuki G, Suwabe K, Hobo T, Takahashi H, Shiono K, Yano K, Tsutsumi N, Nakazono M, Nagamura Y, et al: Comprehensive network analysis of anther-expressed genes in rice by the combination of 33 laser microdissection and 143 spatiotemporal microarrays. PLoS One. 2011, 6 (10): e26162-10.1371/journal.pone.0026162.PubMed CentralView ArticlePubMed
- Afzal AJ, Natarajan A, Saini N, Iqbal MJ, Geisler M, El Shemy HA, Mungur R, Willmitzer L, Lightfoot DA: The nematode resistance allele at the rhg1 locus alters the proteome and primary metabolism of soybean roots. Plant Physiol. 2009, 151 (3): 1264-1280. 10.1104/pp.109.138149.PubMed CentralView ArticlePubMed
- Mochida K, Uehara-Yamaguchi Y, Yoshida T, Sakurai T, Shinozaki K: Global landscape of a co-expressed gene network in barley and its application to gene discovery in Triticeae crops. Plant Cell Physiol. 2011, 52 (5): 785-803. 10.1093/pcp/pcr035.PubMed CentralView ArticlePubMed
- Bencharki B, Boissinot S, Revollon S, Ziegler-Graff V, Erdinger M, Wiss L, Dinant S, Renard D, Beuve M, Lemaitre-Guillier C, et al: Phloem protein partners of Cucurbit aphid borne yellows virus: possible involvement of phloem proteins in virus transmission by aphids. Mol Plant Microbe Interact. 2010, 23 (6): 799-810. 10.1094/MPMI-23-6-0799.View ArticlePubMed
- Zhang C, Shi H, Chen L, Wang X, Lu B, Zhang S, Liang Y, Liu R, Qian J, Sun W, et al: Harpin-induced expression and transgenic overexpression of the phloem protein gene AtPP2-A1 in Arabidopsis repress phloem feeding of the green peach aphid Myzus persicae. BMC Plant Biol. 2011, 11: 11-10.1186/1471-2229-11-11.PubMed CentralView ArticlePubMed
- Krouk G, Mirowski P, LeCun Y, Shasha DE, Coruzzi GM: Predictive network modeling of the high-resolution dynamic plant transcriptome in response to nitrate. Genome Biol. 2010, 11 (12): R123-10.1186/gb-2010-11-12-r123.PubMed CentralView ArticlePubMed
- Bolstad BM, Irizarry RA, Astrand M, Speed TP: A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003, 19 (2): 185-193. 10.1093/bioinformatics/19.2.185.View ArticlePubMed
- Langfelder P, Horvath S: WGCNA: an R package for weighted correlation network analysis. BMC Bioinforma. 2008, 9: 559-10.1186/1471-2105-9-559.View Article
- Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T: Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics. 2011, 27 (3): 431-432. 10.1093/bioinformatics/btq675.PubMed CentralView ArticlePubMed
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.