Comparative Transcriptome Analysis Uncovers Genes and Pathways Relating to Downy Mildew Resistance in Isabgol (Plantago ovata Forsk.)


 BackgroundIsabgol (Plantago ovata Forsk.) downy mildew (DM), caused by obligate oomycete pathogen Peronospora plantaginis Underwood., is the single most damaging disease of Isabgol. However, reports on the genes and pathways involved in mediating resistance are still unknown.ResultsIn the present study, transcriptomic analysis was carried out by next generation sequencing of DM infected and uninfected leaves of resistant and susceptible genotypes to understand the genetic mechanisms underlie the host-plant resistance. A total of 33.47 million reads were generated and de novo assembled and annotated. The assembly using Trinity yielded 38803, 40175, 45451 and 39533 non-redundant transcript contigs, respectively in susceptible infected, Susceptible uninfected, Resistant infected and resistant uninfected samples. More than 90% of coding DNA sequence (CDS) predicted were annotated using BLAST search. Isabgol transcriptome showed the highest similarity to Sesamum indicum (41-59%) followed by Erythranthe guttata (10-13%). Putative CDS encoding Pathogen associated Molecular Pattern (PAMP)-triggered immunity (PTI), Effector-Triggered Immunity (ETI), Cell wall degrading enzymes, Phytohormone signalling and Phenylpropanoid biosynthesis pathways involved in host-pathogen interaction were identified in addition to the identification of several candidate resistance (R) genes enriched in response to DM infection. We identified significantly differentially expressed genes (DEGs) genes; 6928 DEGs in resistant uninfected (RU) vs resistant infected (RI) of DPO-185 (resistant) genotype and 8779 susceptible uninfected (SU) vs susceptible infected (SI) of DPO-14 (susceptible) genotype. Expression of 11 genes involved in plant defense quantified accordingly using RT-qPCR.ConclusionsThis study for the first time provides the glimpse of transcriptional responses to the DM resistance in Isabgol. The genes and pathways will be helpful in further investigating the molecular mechanisms associated with DM disease resistance in Isabgol and to develop control mechanisms accordingly.

