- Research article
- Open Access
Large-scale identification of wheat genes resistant to cereal cyst nematode Heterodera avenae using comparative transcriptomic analysis
BMC Genomics volume 16, Article number: 801 (2015)
Cereal cyst nematode Heterodera avenae, an important soil-borne pathogen in wheat, causes numerous annual yield losses worldwide, and use of resistant cultivars is the best strategy for control. However, target genes are not readily available for breeding resistant cultivars. Therefore, comparative transcriptomic analyses were performed to identify more applicable resistance genes for cultivar breeding.
The developing nematodes within roots were stained with acid fuchsin solution. Transcriptome assemblies and redundancy filteration were obtained by Trinity, TGI Clustering Tool and BLASTN, respectively. Gene Ontology annotation was yielded by Blast2GO program, and metabolic pathways of transcripts were analyzed by Path_finder. The ROS levels were determined by luminol-chemiluminescence assay. The transcriptional gene expression profiles were obtained by quantitative RT-PCR.
The RNA-sequencing was performed using an incompatible wheat cultivar VP1620 and a compatible control cultivar WEN19 infected with H. avenae at 24 h, 3 d and 8 d. Infection assays showed that VP1620 failed to block penetration of H. avenae but disturbed the transition of developmental stages, leading to a significant reduction in cyst formation. Two types of expression profiles were established to predict candidate resistance genes after developing a novel strategy to generate clean RNA-seq data by removing the transcripts of H. avenae within the raw data before assembly. Using the uncoordinated expression profiles with transcript abundance as a standard, 424 candidate resistance genes were identified, including 302 overlapping genes and 122 VP1620-specific genes. Genes with similar expression patterns were further classified according to the scales of changed transcript abundances, and 182 genes were rescued as supplementary candidate resistance genes. Functional characterizations revealed that diverse defense-related pathways were responsible for wheat resistance against H. avenae. Moreover, phospholipase was involved in many defense-related pathways and localized in the connection position. Furthermore, strong bursts of reactive oxygen species (ROS) within VP1620 roots infected with H. avenae were induced at 24 h and 3 d, and eight ROS-producing genes were significantly upregulated, including three class III peroxidase and five lipoxygenase genes.
Large-scale identification of wheat resistance genes were processed by comparative transcriptomic analysis. Functional characterization showed that phospholipases associated with ROS production played vital roles in early defense responses to H. avenae via involvement in diverse defense-related pathways as a hub switch. This study is the first to investigate the early defense responses of wheat against H. avenae, not only provides applicable candidate resistance genes for breeding novel wheat cultivars, but also enables a better understanding of the defense mechanisms of wheat against H. avenae.
Cereal cyst nematode (CCN) Heterodera avenae is an important soil-borne wheat pathogen found worldwide that causes significant annual yield losses [1, 2]. Planting and breeding of resistant wheat cultivars are considered as the most economical, environmentally sustainable, and accessible strategies for control of the disease caused by CCN . In past decades, approximately 13 genes with underlying resistance to CCN were identified based on traditional mapping-based cloning strategy, including nine Cre genes (1–8, R) and four Ha (1–4) genes. Few sources of genetic resistance are documented against CCN in the hexaploid bread wheat Triticum aestivum. The dominant resistance gene Cre1 was characterized in the line Aus 10894/Loros and was extensively used as a source of resistance in commercial breeding programs [4, 5]. Another resistance gene, i.e., Cre8, was also derived directly from bread wheat T. aestivum [6, 7]. Sources of CCN resistance have been found predominately in wild grasses and relatives of bread wheat, i.e., Aegilops ventricosa, Ae. squarrosa, and Ae. triuncialis, as well as in cultivated rye, i.e., Secale cereal [4, 5, 8, 9]. The resistance gene CreR was identified from the long arm of 6R within rye , whereas the other six Cre resistance genes originated from Aegilops spp. . In addition, the four resistance genes Ha1-4 were derived from barley [5, 10]. Moreover, most of these resistance genes conferred resistance to limited or specific pathotypes of CCN [5, 10]. The resistant introgression line H93-8 was generated by transferring the resistance gene Cre2 from wild grass Ae.ventricosa to hexaploid bread wheat and exhibited high resistance to CCN pathotypes . However, applicable resistance genes against H. avenae are still limited, and thus, identification of additional genes that confer resistance to CCN are urgently needed.
Isolation and investigation of only a few limited number of resistance genes is time-consuming via the traditional strategy. Recently, transcriptomics has been developed for high-throughput identification of candidate genes of interest and is used to quickly and globally isolate candidate genes against important plant parasitic nematodes, i.e., the soybean cyst nematode Heterodera glycines and root knot nematode Meloidogyne spp. [11–16]. Many defense-related genes and pathways were identified in incompatible soybean lines during early defense responses to H. glycines via comparative transcriptomic analysis [13, 14]. During the resistant reaction of soybean Glycine max genotype PI 548402 (Peking) to H. glycines, lipoxygenase-9 and lipoxygenase-4 were the most highly induced components. Components of the phenylpropanoid pathway (PAL) (i.e., phenylalanine ammonia lyase, isoflavone reductase, chalcone isomerase, cinnamoyl-CoA reductase and caffeic acid O-methyltransferase) were also remarkably upregulated, which indicated the importance of the jasmonic acid and phenylpropanoid signaling pathways during the early defense response of resistant soybean to H. glycines . A microarray analysis was performed using two genetically related soybean sister lines TN02-226 (resistant) and TN02-275 (susceptible) with the race 2 population of H. glycines, and the results showed that 42 transcripts increased in the resistant line but decreased in the susceptible line. These genes were associated with metabolism, cell-wall modification, signal transduction, transcription, and defense . During the early defense response of the resistant reaction of G. max genotype PI 88788 to H. glycines population NL1-RHg/HG-type 7, the jasmonic acid biosynthesis and 13-lipoxygenase pathways were significantly induced in transcript abundance, similar to results from the previous study . The most highly induced pathways contained components of jasmonic acid biosynthesis, flavonoid biosynthesis, 13-lipoxygenase pathway, and phenylpropanoid biosynthesis, which demonstrated that the jasmonic acid defense pathway was involved in the localized resistant reaction of G. max [PI 88788] to H. glycines . Similarly, global gene expression changes in roots during incompatible and compatible associations with M. incognita were also investigated using comparative transcriptomic analysis, and 217 genes were significantly differentially expressed during the time of M. incognita infection corresponding to establishment of feeding sites, and 58 genes exhibited differential regulation in resistant roots compared with uninfected roots, including the glycosyltransferase. Functional characterization demonstrated that the glycosyltransferase gene was required for Mi-mediated nematode resistance . Until now, no transcriptomic analysis of incompatible defense response to CCN has been reported, and only one report exists on the compatible response of Ae. variabilis against CCN .
In this study, to apply high-throughput identification of additional candidate genes associated with underlying resistance against CCN, we performed comparative transcriptomic analysis using an incompatible wheat cultivar VP1620 and a compatible wheat cultivar WEN19 inoculated with CCN to achieve the corresponding incompatible (_I) and compatible (_C) reactions at early three time points of 24 h, 3 d, and 8 d. The VP1620 used herein allowed CCN to penetrate and migrate normally within its roots but significantly disturbed the transition of developmental stages, eventually resulting in a dramatic reduction in cyst formation. A novel methodology was generated to establish two types of expression profiles and narrow the scopes of candidate resistance genes. A total of 606 candidate resistance genes were identified, including 484 overlapping genes as well as 122 VP1620-specific genes, and functional characterization analysis indicated that phospholipase was likely to play vital roles against CCN as a hub switch involved in many defense-related pathways.
Results and discussion
Resistant wheat VP1620 decreased cyst formation of H. avenae
To globally identify wheat resistance genes against H. avenae, an incompatible wheat cultivar VP1620 was used to perform comparative transcriptomic analyses. The number of cysts formed on VP1620 roots was significantly reduced by approximately 8-fold compared to the compatible host WEN19 (Fig. 1a). In this study, the early response of VP1620 against H. avenae was characterized by choosing three key time points at 24 h, 3 d, and 8 d, and two main reasons for these choices were taken into consideration. First, these three time points were typical time points that represent the penetration stage (24 h), migration stage (3 d) and feeding site establishment and maintenance (8 d) during the early reaction of parasitism and pathogenesis of H. avenae, consistent with the report from a previous study . Second, these three time points allowed H. avenae to maintain the post-parasitic Juvenile 2 (post-J2) stage within both compatible wheat WEN19 and incompatible wheat VP1620, and it was both possible and feasible to rule out the impacts of diverse effectors from different developmental stages of CCN in the comparative analysis because the developmental procedures were not synchronous between WEN19 and VP1620 after 8 days post CCN inoculation (Additional file 1: Figure S1; Additional file 2: Table S1). Time-course examinations were assayed at 24 h, 3 d and 8 d after inoculating VP1620 and WEN19 with freshly hatched J2 CCN. The results showed that the numbers of post-J2 CCN within the VP1620 roots were nearly equivalent to those in WEN19 roots at both 24 h and 3d, whereas the number of post-J2 within VP1620 decreased compared with that in WEN19 at 8 d (Fig. 1b; Table 1). Although VP1620 had no significant effect on the penetration efficiency and early developmental process of the post-J2 stage, it had dramatic impacts on the development of the latter J3 and J4 stages. The numbers of total CCN as well as the corresponding latter J3 and J4 stages within VP1620 roots were significantly decreased at 19 d, 25 d and 33 d compared with those of the compatible control WEN19 (Additional file 2: Table S1). Furthermore, the percentages of J3 and J4 were also obviously delayed and decreased vs. those in WEN19 (Additional file 2: Table S1). The VP1620 failed to hamper the penetration and migration of H. avenae at the early stage but significantly reduced the total number of CCN in its roots and dramatically disturbed the development of the following J3 and J4 stages. These results strongly demonstrated that the nature of the VP1620 underlying resistance to CCN was anti-development.
A novel strategy was established to assemble and analyse transcriptomic data
For the incompatible wheat VP1620, RNA samples of roots with CCN infections at 24 h, 3 d and 8 d and corresponding samples without CCN inoculation were prepared and sequenced by Illumina Hiseq™ 2000. Because the CCNs are typical sedentary endoparasitic nematodes, they remained within the root tissues until cyst formation . Thus, the root samples of both VP1620 and WEN19 collected at 24 h, 3 d and 8 d after CCN infection also contained an amount of post-J2 CCN within these roots. To rule out their negative effects on the following comparative transcriptomic analysis, the corresponding CCN transcripts were removed according to the SOAP method before assembly (Fig. 2a). The samples of 24 h-I_CN, 3d-I_CN and 8d-I_CN eventually yielded 109,831,800, 106,544,764, and 108,397,110 clean reads after filtration, respectively (Table 2). In contrast, the samples of 24 h-I_0, 3d-I_0 and 8d-I_0 did not require removal of CCN transcripts and generated 110,696,256, 105,034,822 and 110,924,676 clean reads, respectively (Table 2). All six RNA-Seq data sets were pooled and assembled using the Trinity program (Fig. 2a). A total of 137,632 unique transcripts (I_transcripts) were generated with a total length of 96,530,157 nt, a mean length of 701 nt and an N50 of 1257 bp (Table 2). Among the 137,632 transcripts, 112,045 transcripts were larger than 200 nt in length with an approximately percentage of 81.41 %.
Similarly, 24 h-C_CN, 3d-C_CN and 8d-C_CN correspondingly produced 106,514,220, 106,119,562 and 106,343,260 clean reads, whereas 24 h-C_0, 3d-C_0 and 8d-C_0 yielded 107,577,970, 109,875,622 and 109,112,640 reads, respectively (Table 2). The transcriptome of compatible wheat was assembled by pooling all six RNA-Seq data sets using the Trinity program (Fig. 2a). There were 142,116 unique transcripts (C_transcripts) with a total length of 97,125,168 nt, a mean length of 683 nt and an N50 of 1214 bp, respectively (Table 2). A total of 115,305 transcripts had a length greater than 200 nt, and the percentage of transcripts that were greater than 200 nt in length was approximately 81.13 %, which was nearly equivalent to the 81.41 % of the I_transcripts. In addition, the C_transcripts had 4484 more unique transcripts than the I_transcripts, and its total length was 595,011 nt longer than that of the I_transcripts. However, both the mean length and N50 of the C_transcripts were somewhat shorter than those of the I_transcripts (Table 2).
Because incompatible (VP1620) and compatible (WEN19) wheat had different genetic backgrounds, it was necessary to first unify their genetic backgrounds before further comparative analysis. The TGI Clustering Tool and BLANTN were applied to unify their genetic backgrounds (Fig. 2a). A total of 12,407 VP1620-specific (I_unique_transcripts), and 13,492 WEN19-specific (C_unique_transcripts) transcripts were produced (Fig. 2b). Additionally, 75,947 transcripts were shared by both cultivars (Overlap_transcripts). In this manner, we studied the remarkable differences in transcription between the compatible and incompatible responses of wheat against H. avenae.
The different responses of compatible or/and incompatible hosts to their plant parasitic nematodes at the early stages were previously investigated [14, 17–19]. However, the transcripts of plant parasitic nematodes were not removed from the raw transcriptomic data, which still contained a certain percentage of endoparasitic nematodes within the roots . In contrast, in this study, we filtered the transcripts of H. avenae from the detected transcriptomics prior to comparative analysis (Fig. 2). This process represents a novel modification that yields precise and clean transcriptomic data and ensures the credibility and accuracy of the transcriptomic data for further analysis.
Upregulated genes in VP1620 induced by H. avenae were identified
To identify upregulated genes of VP1620 as an early defense response of wheat to H. avenae, a comparative transcriptomic analysis of 24 h-I_CN, 3d-I_CN and 8d-I_CN was performed. Groups of 2048, 1924, and 2058 upregulated genes were identified at 24 h, 3 d and 8 d, respectively (Fig. 3a). Among these genes, 889 genes were simultaneously overlapping among 24 h-I_CN, 3d-I_CN and 8d-I_CN and significantly upregulated in at least one of the three time points (Fig. 3a). A total of 1349 genes were common between any two time-point pairs, including 365 genes between 24 h and 3 d, 604 genes between 24 h and 8 d, and 380 genes between 3 d and 8 d. The number of unique genes at 24 h, 3 d and 8 d was 190, 290 and 185, respectively (Fig. 3a).
To explore their “Molecular Function” and “Biological Process”, 889 genes were categorized into subgroups based on GO analysis, and approximately 50 % genes contained GO information. For “Molecular Function”, the “heme binding” subgroup had the largest number with 38 genes (Fig. 3b), a result similar to that of a previous study in which Dap1p (damage resistance protein 1) bound to heme and stabilized the cytochrome P450 protein Erg11p/Cyp51p to result in antifungal resistance . The Dap1p-cytochrome P450 protein pathway is broadly conserved in eukaryotic species [20, 21]. The following two large subcategories were “ammonia-lyase activity” and “oxidoreductase activity” with 23 and 11 genes, respectively (Fig. 3b). The percentage of the top three GO subdivisions accounted for approximately 77.42 %. Defense-related secondary metabolite phenylalanine ammonia-lyase (PAL) activity was significantly induced in response to H. avenae, which was consistent with a previous study in which PAL was remarkably induced in early defense reactions of resistant soybean to H. glycines . The PAL activities were also obviously enhanced within the defense responses to pathogens and aphids [22, 23]. Our results, combined with those of the previous studies, demonstrated that PAL may play vital and general roles against diverse pests. Oxidoreductase activity was dramatically enhanced in response to H. avenae, which was consistent with the responses of soybean cyst nematode and pine wood nematode [24, 25]. An oxidoreductase, i.e., the 2OG-Fe (II) oxygenase family protein with high abundant transcripts, was identified from resistant soybean induced by soybean cyst nematode . Moreover, the “oxidoreductase activity” subcategory was also present with a high percentage in resistant pine infected by pine wood nematode . The high percentage of oxidoreductase strongly suggested that oxidoreductase was likely to play conserved roles in defense responses to diverse plant parasitic nematodes.
The “Biological Process” is primarily composed of diverse catabolic, biosynthetic and metabolic processes. There were 11 significantly enriched subcategories with more than 10 genes in each subcategory (Fig. 3b). A total of 295 upregulated genes were included in these 11 subcategories. It is clear that most of the subcategories were strongly involved in catabolic processes (Fig. 3b), indicating that catabolic processes played important roles in resistance against H. avenae, likely due to sufficient yield and essential precursors for biosynthesis of antimicrobial components. The functional group with the highest number (39 genes) was “organic acid biosynthetic process” (Fig. 3b). The subcategory “benzene-containing compound metabolic process”, which included 25 genes, was also significantly enriched (Fig. 3b), consistent with identification of a similar subcategory within an incompatible reaction of Arabidopsis to Pseudomonas syringae . The phenylpropanoid pathway is considered closely related to plant defense against nematode infection [13, 14, 27]. Identification of the subcategory “phenylpropanoid metabolic process” with 24 genes (Fig. 3b) was consistent with a previous study in which the PAL gene was specifically upregulated in the resistant tomato induced by potato cyst nematode Globodera rostochiensis , suggesting key roles for PAL-mediated metabolism in H. avenae resistance . It is well known that lignin biosynthesis plays vital roles in the defense response . The subcategory “cinnamic acid metabolic process” contained 25 genes, demonstrating that a cinnamic acid metabolic process was likely to participate in H. avenae resistance resulting from stiffening of the cell wall by enhanced lignin biosynthesis.
Monosaccharide transport and syncytium formation are two essential processes for feeding site establishment and maintenance, which serve as the sole nutrient sink in which plant parasitic nematodes develop and parasitize . Two subcategories of “syncytium formation” and “monosaccharide transport” with four and six genes, respectively, were also identified as induced (Fig. 3b), illustrating that VP1620 had a potential ability to feed H. avenae, although it was eventually proven to be highly resistant to CCN, which was consistent with the results mentioned above that VP1620 allowed H. avenae to develop into a few cysts (Fig. 1a).
Among those common 889 genes, 18 genes were dramatically upregulated simultaneously among the three time points (Additional file 3: Figure S2). Most of the 18 vital resistance genes were defense-related components, including one pathogenesis-related protein, one elongation factor-1 α-like protein, one probable mediator of RNA polymerase II transcription subunit, three cell wall-associated hydrolase genes associated with innate immune detection, phytoalexin biosynthesis-related protein premnaspirodiene oxygenase-like, two potent antioxidant protein curcuminoid synthase-like, and lignin biosynthesis involving proteins cationic peroxidase SPC4 and cinnamoyl-CoA reductase 1 (Additional file 4: Table S2).
Candidate resistance genes were identified from overlapping genes
To globally identify resistance genes, expression profiles of all overlapping genes within VP1620 infected by H. avenae (I_CN) were compared with those of the compatible version (C_CN). The expression trends of both I_CN and C_CN displayed nine types at three time points covering the former stage of 24h_vs_3d (from 24 h to 3 d) and the latter stage of 3d_vs_8d (from 3 d to 8 d). Because the time intervals spanned a long period (seven days), a great number of growth and development-related genes were quickly and highly induced and could be detected by comparing C_0 and I_0. To rule out the interfering effects of these genes on identification of resistance genes, we removed these genes before comparing the expression profiles. A total of 67,338 overlapping genes were generated after filtration (Fig. 4). Among the 67,338 overlapping genes, 61,596 genes shared similar expression profiles, as shown outside the brackets on the diagonal line (Fig. 4), which suggested that these genes were likely independent of CCN resistance. However, 5742 genes had uncoordinated expression patterns, as shown outside the brackets beyond the diagonal line (Fig. 4), which indicates that these genes were likely associated with CCN resistance.
To narrow the scopes of the genes that exhibit diverse expression profiles, the genes with expression levels that were obviously simultaneously unchanged among the three time points I_CN compared with VP1620 infected by water (I_0) were wiped off, whereas the genes induced at one of the three time points in I_CN compared with I_0 were retained for further analysis. A total of 1402 genes were ultimately yielded after this procedure (Fig. 4). Among the 1402 genes, 668 genes were found to share similar expression profiles, as shown within the brackets on the diagonal line, whereas 734 genes exhibited uncoordinated expression patterns between I_CN and C_CN, as shown within the brackets beyond the diagonal line (Fig. 4). Obviously, these 734 genes are the most potential to be associated with the wheat resistance to CCN.
It is well known that most of the plant resistance genes were highly induced by diverse pathogens, and thus, in this study, the significantly upregulated genes either in 24h_vs_3d or in 3d_vs_8d were viewed and treated as candidate resistance genes. Among those 734 potential genes described above, a total of 302 genes were dramatically induced with uncoordinated expression profiles between I_CN and C_CN (Fig. 4) and taken as the candidate resistance genes. Among these 302 candidate resistance genes, only 2 genes exhibited continuously enhanced expression trends (U_U) in both the 24h_vs_3d and 3d_vs_8d stages, whereas 254 genes displayed apparently increased expression trends in the 24h_vs_3d stages, including 89 genes with an unchanged expression trend (U_F) and 165 genes with sharply decreased expression tendency (U_D) in the 3d_vs_8d stages; however, 46 genes exhibited dramatically increased expression trends in the 3d_vs_8d stage, including 41 genes that maintained an unaltered expression trend (F_U) and 5 genes that showed obviously decreased expression tendencies (D_U) 24h_vs_3d stages (Fig. 4). Clearly, the number of genes that exhibited significantly increasing expression trends in the 24h_vs_3d stages were increased by approximately 5.5-fold compared to those in the 3d_vs_8d stages, indicating that dramatic transcriptional changes in the early stage were likely to determine the following resistant characteristics.
Phospholipase functions as a central regulator in resistance against H. avenae
To obtain an overview of the processes altered during the early stages of VP1620 in defense as a response to H. avenae infestation, the 302 candidate resistance genes identified above were functionally classified by KEGG analysis. A total of 224 genes out of 302 genes, including 169 genes that exhibited increasing expression trends in the 24h_vs_3d stage and 55 genes that displayed increasing expression trends in the 3d_vs_8d stage, were categorized into 50 pathways (Additional file 5: Table S3), among which 25 pathways had at least two genes (Fig. 5a). Six pathways, i.e., “Lysosome”, “Phagosome”, “ABC transport”, “RNA degradation”, “Cysteine and methionine metabolism” and “Biosynthesis of amino acid”, contained KEGG information only in the 24h_vs_3d stage with a total number of 14 genes (Fig. 5a). In contrast, the two pathways of “Protein processing in endoplasmic” and “Tryptophan metabolism” contained KEGG information only in the 3d_vs_8d stage with a total number of four genes (Fig. 5a). The “Glycerophospholipid metabolism” pathway was enriched with the largest number of members (21 genes), among which were 16 genes in the 24h_vs_3d stage and five genes in the 3d_vs_8d stage. A total of 20 genes were enriched in “Ether lipid metabolism”, including 15 and five genes in the 24h_vs_3d and 3d_vs_8d stages, respectively. Another eight enriched KEGG-pathways contained members beyond 10 genes (Fig. 5a). Our results indicated that the number of genes dramatically enhanced during the former stage (24h_3d) was significantly larger than that in the latter stage (3d_8d), which supports the hypothesis that early reaction of the incompatible defense would determine the subsequent and even the final phenotypes or resistant characteristics.
Endocytosis has been demonstrated as involved in plant immunity through many protein–protein interactions [29, 30]. The “Endocytosis” pathway was enriched with 15 genes in the 24h_vs_3d stage and six genes in the 3d_vs_8d stage (Fig. 5a), indicating that endocytosis is also likely involved in the defense response to CCN infestation. It has been reported that lysosomes are involved in plant defense by participation in autophagy leading and programmed cell death contributing to the defense response . Unlike “Endocytosis”, both “Lysosome” and “Phagosome” had only three genes within the 24h_vs_3d stage (Fig. 5a), suggesting that lysosomes were likely to function during the early defense reaction to CCN attack.
Interestingly, in this study, two genes involved in the “benzoxazinoid biosynthesis” pathway, DIMBOA (2, 4-dihydroxy-7-methoxy-1, 4-benzoxazin-3-one) which is a naturally occurring hydroxamic acid, and benzoxazinoid, were identified (Fig. 5a). DIMBOA is a powerful antibiotic present in maize and related grasses, particularly wheat, and serves as a natural defense against a wide range of pests, including insects, pathogenic fungi and bacteria. The identification of significantly upregulated expression genes involved in the “benzoxazinoid biosynthesis” pathway showed that DIMBOA was also likely to play an important role in defense against plant parasitic nematodes.
The RNA metabolism was likely involved in resistance to H. avenae. The RNA methylation probably plays vital roles in resistance to H. avenae, similar to a new resistance mechanism of resistance against aminoglycosides among gram-negative pathogens according to methylation of ribosomal RNA . In genetic pathways, the “RNA transport”, “mRNA surveillance pathway” “DNA repair and recombination” pathways contained 14, 12 and 7 genes, respectively, in the former 24h_vs_3d stage, and 7, 7 and 1 gene, respectively, in the latter 3d_vs_8d stage (Fig. 5a). In this study, the “RNA methylation” subgroup included 11 genes that were also identified (Fig. 3b), indicating that RNA methylation likely played a vital role in resistance to H. avenae. It has been reported that a candidate resistance gene of maize against Aspergillus flavus was a component of the nuclear pore complexes (NPCs) in the nuclear envelope, and NPCs were involved in RNA transport . The RNA transport pathways are composed of various protein complexes that regulate gene expression and nucleocytoplasmic trafficking. The nucleocytoplasmic trafficking pathways are fundamental for normal cell functions as well as plant defense responses . Previous studies demonstrated that RNA transport pathway genes played direct roles in plant defense systems [33, 35]. The mRNA surveillance pathway, also known as nonsense-mediated mRNA decay (NMD), represents a splicing- and translation-dependent pathway of RNA degradation, which limits the expression of transcripts bearing premature translation termination codons [36–38]. In addition, the DNA repair and recombination pathway was also enriched (Fig. 5a). The DNA repair and recombination pathways have many functions in different species, i.e., maintenance of genomic integrity, chromosome segregation and recombination, and immune system development. The Arabidopsis RAR genes include homologues of many DNA repair genes that are defective in different human diseases [39–41]. To the best of our knowledge, this paper is the first report that RNA and DNA metabolism is likely involved in plant immunity against plant parasitic nematodes.
The KEGG pathway enrichment tests showed that the cAMP signaling pathway, Ras signaling pathway, glycerophospholipid metabolism, and ether lipid metabolism contained large numbers of gene members (Fig. 5a). The cAMP signaling pathway involves the stress response through accumulation of phytoalexin and ethylene production by regulating Ca2+ or K+ flux and thereby initiates Ca2+ and protein kinase signaling cascades . Ras are signaling proteins known to transmit extracellular signals in yeast and animals and have a crucial role in cellular signaling in animals and various lower eukaryotes [43, 44]. Recently, a plant intracellular Ras-group LRR family in Arabidopsis was identified and proven to function in development . Members of the Ras superfamily share several common structural features, including four guanine nucleotide binding domains and an effector binding domain, and potentially function in intracellular defense signal transduction via protein-protein interactions mediated by these typical domains . Glycerophospholipid serves as a structural component of cell membranes, which provides a robust permeability barrier and a first line of defense against different stresses; additionally, the glycerophospholipid pathway was previously revealed as critical for programmed cell death, which is involved in the defense response to attacks from diverse pathogens [46, 47]. Ether lipid metabolism was likely to associate with the lipoxygenase (LOX) pathway to generate antimicrobial compounds such as the divinyl ether fatty acids, leading to defense responses . Furthermore, our results also demonstrated that several vital LOX genes were also significantly upregulated and induced by H. avenae infestation (Fig. 6).
Phospholipase D (PLD) has been shown to be involved in both abiotic and biotic stress signaling [49, 50]. Interestingly, in this study, phospholipase D1/2 coding genes (PLD) were found in many of the vital and large KEGG pathways mentioned previously, i.e., “Endocytosis”, “cAMP signaling pathway”, “Ras signaling pathway”, “Glycerophospholipid metabolism”, and “Ether lipid metabolism” (Fig. 5b). Moreover, PLD localized in the connection position and linked these pathways to form a complex network (Fig. 5b). This finding strongly suggested that phospholipase was likely to be as a central regulator in resistance to cereal cyst nematode H. avenae.
It has been reported that phospholipase positively mediated plant defense response through its role in resistance signaling associated with reactive oxygen species (ROS) [49–51]. To characterize the role of ROS in defense response of I_CN at the early stages, we assayed the ROS contents within VP1620 infected by CCN at 24 h, 3 d and 8 d. The results showed two ROS bursts, including a strong ROS burst at 24 h and an obviously lower ROS burst at 3 d (Fig. 6a). However, no dramatic change of ROS content was detected at 8 d (Fig. 6a). Moreover, 19 ROS-producing genes, including 13 Class III peroxidase (POX) genes and six lipoxygenase (LOX) genes, were identified from overlapping genes between _I and _C (Additional file 6: Table S4). The qRT-PCR results showed that three of 13 POXs were significantly upregulated at 24 h after CCN infestation, and their changes in relative expression were more than two-fold; however, these three genes had no obvious induction of transcript abundance at either 3 d or 8 d (Fig. 6b). Compared with a low percentage (23 %) of POX genes responsible for ROS production, five of the six (83 %) LOX genes identified herein were dramatically upregulated at 24 h after CCN inoculation. In addition, LOX2 and LOX4 were also significantly induced by CCN at 3 d (Fig. 6b). These results showed that significant upregulation of POX combined with LOX were responsible for the former strong ROS production at 24 h, and only the dramatic induction of LOXs contributed to the latter ROS burst. As a whole, lipoxygenase was primarily responsible for ROS bursts in VP1620 after CCN infestation. Our results further supported the roles of phospholipase in wheat defense responses to H. avenae.
Plant-pathogen interaction was significantly enriched
In this study, “plant-pathogen interaction (PPI)” was enriched with the largest numbers (20 genes), including 18 genes in the 24h_vs_3d stage and two genes in the 3d_vs_8d stage, respectively (Fig. 5a). Seven key defense-related genes were dramatically upregulated (Fig. 7). Among these seven genes, five genes were found to be significantly upregulated in the 24h_vs_3d stage, including two WRKY transcription factor genes (WRKY25, WRKY29), one MAPK cascade gene (BAK1), one Ca2+ calmodulin-like protein (CaMC/ML) gene, and one CERK1 gene; however, only one MKK4/5 gene (two redundant MAPKKs) showed dramatically enhanced expression in the 3d_vs_8d stage. Furthermore, one gene RPM1 was remarkably and continuously upregulated in both 24h_vs_3d and 3d_vs_8d stages (Fig. 7).
Chitin is a main component of the cell wall of fungal pathogens [52, 53], the exoskeletal cuticle and gut lining of insects [54, 55], and the egg shell of the nematode [56, 57]. Previous studies indicated that CEBiP (chitin elicitor-binding protein) and its partner LysM-RLK (receptor-like kinase) chitin elicitor receptor kinase1 (CERK1) are involved in plant immunity by recognizing the chitin oligosaccharides of the cell wall of fungal pathogens through extracellular LysM domain [58, 59]. In this study, ten genes were identified with LysM domain, and only one gene CERK1, which not only had orthologs to both AtCERK1 and OsCERK1 but also had an extra tyrosine kinase domain, was remarkably enhanced in the incompatible response against CCN infestation (both 3dpi and 8dpi) (Fig. 7), while the other nine genes had no orthologs to AtCERK1 or OsCERK1. One of these nine genes was obviously down-regulated (3dpi), while the other eight genes had no significant changes in their expression level. These results indicated CERK1 played potential roles in defense responses to CCN. Recently, it was reported that chitin-binding protein was also associated with plant defense against the insect Hessian fly . Our finding of CERK1 gene associated with defense response to H. avenae, combined with previous studies on similar functions with fungal pathogens and insect, indicated that CERK1 activated defense-related genes by recognizing chitin not only from fungal pathogens but also from insect and plant parasitic nematode, which was treated as a supplemental view for current perspectives. It has been shown that the receptor-like kinase (RLK) pattern recognition receptor (PRR) protein CERK1 perceives the fungal pathogen-activated molecular patterns (PAMPs) chitin to trigger a PAMP-triggered immunity (PTI) response . The PTI signaling requires its dimerization with the receptor-like protein CEBiP1 through direct binding with chitin [61, 62]. In this study, it is likely that putative CERK1 was also involved in PTI response of incompatible wheat to CCN, possibly through recognition of the chitin component of CCN. Plant CaM/CML induced downstream nitric oxide (NO) synthesis as an intermediary sites in a pathogen perception signaling cascade, eventually leading to innate plant immune responses . Plant MAPKs (mitogen-activated protein kinases) cascades play pivotal roles in regulating plant development and signaling responses to a variety of stress stimuli. Activation of MAPKs is performed by their upstream kinases, i.e., MAPK kinases (MAPKKs), and MKK4/MKK5 as a type of MAPKKs acted upstream of MPK3/6 in regulating plant development and defense responses by regulating jasmonic acid (JA) signal transduction . In this study, pathogenesis-related (PR) genes orthologous to PR4 and PR10, which involves the defense meditate by JA pathways [64, 65], were significantly up-regulated at different times. However, most of the other PR genes were not obviously changed in their expression level in the event of CCN attack. These results indicated that JA pathway may play a pivotal role in wheat to CCN similar to that of rice defense against root knot nematodes .
The BAK1 (brassinosteroid insensitive 1 associated kinase 1) is a positive regulator of the brassinosteroid (BR) hormonal signaling pathway . The BAK1 was significantly up-regulated (Fig. 7), while both FLS2 (FLAGELLIN SENSITIVE2) and EFR (elongation factor Tu (EF-Tu) receptor) were not obviously changed. FLS2 and EFR are two vital receptors to recognize the specific components of bacterial pathogens to activate PAMP-triggered immunity [66, 67]. So, the BAK1 is a ligand-independent co-receptor of RLKs such as FLS2 and EFR and acts as a central regulator of PTI triggered by diverse PAMPs . The enhanced expression of both MKK4/MKK5 and BAK1 in this study resulted in PTI of incompatible wheat to CCN, likely according to association with the plant hormone network. Many WRKY genes were proven as involved in innate plant immune responses and showed increased expression in response to diverse pathogens, i.e., fungal and bacterial pathogens [61, 62], soybean cyst nematode  and root knot nematode , as well as cereal cyst nematode H. avenae in this study. The RPM1 encodes a peripheral membrane protein with LZ-NBS-LRR domains and confers resistance to specific Pseudomonas syringae strains . In addition, the locus of leaf rust resistance gene Lr10 was similar to the Arabidopsis RPM1 locus . Resistant complex components of RAR1, SGT1 and HSP90 played vital roles in defense response to diverse pests. RAR1 was originally isolated from barley Hordeum vulgare and is present in eukaryotes except for yeast Saccharomyces cerevisiae. RAR1 was required for a subset of R-gene-mediated resistance responses in monocot and dicot plant species [72, 73]. SGT1 interacting with RAR1 also contributed to R-gene-mediated resistance . The complex of both RAR1 and SGT1 interacting with chaperon HSP90 functioned in various signal transduction networks [74, 75]. However, the complex of HSP90 and SGT1 (except for RAR1) was required for Mi-1-mediated aphid and nematode resistance , providing further evidence for common components in early resistance defense signaling against diverse pathogens and pests. In this study, the PRM1 was significantly upregulated in both the 24h_3d and 3d_8d stages, strongly suggesting that PRM1 likely served as a constitutive defense gene involved in the resistance to H. avenae. Whether the total or partial complex components of SGT1, RAR1 and HSP90 are required for resistance to H. avenae must be further investigated.
Additional rescued candidate resistance genes were identified
For genes that exhibit similar expression profiles, their scales of altered expression levels could remarkably differ between I_CN and C_CN. To identify additional resistance genes from the genes that exhibited similar expression profiles, the expression profiles were analyzed using the scales of altered expression levels. Among the 668 genes exhibiting similar expression profiles (Fig. 4), 182 genes were rescued and isolated with uncoordinated expression patterns, and 486 genes were identified with similar expression profiles using the scales of altered expression levels as a standard (Additional file 7: Table S5). Among the overlapping genes, most of the candidate resistance genes were identified from uncoordinated expression patterns according to transcript abundance as a standard. Additional resistance genes were subsequently rescued from similar expression profiles with the scales of transcriptional changes as a standard. In this study, we combined both strategies to develop a further comprehensive methodology to identify additional candidate resistance genes not only from uncoordinated expression profiles as the main section but also from similar expression patterns as a supplementary component.
The KEGG pathway analysis revealed that “phenylalanine metabolism”, “phenylpropanoid biosynthesis” and “diterpenoid biosynthesis” involving phytoalexin biosynthesis were enriched with four or five genes, respectively, and “RNA transport” and “mRNA surveillance pathway” were also dramatically enriched with 14 and six genes, respectively. In contrast, only one member of the LRR receptor-like serine/threonine-protein kinase FLS2 was identified with typical characteristics of a resistance gene in “plant-pathogen interaction pathway” (Additional file 8: Table S6). Moreover, defense-related pathways “Endocytosis”, “cAMP signaling pathway”, “Ras signaling pathway”, “Glycerophospholipid metabolism”, and “Ether lipid metabolism” were also significantly enriched, and these pathways were cross-linked to form an extraordinary network by the hub switch phospholipase D1/2, similar to findings mentioned previously (Fig. 7).
VP1620-specific resistance genes were also identified
To identify VP1620-specific defense-related genes as the candidate resistant genes, incompatible unique (I_unique) transcripts were classified using expression profiles with transcript abundance as a standard. A total of 150 genes were significantly upregulated at one of the three time points in I_CN compared with I_0 (Additional file 9: Figure S3). Therefore, these 122 genes can be taken as the candidate resistance genes. Among these genes, 122 genes had uncoordinated expression profiles, as shown beyond the diagonal line, and 28 genes shared similar expression profiles, as shown on the diagonal line (Additional file 9: Figure S3). The number of genes exhibiting uncoordinated expression profiles was approximately four times that of similar expression profiles. It is clear that 88 % of the 150 genes (132 genes) had an F_F expression patterns in I_0 vs. those that had diverse expression profiles in I_CN. For these 132 genes, except for 28 genes that showed similar expression trends between I_CN and I_0, 63 genes out of the left 104 genes displaying uncoordinated expression profiles were significantly upregulated in at least one of the three time points (Additional file 9: Figure S3). A total of 17 vital VP1620-specific defense related genes were identified (Table 3). Among these 17 genes, seven protein kinase coding genes, including three proline-rich receptor-like protein kinases, one G-type lectin S-receptor-like serine/threonine-protein kinase, one protein-kinase-domain-containing protein, one cysteine-rich receptor-like protein kinase and one probable LRR receptor-like serine/threonine-protein kinase (except for one proline-rich receptor-like protein kinase) exhibited similar expression profiles U_D in I_CN (Table 3). Among the three transcription factor genes, one MADS-box transcription factor displayed a U_D expression pattern, whereas the other two ethylene-responsive transcription factor ABI4 exhibited U_F in I_CN (Table 3). Among the seven resistance genes, two leucine-rich repeat (LRR) genes and one resistance RGA2 gene exhibited U_D expression, and two hydroxyproline-rich glycoprotein genes displayed the corresponding F_U and U_D expression profiles in I_CN (Table 3). Moreover, the results of KEGG pathway enrichment analysis also indicated that phospholipase in partially VP620-specific components was involved in many defense-related pathways similar to those mentioned previously.
This study, compared to identification of resistance genes against other plant parasitic nematodes, is the first to present large-scale identification of wheat resistance genes against cereal cyst nematode H. avenae using comparative transcriptomic analysis of incompatible and compatible reactions. In total, 606 candidate resistance genes including 302 primarily identified overlapping genes, 184 additional rescued genes, and 122 VP1620-specific genes, were identified on a large scale and functionally classified into diverse defense-related pathways. Phospholipases might play vital roles in early defense responses to H. avenae by involving diverse defense-related pathways as a hub switch. In addition, the plant-pathogen interaction pathway was significantly enriched and also played key roles in defense against H. avenae. These data not only provide applicable candidate resistance genes for breeding novel wheat cultivars but also strengthen new insights into defense responses to H. avenae, especially a better understanding of the mechanism underlying early defense response of cereal crops to cereal cyst nematode H. avenae.
Materials and methods
Nematode infection manipulation
Cereal cyst nematode H. avenae was maintained on a compatible wheat WEN19 at 16 °C in a greenhouse, mature cysts were collected, and pre-J2 were hatched as previously described [2, 8]. Both WEN19 and VP1620 were maintained at 16 °C in a greenhouse, and one-week-old wheat seedlings were inoculated with freshly hatched H. avenae as described previously . To block continuous infection of H. avenae, the roots of both groups of infected plants were rinsed with double distilled water, and the rinsed plants were transferred to new pots for further growth . The infected root samples of both WEN19 and VP1620 were harvested at 24 h, 3 d and 8 d, and the corresponding root samples without CCN inoculation were also collected as controls. Among these root samples, a subset were used to examine the number of different developmental CCN within the roots by staining the CCN with 0.01 % acid fuchsin solution as previously described [76, 77]; the others were subjected to RNA isolation with TRIzol reagent following the manufacturer’s instructions (Invitrogen).
Total RNAs were subjected to RNA-sequencing by Illumina Hiseq™ 2000. To remove CCN transcripts from raw transcriptomes, RNA-Seq data was mapped to CCN transcripts using SOAP , and unmapped reads from each sample were used for assembly in the next step. Transcriptome assemblies were performed with Trinity . The RNA-Seq data from each sample of _I or _C was assembled together. Transcripts from _I and _C were used in a further process for redundancy removal using the TGI Clustering Tool (http://sourceforge.net/projects/tgicl/) and BLASTN. The acquired overlap and unique transcripts were used for further analyses. The BLASTX alignments (e-value < 0.00001) between transcripts and protein databases (NR, Swiss-Prot, KEGG and COG) were performed. The best hits were used to decide on sequence direction of transcripts. If hits of different databases conflicted with each other, a priority order of NR, SwissProt, KEGG and COG was followed for deciding the sequence direction of transcripts. For transcripts with no hits in these databases, ESTScan  was introduced to decide their sequence directions. To annotate functions of transcripts, they were searched against databases NT (BLASTN, e-value < 1e-5), NR, Swiss-Prot, KEGG and COG (BLASTX, e-value < 1e-5). The Blast2GO program  was used to obtain Gene Ontology annotation of transcripts, and metabolic pathways of transcripts were analyzed by Path_finder . The expressions of genes were measured by Fragments Per kb per Million Fragments (FPKM) . The algorithm used to identify differentially expressed genes was similar to that previously reported . The calculated p-values were corrected using the Bonferroni method. A gene was considered to be differentially expressed between two samples if the corrected p-value (FDR) ≤ 0.001 and the fold change ≥ 2.
Measurement of ROS level
The ROS levels within root samples were determined via luminol-chemiluminescence assay with minor modifications . In a brief, root samples of both WEN19 and VP1620 collected at 24 h, 3 d and 8 d were cut and dropped into double distilled water for 4 h. An amount of 10 mg root per sample was placed in Eppendorf tubes containing 100 μL of luminal (Bio-Rad Immun-Star horseradish peroxidase substrate) and 1.0 μL of horseradish peroxidase (Jackson ImmunoResearch). Luminescence was immediately measured at 1 min intervals for 30 min with a Glomax 20/20 luminometer (Promega) . Three replicates were performed for each sample and treatment. Standard errors were calculated for each treatment.
A total of 12 RNAs were treated with DNase to remove contaminating gDNA. First-strand cDNA was synthesized using the Invitrogen Superscript III Kit [52, 86], and qRT-PCR assay was performed using the iQ5 real-time PCR detection system (Bio-Rad). The actin gene was used as endogenous reference , and relative expression was determined as previously described . Mean and standard errors were determined with data from three independent replicates. The primers used in this study are listed in Additional file 6: Table S4.
Cereal cyst nematode
Chitin elicitor-binding proteins
Elongation factor Tu (EF-Tu) receptor
Mitogen-acitivated protein kinases
Pathogen-activitated molecular patterns
Post-parasitic Juvenile 2
Class III peroxidase
Reactive oxygen species
Kumar M, Gantasala NP, Roychowdhury T, Thakur PK, Banakar P, Shukla RN, et al. De novo transcriptome sequencing and analysis of the cereal cyst nematode, Heterodera avenae. PLoS ONE. 2014;9(5):e96311.
Long HB, Peng DL, Huang WK, Peng H, Wang GF. Molecular characterization and functional analysis of two new β-1,4-endoglucanase genes (Ha-eng-2, Ha-eng-3) from the cereal cyst nematode Heterodera avenae. Plant Pathol. 2013;62:953–60.
Simonetti E, Alba E, Montes MJ, Delibes A, Lopez-Brana I. Analysis of ascorbate peroxidase genes expressed in resistant and susceptible wheat lines infected by the cereal cyst nematode, Heterodera avenae. Plant Cell Rep. 2010;29(10):1169–78.
Delibes A, Romero D, Aguaded S, Duce A, Mena M, Lopez-Brana I, et al. Resistance to the cereal cyst nematode (Heterodera avenae Woll.) transferred from the wild grass Aegilops ventricosa to hexaploid wheat by a "stepping-stone" procedure. Theor Appl Genet. 1993;87(3):402–8.
Seah S, Miller C, Sivasithamparam K, Lagudah ES. Root responses to cereal cyst nematode (Heterodera avenae) in hosts with different resistance genes. New Phytol. 2000;146:527–33.
Safari E, Gororo NN, Eastwood RF, Lewis J, Eagles HA, Ogbonnaya FC. Impact of Cre1, Cre8 and Cre3 genes on cereal cyst nematode resistance in wheat. Theor Appl Genet. 2005;110(3):567–72.
Majnik MJ, Ogbonnaya FC, Moullet O, Lagudah ES. The Cre1 and Cre3 nematode resistance genes are located at homeologous loci in the wheat genome. Mol Plant Microbe Interact. 2003;16(12):1129–34.
Xu DL, Long H, Liang JJ, Zhang J, Chen X, Li JL, et al. De novo assembly and characterization of the root transcriptome of Aegilops variabilis during an interaction with the cereal cyst nematode. BMC Genomics. 2012;13:133.
Taylor C, Shepherd KW, Langridge P. A molecular genetic map of the long arm of chromosome 6R of rye incorporating the cereal cyst nematode resistance gene, CreR. Theor Appl Genet. 1998;97:1000–12.
Kretschmer JM, Chalmers KJ, Manning S, Karakousis A, Barr AR, Islam AK, et al. RFLP mapping of the Ha 2 cereal cyst nematode resistance gene in barley. Theor Appl Genet. 1997;94(2):1060–4.
Schaff JE, Nielsen DM, Smith CP, Scholl EH, Bird DM. Comprehensive transcriptome profiling in tomato reveals a role for glycosyltransferase in Mi-mediated nematode resistance. Plant Physiol. 2007;144(2):1079–92.
Klink VP, Overall CC, Alkharouf NW, MacDonald MH, Matthews BF. A time-course comparative microarray analysis of an incompatible and compatible response by Glycine max (soybean) to Heterodera glycines (soybean cyst nematode) infection. Planta. 2007;226(6):1423–47.
Mazarei M, Liu W, Al-Ahmad H, Arelli PR, Pantalone VR, Stewart CJ. Gene expression profiling of resistant and susceptible soybean lines infected with soybean cyst nematode. Theor Appl Genet. 2011;123(7):1193–206.
Klink VP, Hosseini P, Matsye PD, Alkharouf NW, Matthews BF. Syncytium gene expression in Glycine max([PI 88788]) roots undergoing a resistant reaction to the parasitic nematode Heterodera glycines. Plant Physiol Biochem. 2010;48(2–3):176–93.
Beneventi MA, da Silva OJ, de Sa ME, Firmino AA, de Amorim RM, Albuquerque EV, et al. Transcription profile of soybean-root-knot nematode interaction reveals a key role of phythormones in the resistance reaction. BMC Genomics. 2013;14:322.
Wan J, Vuong T, Jiao Y, Joshi T, Zhang H, Xu D, et al. Whole-genome gene expression profiling revealed genes and pathways potentially involved in regulating interactions of soybean with cyst nematode (Heterodera glycines Ichinohe). BMC Genomics. 2015;16(1):148.
Klink VP, Hosseini P, Matsye P, Alkharouf NW, Matthews BF. A gene expression analysis of syncytia laser microdissected from the roots of the Glycine max (soybean) genotype PI 548402 (Peking) undergoing a resistant reaction after infection by Heterodera glycines (soybean cyst nematode). Plant Mol Biol. 2009;71(6):525–67.
Szakasits D, Heinen P, Wieczorek K, Hofmann J, Wagner F, Kreil DP, et al. The transcriptome of syncytia induced by the cyst nematode Heterodera schachtii in Arabidopsis roots. Plant J. 2009;57(5):771–84.
Barcala M, Garcia A, Cabrera J, Casson S, Lindsey K, Favery B, et al. Early transcriptomic events in microdissected Arabidopsis nematode-induced giant cells. Plant J. 2010;61(4):698–712.
Craven RJ, Mallory JC, Hand RA. Regulation of iron homeostasis mediated by the heme-binding protein Dap1 (damage resistance protein 1) via the P450 protein Erg11/Cyp51. J Biol Chem. 2007;282(50):36543–51.
Mallory JC, Crudden G, Johnson BL, Mo C, Pierson CA, Bard M, et al. Dap1p, a heme-binding protein that regulates the cytochrome P450 protein Erg11p/Cyp51p in Saccharomyces cerevisiae. Mol Cell Biol. 2005;25(5):1669–79.
Chaman ME, Copaja SV, Argandona VH. Relationships between salicylic acid content, phenylalanine ammonia-lyase (PAL) activity, and resistance of barley to aphid infestation. J Agric Food Chem. 2003;51(8):2227–31.
Kim DS, Hwang BK. An important role of the pepper phenylalanine ammonia-lyase gene (PAL1) in salicylic acid-dependent signalling of the defence response to microbial pathogens. J Exp Bot. 2014;65(9):2295–306.
Hewezi T, Baum TJ. Manipulation of plant cells by cyst and root-knot nematode effectors. Mol Plant Microbe Interact. 2013;26(1):9–16.
Hirao T, Fukatsu E, Watanabe A. Characterization of resistance to pine wood nematode infection in Pinus thunbergii using suppression subtractive hybridization. BMC Plant Biol. 2012;12:13.
Frei dit Frey N, Garcia AV, Bigeard J, Zaag R, Bueso E, Garmier M, et al. Functional analysis of Arabidopsis immune-related MAPKs uncovers a role for MPK3 as negative regulator of inducible defences. Genome Biol. 2014;15(6):R87.
Uehara T, Sugiyama S, Matsuura H, Arie T, Masuta C. Resistant and susceptible responses in tomato to cyst nematode are differentially regulated by salicylic acid. Plant Cell Physiol. 2010;51(9):1524–36.
Tronchet M, Balague C, Kroj T, Jouanin L, Roby D. Cinnamyl alcohol dehydrogenases-C and D, key enzymes in lignin biosynthesis, play an essential role in disease resistance in Arabidopsis. Mol Plant Pathol. 2010;11(1):83–92.
Robatzek S. Vesicle trafficking in plant immune responses. Cell Microbiol. 2007;9(1):1–8.
Engelhardt S, Boevink PC, Armstrong MR, Ramos MB, Hein I, Birch PR. Relocalization of late blight resistance protein R3a to endosomal compartments is associated with effector recognition and required for the immune response. Plant Cell. 2012;24(12):5142–58.
Dickman MB, Fluhr R. Centrality of host cell death in plant-microbe interactions. Annu Rev Phytopathol. 2013;51:543–70.
Doi Y, Arakawa Y. 16S ribosomal RNA methylation: emerging resistance mechanism against aminoglycosides. Clin Infect Dis. 2007;45(1):88–94.
Kelley RY, Williams WP, Mylroie JE, Boykin DL, Harper JW, Windham GL, et al. Identification of maize genes associated with host plant resistance or susceptibility to Aspergillus flavus infection and aflatoxin accumulation. PLoS ONE. 2012;7(5):e36892.
Asters MC, Williams WP, Perkins AD, Mylroie JE, Windham GL, Shan X. Relating significance and relations of differentially expressed genes in response to Aspergillus flavus infection in maize. Sci Rep. 2014;4:4815.
Cheng YT, Germain H, Wiermer M, Bi D, Xu F, Garcia AV, et al. Nuclear pore complex component MOS7/Nup88 is required for innate immunity and nuclear accumulation of defense regulators in Arabidopsis. Plant Cell. 2009;21(8):2503–16.
Tian DQ, Pan XY, Yu YM, Wang WY, Zhang F, Ge YY, et al. De novo characterization of the Anthurium transcriptome and analysis of its digital gene expression under cold stress. BMC Genomics. 2013;14:827.
Bai TT, Xie WB, Zhou PP, Wu ZL, Xiao WC, Zhou L, et al. Transcriptome and expression profile analysis of highly resistant and susceptible banana roots challenged with Fusarium oxysporum f. sp. cubense tropical race 4. PLoS ONE. 2013;8(9):e73945.
Wu J, Kang JH, Hettenhausen C, Baldwin IT. Nonsense-mediated mRNA decay (NMD) silences the accumulation of aberrant trypsin proteinase inhibitor mRNA in Nicotiana attenuata. Plant J. 2007;51(4):693–706.
Singh SK, Choudhury SR, Roy S, Sengupta DN. Understanding DNA repair and recombination in higher plant genome: information from genome-wide screens in Arabidopsis and rice. Plant Signal Behav. 2011;6(1):120–2.
Initiative AG. Analysis of the genome sequence of the flowering plant Arabidopsis thaliana. Nature. 2000;408(6814):796–815.
Singh SK, Roy S, Choudhury SR, Sengupta DN. DNA repair and recombination in higher plants: insights from comparative genomics of Arabidopsis and rice. BMC Genomics. 2010;11:443.
Zhao J, Guo Y, Kosaihira A, Sakai K. Rapid accumulation and metabolism of polyphosphoinositol and its possible role in phytoalexin biosynthesis in yeast elicitor-treated Cupressus lusitanica cell cultures. Planta. 2004;219(1):121–31.
Forsthoefel NR, Cutler K, Port MD, Yamamoto T, Vernon DM. PIRLs: a novel class of plant intracellular leucine-rich repeat proteins. Plant Cell Physiol. 2005;46(6):913–22.
Forsthoefel NR, Klag KA, Simeles BP, Reiter R, Brougham L, Vernon DM. The Arabidopsis plant intracellular Ras-group LRR (PIRL) family and the value of reverse genetic analysis for identifying genes that function in gametophyte development. Plants. 2013;2:507–20.
Yang Z. Small GTPases: versatile signaling switches in plants. Plant Cell. 2002;14(Suppl):S375–88.
Reina-Pinto JJ, Yephremov A. Lipid determinants of cell death. Plant Signal Behav. 2009;4(7):625–8.
Dalebroux ZD, Matamouros S, Whittington D, Bishop RE, Miller SI. PhoPQ regulates acidic glycerophospholipid content of the Salmonella Typhimurium outer membrane. Proc Natl Acad Sci U S A. 2014;111(5):1963–8.
Kolomiets MV, Chen H, Gladon RJ, Braun EJ, Hannapel DJ. A leaf lipoxygenase of potato induced specifically by pathogen infection. Plant Physiol. 2000;124(3):1121–30.
Wang X. Regulatory functions of phospholipase D and phosphatidic acid in plant growth, development, and stress responses. Plant Physiol. 2005;139(2):566–73.
Pinosa F, Buhot N, Kwaaitaal M, Fahlberg P, Thordal-Christensen H, Ellerstrom M, et al. Arabidopsis phospholipase Dδ is involved in basal defense and nonhost resistance to powdery mildew fungi. Plant Physiol. 2013;163(2):896–906.
Testerink C, Munnik T. Phosphatidic acid: a multifunctional stress signaling lipid in plants. Trends Plant Sci. 2005;10(8):368–75.
Kong LA, Yang J, Li GT, Qi LL, Zhang YJ, Wang CF, et al. Different chitin synthase genes are required for various developmental and plant infection processes in the rice blast fungus Magnaporthe oryzae. PLoS Pathog. 2012;8(2):e1002526.
Lenardon MD, Munro CA, Gow NA. Chitin synthesis and fungal pathogenesis. Curr Opin Microbiol. 2010;13(4):416–23.
Zhu YC, Specht CA, Dittmer NT, Muthukrishnan S, Kanost MR, Kramer KJ. Sequence of a cDNA and expression of the gene encoding a putative epidermal chitin synthase of Manduca sexta. Insect Biochem Mol Biol. 2002;32(11):1497–506.
Tellam RL, Vuocolo T, Johnson SE, Jarmey J, Pearson RD. Insect chitin synthase cDNA sequence, gene organization and expression. Eur J Biochem. 2000;267(19):6025–43.
Zhang Y, Foster JM, Nelson LS, Ma D, Carlow CK. The chitin synthase genes chs-1 and chs-2 are essential for C. elegans development and responsible for chitin deposition in the eggshell and pharynx, respectively. Dev Biol. 2005;285(2):330–9.
Fanelli E, Di Vito M, Jones JT, De Giorgi C. Analysis of chitin synthase function in a plant parasitic nematode, Meloidogyne artiellia, using RNAi. Gene. 2005;349:87–95.
Miya A, Albert P, Shinya T, Desaki Y, Ichimura K, Shirasu K, et al. CERK1, a LysM receptor kinase, is essential for chitin elicitor signaling in Arabidopsis. Proc Natl Acad Sci U S A. 2007;104(49):19613–8.
Tanaka S, Ichikawa A, Yamada K, Tsuji G, Nishiuchi T, Mori M, et al. HvCEBiP, a gene homologous to rice chitin receptor CEBiP, contributes to basal resistance of barley to Magnaporthe oryzae. BMC Plant Biol. 2010;10:288.
Giovanini MP, Saltzmann KD, Puthoff DP, Gonzalo M, Ohm HW, Williams CE. A novel wheat gene encoding a putative chitin-binding lectin is associated with resistance against Hessian fly. Mol Plant Pathol. 2007;8(1):69–82.
Liu W, Liu J, Triplett L, Leach JE, Wang GL. Novel insights into rice innate immunity against bacterial and fungal pathogens. Annu Rev Phytopathol. 2014;52:213–41.
Meng X, Zhang S. MAPK cascades in plant disease resistance signaling. Annu Rev Phytopathol. 2013;51:245–66.
Ma W, Smigel A, Tsai YC, Braam J, Berkowitz GA. Innate immunity signaling: cytosolic Ca2+ elevation is linked to downstream nitric oxide generation through the action of calmodulin or a calmodulin-like protein. Plant Physiol. 2008;148(2):818–28.
Nahar K, Kyndt T, De Vleesschauwer D, Hofte M, Gheysen G. The jasmonate pathway is a key player in systemically induced defense against root knot nematodes in rice. Plant Physiol. 2011;157(1):305–16.
Hamamouch N, Li C, Seo PJ, Park CM, Davis EL. Expression of Arabidopsis pathogenesis-related genes during nematode infection. Mol Plant Pathol. 2011;12(4):355–64.
Wang S, Sun Z, Wang H, Liu L, Lu F, Yang J, et al. Rice OsFLS2-Mediated Perception of Bacterial Flagellins Is Evaded by Xanthomonas oryzae pvs. oryzae and oryzicola. Mol Plant. 2015;8(7):1024–37.
Lu F, Wang H, Wang S, Jiang W, Shan C, Li B, et al. Enhancement of innate immune system in monocot rice by transferring the dicotyledonous elongation factor Tu receptor EFR. J Integr Plant Biol. 2015;57(7):641–52.
Klink VP, Overall CC, Alkharouf NW, MacDonald MH, Matthews BF. Laser capture microdissection (LCM) and comparative microarray expression analysis of syncytial cells isolated from incompatible and compatible soybean (Glycine max) roots infected by the soybean cyst nematode (Heterodera glycines). Planta. 2007;226(6):1389–409.
Bhattarai KK, Atamian HS, Kaloshian I, Eulgem T. WRKY72-type transcription factors contribute to basal immunity in tomato and Arabidopsis as well as gene-for-gene resistance mediated by the tomato R gene Mi-1. Plant J. 2010;63(2):229–40.
Boyes DC, Nam J, Dangl JL. The Arabidopsis thaliana RPM1 disease resistance gene product is a peripheral plasma membrane protein that is degraded coincident with the hypersensitive response. Proc Natl Acad Sci U S A. 1998;95(26):15849–54.
Tornero P, Chao RA, Luthin WN, Goff SA, Dangl JL. Large-scale structure-function analysis of the Arabidopsis RPM1 disease resistance protein. Plant Cell. 2002;14(2):435–50.
Bhattarai KK, Li Q, Liu Y, Dinesh-Kumar SP, Kaloshian I. The Mi-1-mediated pest resistance requires Hsp90 and Sgt1. Plant Physiol. 2007;144(1):312–23.
Scofield SR, Huang L, Brandt AS, Gill BS. Development of a virus-induced gene-silencing system for hexaploid wheat and its use in functional analysis of the Lr21-mediated leaf rust resistance pathway. Plant Physiol. 2005;138(4):2165–73.
Liu Y, Burch-Smith T, Schiff M, Feng S, Dinesh-Kumar SP. Molecular chaperone Hsp90 associates with resistance protein N and its signaling proteins SGT1 and Rar1 to modulate an innate immune response in plants. J Biol Chem. 2004;279(3):2101–8.
Takahashi A, Casais C, Ichimura K, Shirasu K. HSP90 interacts with RAR1 and SGT1 and is essential for RPS2-mediated disease resistance in Arabidopsis. Proc Natl Acad Sci U S A. 2003;100(20):11777–82.
Long HB, Peng H, Huang WK, Wang GF, Gao BL, Moens M, et al. Identification and molecular characterization of a new β-1,4-endoglucanase gene (Ha-eng-1a) in the cereal cyst nematode Heterodera avenae. Eur J Plant Pathol. 2012;134:391–400.
Thies JA, Merrill SB, Corley EL. Red food coloring stain: new, safer procedures for staining nematodes in roots and egg masses on root surfaces. J Nematol. 2002;34(2):179–81.
Li R, Yu C, Li Y, Lam TW, Yiu SM, Kristiansen K, et al. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics. 2009;25(15):1966–7.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.
Iseli C, Jongeneel CV, Bucher P. ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences. Proc Int Conf Intell Syst Mol Biol. 1999:138–148.
Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.
Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36(Database issue):D480–4.
Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5(7):621–8.
Audic S, Claverie JM. The significance of digital gene expression profiles. Genome Res. 1997;7(10):986–95.
Park CH, Chen S, Shirsekar G, Zhou B, Khang CH, Songkumarn P, et al. The Magnaporthe oryzae effector AvrPiz-t targets the RING E3 ubiquitin ligase APIP6 to suppress pathogen-associated molecular pattern-triggered immunity in rice. Plant Cell. 2012;24(11):4748–62.
Wang GF, Peng DL, Gao BL, Huang WK, Kong LA, Long HB, et al. Comparative transcriptome analysis of two races of Heterodera glycines at different developmental stages. PLoS ONE. 2014;9(3):e91634.
Simonetti E, Veronico P, Melillo MT, Delibes A, Andres MF, Lopez-Brana I. Analysis of class III peroxidase genes expressed in roots of resistant and susceptible wheat lines infected by Heterodera avenae. Mol Plant Microbe Interact. 2009;22(9):1081–92.
This project was supported by grants from the 973 project (2013CB127502) and the Natural Science Foundation of China (31301645,31171827).
The authors declare that they have no competing interests.
KLA carried out the biological and molecular studies, participated in the transcriptomic analysis and drafted the manuscript. WDQ and CJK participated in the biological studies. HWK, PH and WGF participated in the design of the study. LZG and YJ participated in transcriptomic data analysis. LSM helped to draft the manuscript. KLA and PDL conceived of the study, and participated in its design and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.
VP1620 disturbed developmental transitions of H. avenae. H. avenae at developmental post-J2, J3 and J4 stages within VP1620 and WEN19 were examined by staining with acid fuchsin at 14 d, 19 d, 25 d, and 33 d. The developmental transitions within VP1620 were significantly affected compared to those in WEN19 (scale bar = 500 μm). (PPTX 540 kb)
Number and percentage of H. avenae juveniles in wheat roots of both compatible and incompatible wheat lines at different time-course points. (DOCX 16 kb)
Identification of dramatically upregulated genes within VP1620 at three time points. Venn diagram showing the numbers of genes significantly upregulated at 24 h, 3 d and 8 d, respectively. Overlapping numbers indicate the amount of overlapping genes simultaneously and significantly upregulated among two or three time points. (PPTX 42 kb)
Annotation of 18 transcripts with dramatically and simultaneously induction among three time points. (DOCX 17 kb)
KEGG pathways with candidate resistance genes. (DOCX 21 kb)
Identification of ROS-producing genes and primers used for qRT-PCR assay. (XLSX 11 kb)
Additional rescued candidate resistance genes. (XLSX 11 kb)
KEGG enrichment analysis of the rescued genes exhibiting uncoordinated expression profiles. (XLSX 16 kb)
Expression profiles of VP1620-specific genes. Cross-table representation of the expression profiles of VP1620-unique genes within the I_CN and I_0 groups. Numbers representing significantly upregulated genes induced by CCN infestation in each cluster combination within the dataset are indicated in each square. The top and left squares show trend lines for changes in expression pattern across three time points. Abbreviation: U, up-forward trend; D, down-forward trend; F, unchanged trend. a Total numbers of genes displaying expression profiles of each I_CN combined with nine I_0; b Total numbers of genes displaying similar expression profiles on the diagonal line; c Total numbers of genes displaying uncoordinated expression profiles outside the diagonal line. (PPTX 101 kb)
About this article
Cite this article
Kong, LA., Wu, DQ., Huang, WK. et al. Large-scale identification of wheat genes resistant to cereal cyst nematode Heterodera avenae using comparative transcriptomic analysis. BMC Genomics 16, 801 (2015). https://doi.org/10.1186/s12864-015-2037-8
- Cereal cyst nematode
- Heterodera avenae
- Triticum aestivum
- Transcripts filtration
- Expression profile
- Resistance genes
- Defense response