ubiquitin mediated proteolysis and polyprionid biosynthesis involved in the defence response to the pathogens in plants, too [10]. Defense responses include strengthening the cell walls [11], the synthesis of pathogenesis related (PR) proteins and antimicrobial compounds such as phytoalexins [12], and the hypersensitive response (HR), in which cells undergo programmed cell death in the infected region to block further spreading of the pathogen [13].
The oomycete pathogen Peronospora plantaginis (PP) causes Downy mildew (DM) disease of Isabgol is a typical obligate biotroph [14][15]. The pathogen releases motile zoospores on the host plant and enters leaves through stomata [16]. Like other DM pathogens, the colonization of PP fungi is achieved by germinating zoospores, intercellar mycelia growth in the host cells and formation of haustoria. DM of isabgol produces symptoms in leaves as well as floral parts [17]. Disease affected leaves show characteristic ash-coloured downy growth, and as the disease progresses, leaves turn yellowish due to loss of chlorophyll, ultimately reducing effective photosynthetic area [14]. DM outbreak coincides with flowering time and spreads considerably during anthesis and grain filling phases [18]. Floral infection includes systemic infection resulted in long spikes bearing weak and sterile florets, which later turned black due to saprophytic growth [15]. DM disease adversely affects seed yield (73.45% loss) and quality of Isabgol [15].
Controlling the disease through chemical sprays not economical to the farmers besides being environmentally unsustainable. Therefore, development DM resistant varieties is one of the objectives of Isabgol breeding, especially in the warm and humid Isabgol-growing regions where the DM disease occurs regularly. It has also been recognized as one of the most economical, environmentally safe and effective strategies for disease control [19]. The available varieties of Isabgol are susceptible to DM, with poor to moderate yielders. Development of DM resistant varieties becomes a challenge for Isabgol breeders. Characterization of the mechanism underlying DM resistance and the breeding new varieties with resistance and good agronomic traits have been very difficult because resistance to DM is a complex trait and controlled by polygenes [20][21]. Furthermore, current understanding of DM resistance mechanism in Isabgol at the enzyme, gene and regulatory levels is incomplete. Recently, genome-wide studies of model plant species have resulted in an increase in our knowledge of, and capacity to understand, basic biological processes. More specifically, transcriptome analysis has great promise in identification of various genes involved in plant resistance. Transcriptome analysis gives the gene index of an organism quickly, is used to identify the genes involved in specific resistance pathways [22]. Therefore, transcriptome analysis has been successfully used to investigate the specific resistance genes and pathways against DM of grapes [23][24][25][26][27][28][29] and Arabidopsis [30][31], suggesting that plant defense in response to DM infection is regulated through a complex network involving a large number of genes. Many candidate genes such as RPP, PMR4, RabhA4c, RPW8.1, RPW8.2, PRs, NPRs and WRKY shown to play an important role in resistance and response to DM in Arabidopsis [32][33].
With respect to Isabgol, only a few studies are on the transcriptomic resources: [34] characterized the mucilage biosynthesis pathways in P. ovata using 46,955 unique transcripts, [35] characterized genes involved in xylan biosynthesis using 1,099158 ESTs and [36] generated 250 Mb whole genome sequence to develop simple sequence repeat (SSR) markers for breeding applications. [37][38] sequenced 5,900 ESTs and 3,247 independent mRNAs from the leaf vascular tissue of Plantago major, a relative of Plantago ovata and no complete transcriptome analysis is available in the public domain to understand DM resistance in Isabgol.
Thus, the major objective of the present study is to uncover the genes and pathway relating to hostplant resistance against fungal pathogen P. plantaginis responsible for DM disease of Isabgol by using RNA sequencing technologies. Gene expression profiles were compared between a DM resistant (DPO-185) and DM susceptible (DPO-14) genotypes of Isabgol. Thus, our study will lead to better understanding of the molecular mechanisms regulating DM resistance in Isabgol.

Resistance and susceptible reaction of genotypes
Two Isabgol genotypes DPO-14 and DPO-185 were evaluvated for DM resistance during 2014-2017 by artificial inoculation of PP spores under field condition. The genotype DPO-14 recorded mean score of 4.8 with the score ranging from 4.5 to 5.0 across seasons (Fig. 1). While the genotype DPO-185 recorded mean score of 1.1 with score ranging from 1.0 to 1.5 across seasons indicating the resistance or susceptibility of genotypes to the incidence of DM disease. DPO-14 has most prominent ash-coloured downy growth on the leaves, leaves turn yellowish at the time of flowering, ultimately leading to plant mortality as the disease progresses. A sporadic presence of ash-coloured downy growth on the leaves were found in DPO-185 at the time of flowering.    (Fig. 3). In the molecular function category, "catalytic activity" (GO: 0003824) and "binding activity" (GO: 0005488) were most abundantly represented. Under the cellular component category, the highest number of CDS were associated with "membrane" (GO:0019898) and "cell" (GO: GO:0005623).

KEGG mapping
Ortholog assignment and mapping of the CDS to the biological pathways were performed using KEGG automatic annotation server (KAAS) with default score. A total of 6247, 5266, 5016 and 6278 CDS were enriched into 41 different functional KASS pathway categories respectively in SI, SU, RI and RU (Additional Table 1). The mapped CDS represented the genes involved in metabolism, genetic information processing, environmental information processing, cellular processes organizational systems and human disease. Metabolic pathways of major biomolecules such as carbon, carbohydrates, lipids, nucleotides, amino acids, glycans, cofactors, vitamins, terpenoids, polyketides and others were mapped in the transcriptome of Isabgol.
To further predict the function, CDS were subjected to classification into different protein families We used InterProscan to see protein similarity at domain level, in total, 25274, 15652, 14134 and 21828 transcripts in SI, SU, RI and RU respectively were annotated against the Pfam domains ( Table   2). Pkinase (PF00069) domains represented the most which was followed by PPR (PF01535) and Pkinase_Tyr (PF07714) domains in the Isabgol leaf transcriptome indicating strong signal transduction mechanisms. (Additional Fig. S1, Additional file S1 deep analysis of these genes may shed light on the resistance mechanism in Isabgol. However, some of the changes in gene expression may be due to the difference in genetic background of these genotypes. Furthermore, the functional classification of DEGs in leaf transcriptome was analysed using GO (Gene ontology). A total of 1519 and 3455 DEGs respectively in RU vs RI and SU vs SI were enriched to at least one GO term ( Fig. 6; Additional file S4). Among the biological processes, where the term occurrences took place, metabolic process (11.06%) followed by oxidation-reduction process (6.30%) (Fig. 6) recorded maximum terms enriched. In addition, DEGs involved in biological processes like autophagy, response to stimulus and immune system were recognised through GO annotation. In cellular components, GO terms enriched were "integral to membrane" (37.29%) followed by "membrane" (7.97%). In molecular function, GO terms enriched were structural constituent of "ribosome" (2.61%), "protein kinase activity" (2.70%) and "protein serine/threonine kinase activity" (2.21%).

Mining genes involved in Host-pathogen interaction pathway
Host plants perceive pathogens by different recognition systems which are accompanied by a set of induced defense to repel pathogen attack. Activation of pattern recognition receptors (PRRs) by pathogen-associated molecular patterns (PAMPs) at the surface of plant elicits a defense programme known as PAMP-triggered immunity (PTI). In plants, there are 25 genes identified to be involved in PTI [6][7]. In the present study, CDS were found for 18

Mining genes involved in Cell wall modifying enzymes
Cell wall forms a dynamic structure that determines the outcome of the interactions between host and pathogens. When pathogens start degrading the plant cell wall components, plants are capable of perceiving the loss of cell wall integrity and subsequently activate the defense signalling pathways.
Pathogens try to escape the plant defense and sometimes take advantage of the host cell wall metabolism to facilitate their entry into the tissue. There were five [Lipid transfer protein 2 (LPT2),

Glutathione S-transferase (GST), Callose synthase (CS), Cinnamyl alcohol dehydrogenase (CAD) and
Lignin-forming anionic peroxidase (LAP)] genes encoding the enzymes involved in modifying the cell wall during host-pathogen interaction. In the present study, CDSs were found for LPT2, GST, CS and CAD and LAP genes in the leaf transcriptome indicating their role in host pathogen interaction (Additional Table 6

Genes involved in Phytohormone signalling
Phytohormones emerged as cellular signal molecules with key functions in the regulation of immune responses to microbial pathogens and insects. Phytohormones, such as abscisic acid (ABA), ethylene (ET), jasmonic acid (JA), salicylic acid (SA), brassinosteroids (BRs), gibberellic acids (GAs), auxins and cytokinins (CK) act as signalling molecules to triggering immune responses to pathogens in plants.
There are 43 genes reported to be involved in phytohormone signalling in plants through various phytohormones (Additional Table 8

Genes involved in phenylpropanoid biosynthetic pathway
Phenylpropanoid biosynthetic pathway is one of the important secondary metabolic pathways which plays important role in plant defense [39]. There were 14 genes involved in phenylpropanoid biosynthetic pathway [39][40]. KEGG annotation of the phenylpropanoid biosynthetic pathway in response to PP infection is presented in Fig. 9. In the present study, CDS were detected for 13 of the 14 genes involved in phenylpropanoid biosynthetic pathway in the leaf transcriptome indicating their role in DM resistance (Additional Table 10

Enriched R-Genes in response to DM infection
According to the gene-for-gene hypothesis, resistance (R) proteins recognize and interact with effectors and trigger ETI [7].  [42]. In the present study, CDS were found for 58 candidate R genes in the leaf transcriptome indicating their role in conferring DM resistance (Additional Table 12

Realtime Gene expression analysis
To validate the gene expression profiles from the high throughput RNA sequencing, the transcript levels of 11 DEGs were selected to profile their expression upon DM infection in resistant (DPO-185) and susceptible (DPO-14) genotypes (Table 5)  Understanding the molecular mechanism will expedite resistance breeding. The cDNA from infected (I) and uninfected (U) leaves of two mutant genotypes were sequenced using Illumina Hi-seq platform.
We obtained 9.738 GB of raw sequence data which is consistent with the characteristic throughput of Illumina Hi-seq. In general, clean and high-quality reads are essential for good assembly and gene expression analysis [34][35][36][37][38][39][40][41][42][43][44][45], and hence pre-processing was performed to remove low complexity, low length, duplicate, and rRNA reads and de novo assembled into contigs and further to transcripts and CDS. Trinity assembler was used in the present study. Trinity is the most extensively used to assembler for illumina sequence data [24][25]45]. Trinity uses a novel method for efficient and robust de novo assembly and combines three independent software modules: Inchworm, Chrysalis, and Butterfly, applied sequentially to process large volumes of RNA-seq reads [46]. N50 is an important measure to describe the quality of assembly. Higher the N 50 better will be the quality of the assembly. The N50 value found in this study was higher than the previous reports in Isabgol [34], Raphanus sativus (773bp) [45], Camelia sinensis (521bp) [47] and Codonopsis pilosula (1,243bp) [48] infers fairly good assembly of transcriptome of the study.  [57]. In the biological process category, most of the CDS were associated with "metabolic processes" followed by "cellular process" which may allow for the identification of novel genes and pathways involved in the secondary metabolite pathways leading to host plant resistance. Earlier, maximum GO terms associated for "metabolic process" was reported in the transcriptome of Asparagus racemosus [57]. Under the molecular function category, the largest number of CDS were grouped in the "catalytic activity" followed by "binding activity" and "transporter activity" indicates the dominance of gene regulation, signal transduction and enzymatically active processes in the cell.
Earlier, maximum GO terms associated for "catalytic activity" was also reported in Glycyrrhiza uralensis transcriptome [58].
Using KEGG mapping of the best hit CDS, we have identified large number of CDS involved in metabolism, genetic information processing, environmental information processing, cellular processes and organizational systems. Metabolic pathways of major biomolecules such as carbon, carbohydrates, lipids, nucleotides, amino acids, glycans, cofactors, vitamins, terpenoids, polyketides and others were mapped in the transcriptome of Isabgol. These pathways lead to the synthesis of several molecules, which have been reported to confer plant signalling [4], defense [23][24][25][26][27] and resistance to DM disease [23][24][25][26][27][30][31]. These pathways were also reported to be involved in the biological activity of Isabgol [34,59]. The COG annotated putative proteins were distributed functionally into 26 protein families. Of these cluster for "Signal transduction mechanisms (T)" represented the largest group, followed by "General function prediction only (R)", "Posttranslational modification, protein turnover, chaperones (O)" in transcriptome of Isabgol. The results indicate that genes of signal transduction mechanisms are activated by DM infection. Earlier, Zhou et al. [60] reported the largest cluster for "Signal transduction mechanisms (T)" in plants. The cell wall determines the outcome of the interactions between plants and pathogens and pathogens need to breach the cell wall to colonize the plant tissue [69]. Alterations of plant cell wall and degradation of the components are common features of infection by pathogens [70][71].
Alteration in the plant cell wall induces defense responses against invading pathogens [72]. Earlier during host-pathogen interaction was reported for plant pathogenic oomycetes [74]. Understanding the genes involved in cell wall degradation during host-pathogen interaction will pave way for developing resistance phenotypes with altered cell wall composition.
Phytohormones play key roles in the complex signalling cascades involved in growth, development and defense responses in plants. Plant defense is triggered by synthesis of phytohormones such as jasmonic acid [75], salicylic acid [76], ethylene [75], brassinosteroids [77], gibberellic acids [78], abscisic acid [79], auxins [78] and cytokinins [80] act as signal molecules. There are 43 genes involved in phytohormone signalling in plants through various phytohormones [78] of which CDS were detected for 30 genes in our study which are important candidate genes to uncover the role of Further studies are needed to validate their functions in DM resistance.
The phenylpropanoid pathway serves as a rich source of metabolites in plants, being required for the biosynthesis of lignin, and serving as a starting point for the production of many other important compounds, such as the flavonoids and coumarins which plays important role in plant defense [82][83]. There were 14 genes involved in phenylpropanoid biosynthetic pathway in plants [84]. In the present study, CDS were detected for 13 of the 14 genes involved in phenylpropanoid biosynthetic pathway in the leaf transcriptome. Several differentially expressed genes involving phenylpropanoid biosynthetic pathway were identified (Fig. 9), suggesting their role in conferring DM resistance. The similar results have been reported in many plant-pathogen interactions [40,83]. The actual role of these genes in conferring the defense mechanism against DM of Isabgol remains to be explored.
Plants have resistance (R) genes whose products mediate resistance to specific microbes or pathogens. The product of R gene is R protein that allows recognition of specific pathogen effectors, either through direct binding or by recognition of the effector's alteration of a host protein [8][9]. The R genes are regarded as the second mechanism of disease resistance in plants [8][9]. R genes occur ubiquitously in the plant kingdom are involved in ETI [41]. More than 70 different R genes showing resistance to major plant pathogens has been isolated and characterized in different crop plants [41].
However, till date, no reports are on R-genes specific to DM resistance in Isabgol. We found, CDSs

Raw data processing and Annotation
Next Generation Sequencing (NGS) using Illumina generated raw data in FASTQ format. Quality of the raw data was validated based on Phred quality score (Q). The raw data was filtered and trimmed for low quality score reads, adaptor and primer sequences were removed and reads less than 40 bp were removed using Trimmomatic v0.30 [87]. Sequencing data with Phred quality score Q ≥20 was further used for assembly using Trinity RNA-Sequence assembler (Version 2013) [46] using optimized parameters and K-mer size set to 25. To evalulate the assembly quality, all the paired end reads were aligned back to these contigs by using Bowtie2 program [88]. In-house Perl scripts were used to assess the length distribution of the transcripts, N 50 number, average length, maximum length and number of contigs in different intervals. Coding DNA sequences (CDS) from assembled transcript contigs were predicted using Transdecoder (rel16JAN2014) software with default parameters [89].
CDS were functionally annotated by aligning to green plant database of NCBI using basic local alignment search tool (BLASTX) [90] with an E-value threshold of 1e -06 . Gene ontology (GO) enrich was calculated using WEGO analysis.
All predicted CDS were annotated against protein database to assign putative function after

Mining genes involved in DM resistance
To identify putative CDS/transcripts involved in disease resistance, the literature survey was made on functional genes, transcription factors and enzymes involved directly or indirectly in plant disease resistance. Orthologous genes were identified in the annotated CDS using local BLAST search using BioEdit software [91]. Putative R genes were mined from PRGdb (www.prgdb.org) by Blast search. The differentially expressed genes were mapped to biological pathways of plant pathogen interaction pathways: PAMP-triggered immunity (PTI) and Effector-triggered immunity (ETI), cell wall modifying enzymes, phytohormone signaling pathways and phenylpropanoid biosynthetic pathway available in the Kyoto Encyclopaedia of Genes and Genomes (KEGG) database.

Identification of differentially expressed genes
The relative expression level of commonly expressed CDS based on common 'nr' blast hit in terms of

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.