A Global Transcriptome and Co-expression Analysis Reveals Robust Host Defense Pathway Reprogramming and Identifies Key Regulators of Early Phases of Cicer-Ascochyta Interactions

Plants are constantly confronted by various pathogens, and the outcome of this battle directly interferes with agricultural production. To combat pathogens, plants have evolved an array of sophisticated mechanisms consisting of constitutive and inducible immune systems (Jones and Dangl 2006). On the other hand, phytopathogens employ diverse strategies for their survival, fitness, and to circumvent plant defense (Glazebrook 2005; Toruño et al. 2016). Pathogen strategies vary with their lifestyle, biotrophs need the host alive and necrotrophs want it to be dead. Evidence supports the existence of a short biotrophic phase during necrotrophic infection as well (Chowdhury et al. 2017; Spanu 2012). The host immunity modulation by secreted effector proteins further challenges our understanding of plant responses to pathogens. Therefore, understanding the molecular basis of multifaceted plant-pathogen interactions remains prime focus to avert the crop loss.

The multilayered plant immune system consists of preformed barriers, pathogen detection, and signal transduction components that coordinate the defense responses (Corwin and Kliebenstein 2017; Tsuda and Katagiri 2010). These pathogen-induced responses involve dynamic transcriptional reprogramming of multiple host processes leading to resistance or susceptibility. For instance, extensive degradation of cell-wall polysaccharides is recognized by various receptor kinases, which, in turn, trigger defense responses, such as cell-wall reinforcement at the pathogen penetration site, by lignin and callose deposition (Brutus et al. 2010; Eynck et al. 2012; Flors et al. 2008). These responses activate reactive oxygen species (ROS) production, pathogenesis-related (PR) gene expression, as well as the accumulation of plant defense hormones like salicylic acid (SA) and jasmonic acid (JA) (de Araújo et al. 2019). SA and JA are pivotal to combat biotic stress; however, emerging evidence indicates antagonistic or synergistic roles of abscisic acid (ABA), ethylene (ET), auxin (Aux), and brassinosteroid (BR) in immunity (Catinot et al. 2015; Glazebrook 2005; Pieterse et al. 2012; Verma et al. 2016a).

Chickpea (Cicer arietinum L.) is one of the most widely grown legumes across the world, especially in arid and semi-arid climates (Singh et al. 2021a). A necrotrophic fungus, Ascochyta rabiei, that causes Ascochyta blight (AB) disease is a major constraint for global chickpea production. It infects all the aerial parts of the plant and causes up to 100% crop loss under conducive conditions (Knights and Hobson 2004; Sharma and Ghosh 2016). Further, the existence of a teleomorph stage in A. rabiei accelerates their evolution, thus, actively defying the plant defense strategies (Mehmood et al. 2017). Consequently, A. rabiei isolates have developed resistance against quinone outside inhibitor (QoI) fungicides in Montana and North Dakota (Delgado et al. 2013; Owati et al. 2017). Importantly, resistance erosion has also been observed in several formerly known ‘AB-resistant’ chickpea cultivars (Sambasivam et al. 2020). This necessitates the identification of key genetic regulators of AB resistance in chickpea, which can be exploited for developing resistant plants. Previous expression studies provided limited understanding of chickpea response against A. rabiei infection (Coram and Pang 2006; Jaiswal et al. 2012). An integrated transcriptome and degradome sequencing analysis uncovered numerous microRNAs and their targets that coordinate AB resistance in chickpea cultivars (Garg et al. 2019). However, the early-stage host transcriptional modulation that decides the fate of the battle during Cicer-Ascochyta interactions remains poorly understood. A. rabiei starts penetrating host tissue at 24 hpi, followed by appearance of yellow spots on the leaves at 72 hpi (Ilarslan and Dolar 2002). These stages are crucial to understand how AB-resistant and -susceptible plants respond to A. rabiei perception and penetration. Hence, in-depth investigations are required to appropriately comprehend the molecular components involved in resistance- and susceptibility-related processes.

In this study, comparative transcriptome analysis of AB-resistant (FLIP84-92C) and -susceptible (PI359075) chickpea genotypes was performed at the early stages (24 and 72 hpi) of A. rabiei infection. We identified differential gene expression patterns for defense-related components, such as receptor kinases, phytohormone pathways, and transcription factors (TFs), between both genotypes. In addition, the incorporation of weighted gene co-expression network analysis (WGCNA) revealed resistant or susceptible genotype–specific gene expression modules and identified hub genes governing contrasting traits. Overall, we provide new insights into the complex Cicer-Ascochyta interactions, delineate AB-responsive host-physiological pathways, and deliver potential resistance genes for further characterization and downstream applications in chickpea improvement programs.

RESULTS

Transcript profiling of AB-resistant and -susceptible chickpea genotypes at early stages of Ascochyta infection.

To understand the host transcriptional dynamics upon A. rabiei infection, transcriptome profiling of resistant (FLIP84-92C) and susceptible (PI359075) chickpea genotypes under control and AB-stress conditions was performed at 24 and 72 hpi. The details of raw, filtered, and mapped reads on the reference genome are given in Supplementary Dataset S1. To identify the differentially expressed genes (DEGs), seven pair-wise comparisons were performed between control, 24-, and 72-hpi timepoints and genotypes. Among them, four sets represent transcriptional reprogramming within the genotype in response to A. rabiei (CPIvs24PI, CFLvs24FL, CPIvs72PI, CFLvs72FL; samples are designated as ‘condition genotype’ [CPI]). For instance, in the CPIvs24PI dataset, CPI represents a PI359075 genotype sample at control condition and 24PI depicts PI359075 genotype sample at 24 hpi. Three sets (CPIvsCFL, 24PIvs24FL, 72PIvs72FL) show differential responses between the susceptible and resistant genotypes. Overall fold change distribution of identified genes is shown in Supplementary Fig. S1.

In response to AB stress, a total of 332 (177 up- and 155 downregulated) and 980 (625 up- and 355 downregulated) genes were differentially expressed in susceptible plants at 24 and 72 hpi, respectively. However, in the resistant genotype, 441 (203 up- and 238 downregulated) and 8,908 (4602 up- and 4,306 downregulated) genes showed differential expression (Fig. 1A). Both genotypes showed more DEGs at 72 hpi, though the number was much higher in the resistant genotype, suggesting rapid and dramatic transcriptional changes. Although, 3.88% (24 hpi) and 5.73% (72 hpi) AB stress–specific DEGs were shared between the resistant and susceptible genotype, their expression levels were not correlated (r < 0.4). Among them, 39 showed an opposite expression pattern in both the genotypes (Supplementary Fig. S2). In contrast, a total of 93.1 and 93.6% DEGs were uniquely expressed in the resistant genotype, while 90.9 and 42.1% DEGs were exclusive to the susceptible genotype at 24 and 72 hpi, respectively. Further, the shared DEGs between 24 and 72 hpi in resistant (3.21%) and susceptible (10.53%) genotypes represent constitutively activated defense mechanisms during various infection stages (Fig. 1B).

Download as PowerPoint

Similarly, when resistant and susceptible genotypes were compared, the number of DEGs varied from 856 (450 up- and 406 downregulated; 24 hpi) to 5,102 (2,672 up- and 2,430 downregulated; 72 hpi) (Fig. 1A). Interestingly, only 10.9% overlapped between the control and AB-stress conditions, thus showing a major proportion of AB-induced DEGs. Moreover, 4.4 and 63% of the DEGs were exclusive to 24 and 72 hpi, respectively (Fig. 1C). DEGs with specific induction or suppression upon infection could mediate the variable defense responses exhibited by these contrasting genotypes. Overall, results illustrated robust transcriptional reprogramming at early stages of infection in chickpea; however, the host response differs between incompatible (resistant) and compatible (susceptible) interactions.

Differential transcriptional modulation of host defense pathways upon pathogen attack.

To gain insights into A. rabiei–triggered pathways in these chickpea genotypes, functional classification and pathway assignment of DEGs were performed. In both genotypes, the upregulated DEGs were mostly associated with defense, response to biotic stimulus, hormone signaling, oxidation-reduction, protein phosphorylation, and modification processes (Supplementary Dataset S2). However, the total number of DEGs present in each of these gene ontology (GO) categories was higher in the resistant genotype compared with susceptible. This suggests the robust activation of defense signaling pathways in the resistant genotype. Upregulated genes in FLIP84-92C were mainly associated with defense response by cell-wall thickening, callose deposition, and regulation of fatty acid oxidation, while developmental and reproductive process were down-regulated. In contrast, upregulated genes in the PI359075 genotype were mainly involved in toxin, phytoalexin, camalexin, and flavonoid biosynthetic processes. Notably, when the resistant genotype was compared with susceptible, ‘photosynthesis’, photosystem II assembly’, ‘photosynthetic electron transport chain’, and ‘photosystem II oxygen evolving complex assembly’ processes showed upregulation in resistant plants, indicating differential regulation of photosynthetic processes between both the genotypes. Moreover, enrichment of resistant genotype-specific DEGs suggests their association with translation, ribosome biogenesis, and lipid biosynthesis processes; however, susceptible genotype-specific DEGs belong to cellular carbohydrate metabolic, small molecule, and vitamin biosynthetic process (Supplementary Fig. S3). Such variation in gene regulation upon Ascochyta infection reflects the different strategies tailored by these genotypes.

Further, categorization of DEGs into a MapMan biotic stress pathway revealed minimal transcriptional alterations at 24 hpi in both the genotypes compared with 72 hpi, though a relatively large number of genes was differentially regulated in the resistant genotype (Supplementary Fig. S4). When the resistant genotype was compared with susceptible, considerable upregulation of redox state, glutathione-S-transferase (GST), signaling, secondary metabolites, PR proteins, and proteolysis pathways was observed at 24 hpi. Consistently, signaling and proteolysis pathways were differentially regulated at 72 hpi but showed slight changes in susceptible genotype (Supplementary Fig. S5A and B). Peroxidases and β-glucanase were mostly down-regulated. Mitogen-activated protein kinases (MAPKs) along with WRKY and MYB TFs were also up-regulated in both the genotypes, although the number of genes was substantially higher in resistant plants. ET-responsive factor (ERF) and DOF (DNA binding with one finger) TFs were predominantly down-regulated but bZIP (basic leucine zipper) TFs were induced in resistant plants, whereas no obvious changes were observed in the susceptible genotype. In addition, genes associated with SA, JA, ET, and ABA were up-regulated in the resistant genotype, while BR and Aux-related genes mainly showed downregulation. On the other hand, only a few hormone-related genes were differentially expressed in the susceptible genotype. These observations are concordant with the enrichment analysis (Supplementary Dataset S2).

Moreover, various abiotic stress–related DEGs were also identified that could be common mediators of biotic and abiotic stress responses. Intriguingly, the resistant genotype showed similar differentially regulated pathways in both the comparisons, i.e., CFLvs72FL and 72PIvs72FL (Supplementary Figs. S5A and S6A), indicating the minimal response of the susceptible genotype against A. rabiei. Differential expression of JA-, ET-, and ABA-related genes along with WRKY, ERF, and MYB TFs was observed in genotype-level comparisons under the control condition, implying pre-activated basal defense mechanisms in the resistant genotype (Supplementary Fig. S6B). Further, these RNA-seq data were validated by quantitative real time-PCR (qRT-PCR) analysis (Supplementary Fig. S7; Supplementary Dataset S3). Overall, these analyses demonstrate robust activation of diverse defense mechanisms in the resistant genotype and highlight the importance of signal transduction and gene regulatory machinery in coordinating high-amplitude transcriptional reprogramming against pathogen.

Protein kinases (PKs) and TFs are extensively implicated in host defense.

PKs and TFs are important players in perception and transduction of stress signals for coordinated defense. In our data, 520 differentially expressed protein kinases (DEPKs), and 905 differentially expressed TFs or regulators (DETFs or DETRs) were identified, which were most abundant in the resistant genotype at 72 hpi (Fig. 2A; Supplementary Dataset S4). Among them, 35 DETF and DETR and 63 DEPK families and subfamilies were uniquely expressed in the resistant genotype upon AB stress (Supplementary Figs. S8A and S9A). Between-genotype comparisons showed 17 unique AB stress–specific DETF and DETR families along with 37 DEPKs in the resistant genotype (Supplementary Figs. S8B and S9B). These families include LIM, MADS-M type, C2C2-LSD, C2C2-CO-like, E2F-DB, and double-B box TFs along with LUG, SWI/SNF-SW13, and Pseudo-ARR-B TRs. PKs include receptor-like kinase (RLK)-Pelle_LRR (leucine-rich repeat) subfamilies V, IV, Xb1, and XIIIa, CMGC_GSK, and TTK.

Download as PowerPoint

Further analysis showed significant enrichment of 20 TF, eight TR (transcription regulator), and 12 PK families, at least, in one of the datasets (Fig. 2B to D). Major TRs (Aux/IAA [indole 3-acetic acid], TRAF, SNF2, Pseudo-ARR-B, HMG) and PKs (RLK-Pelle_LRR, CMGC_MAPK, CMGC_CDK) were enriched in the resistant genotype. WRKY, NAC (NAM, ATAK, and CUC), PLAZ, and FAR1 TFs were enriched in both the genotypes under AB stress, whereas MYB, B3, TCP, AP2/ERFs, HSF, GARP-G2, MADS-MIKC, B3, C2C2-DOF, homeodomain (HD)-HD-ZIP, and zf-HD families were exclusively enriched in the resistant genotype. Moreover, hierarchical clustering identified six expression trend clusters of DETFs, DETRs, and DEPKs (Fig. 2E). Cluster 1 shows upregulated expression in the resistant genotype at 72 hpi, which mainly contains WRKY and B3 TFs, and DLSV (DUF26, SD-1, LRR-VIII, and VWA), LRR, receptor-like cytoplasmic kinase, and wall-associated kinase PKs. Lr10 was the most dominant member in this cluster. Resistant genotype-specific cluster 2 includes downregulated Aux/IAA, bHLH, MYB, HB-HD-ZIP TFs and RLK-Pelle_LRR PKs. Cluster 3 is specific to 72 hpi, showing induced expression in both the genotypes. This cluster harbors bHLH and NAC TFs along with RLK-Pelle_DLSV PKs. Interestingly, cluster 4 demonstrates downregulation at 72 hpi but upregulation under control conditions in the resistant genotype (Fig. 2E). In this cluster, bHLH, GATA, and AP2/ERF-ERF TFs along with RLK-Pelle_LRR and serine/threonine kinase were present. Clusters 5 and 6 represent susceptible and resistant genotype-specific upregulation, respectively. In cluster 5, mainly bHLH, AP2/ERF-ERF, C2H2 TFs and probable inactive receptor kinases were identified, whereas cluster 6 contains a large number of WRKY TFs and RLK-Pelle PKs (Supplementary Dataset S5). Noticeably, most of the TFs and PKs were exclusive to one genotype, while others showed similar expression patterns in both genotypes, implicating both common and specific responses to Ascochyta infection in resistant and susceptible hosts. As PKs and TFs regulate multiple signaling pathways, the high abundance of DEPKs and DETFs between both the genotypes might be the main driving force for the differential expression of downstream defense-related genes leading to the distinct response of resistant and susceptible genotypes against AB stress.

The gene expression dynamics are tightly controlled by TFs, which bind to specific cis-regulatory regions of target genes. In our data, targets of MADS, ERF, MYB, NAC, and WRKY TFs were over-represented upon infection, suggesting their role in plant immunity (Supplementary Dataset S6). MYB14 and MYB36 TFs enriched in our datasets (CFLvs72FL and 72PIvs72FL) target phytohormone-related genes, thus suggesting the intricate role of hormones during Cicer-Ascochyta interactions. Along with PR gene transcription activator (PTI6), targets for MADS-box protein, SUPPRESOR OF OVEREXPRESSION OF CONSTANS1 (SOC1) and BASIC-PENTACYSTEINE (BPC), were markedly enriched among the DEGs. The SOC1 and BPC regulate flower and leaf development-related genes in Arabidopsis (Lee et al. 2008; Shanks et al. 2018). However, defense-related targets of these TFs including disease resistance protein RPS5, thaumatin-like protein PR5a, PR10, and WRKY TFs along with hormone-related genes were also abundant in our datasets, indicating their potential role in plant immunity, which remains to be explored. Interestingly, growth and development pathways were down-regulated in the resistant genotype upon A. rabiei exposure (Supplementary Dataset S2). Thus, these genes might be crucial players in the growth-immunity trade-off.

Modulation of hormonal pathway genes highlights the central role of JA, ET, and ABA signaling in response to blight fungus.

The phytohormone crosstalk plays a significant role in defense against invading pathogens. In our datasets, we identified seven, 33, 34, 41, 54, and 90 DEGs associated with SA, ABA, BR, JA, ET, and Aux pathways, respectively (Supplementary Dataset S7). JA biosynthesis genes such as lipoxygenase (LOX), allene oxide synthase (AOS), allene oxide cyclase (AOC), and 12-oxophytodienoate reductase (OPR) were differentially expressed in the resistant genotype (Fig. 3A). LOX and OPRs were up-regulated, whereas AOS and AOC were down-regulated upon A. rabiei infection in the resistant genotype. Moreover, 1-aminocyclopropane-1-carboxylate (ACC) synthase and ACC oxidase were up-regulated, suggesting induced biosynthesis of ET (Fig. 3B). In ABA biosynthesis, zeaxanthin epoxidase was down-regulated, while short-chain alcohol dehydrogenase/reductase (SDR) and abscisic aldehyde oxidase were induced in both the genotypes (Fig. 3C). SA, BR, and Aux biosynthesis genes were not differentially expressed. Furthermore, ABA-8′-hydroxylase was down-regulated in both the genotypes, while jasmonate methyltransferase (JMT) and salicylate carboxymethyltransferase (SAMT) were exclusively induced in the resistant genotype at 72 hpi (Fig. 3D). This suggests the activated biosynthesis of methyl jasmonate (MeJA) and SA-methyl ester (MeSA) for JA- and SA-mediated signaling during plant defense.

Download as PowerPoint

Moreover, DEGs associated with JA, SA, ABA, and ET regulatory pathways were up-regulated, whereas Aux and BR were down-regulated in the resistant genotype at 72 hpi (Supplementary Fig. S10; Supplementary Dataset S7). At 24 hpi, fewer changes were observed. SA response-associated DEGs such as methylxanthosine synthase, pentatricopeptide repeat-containing protein, and glycosyltransferase were exclusively up-regulated at 72 hpi in the resistant genotype. Intriguingly, many DEGs were shared in SA and JA pathways, which is suggestive of possible crosstalk. In the resistant genotype, ERFs were differentially expressed, whereas CBS domain–containing protein, ET receptor, and codine-O-demethylase were clearly overexpressed. In addition, ABA-responsive genes like those of GEM protein, HVA22, and serine/threonine phosphatases were also up-regulated, whereas, Aux-responsive (SAUR and Aux/IAA) and BR-related genes were down-regulated in resistant plants. Further, Aux transporters were differentially regulated, indicating the altered Aux levels during fungal infection. The elevated Aux levels are correlated with Fusarium head blight susceptibility in wheat (Wang et al. 2018).

These observations were also supported by the GO enrichment analysis. In the susceptible genotype, ‘response to JA’ and ‘response to SA’ terms were over-represented among upregulated DEGs, while ‘response to JA’, ‘response to ABA’, ‘SA-mediated signaling pathway’, ‘ET-activated signaling pathway’, and ‘response to ET’ were enriched in the resistant genotype. Among the hormone-responsive terms, mainly WRKY40, MYB33, MYB48, ARR2 TFs, PR proteins, and class I and V chitinases were identified in the resistant genotype. On the other hand, WRKY43, NAC05, and MYB48 TFs were present in the susceptible genotype. Downregulated terms include ‘Aux-activated signaling pathway’ and ‘Aux polar transport’ in the resistant genotype, whereas ‘response to BR’ was commonly represented in both the genotypes (Supplementary Dataset S2). Among the Aux and BR response terms, only AP2/ERF-AIL5 and NAC32 TFs were present, respectively. Identification of these TFs and defense genes under the hormone response category implicates their role in downstream activation of host processes. Collectively, our findings demonstrate induced biosynthesis and signaling of JA, ET, and ABA in the resistant genotype. The SA, Aux, and BR biosynthesis genes remain unchanged. Moreover, the regulatory genes of SA showed upregulation, while Aux and BR were down-regulated during A. rabiei infection. This implicates JA, ET, ABA, and their plausible crosstalk as central regulators of chickpea defense against AB stress. Further, these phytohormones modulate various genes through signaling cascades, which is indicative of complex co-expression programs that work together to defend the host plant.

AB resistance and susceptibility–associated co-expression modules in chickpea.

Further, we performed WGCNA analysis to divulge the modules and gene regulatory networks (GRNs) associated with particular genotypes and timepoints. Correlation of module eigengene with sample trait (i.e., resistance and susceptibility) identified 12 co-expression modules (comprised of 39 to 5,425 genes) associated with particular conditions (Fig. 4A). For further analysis, highly correlated (r ≥ 0.5 and ≤ −0.5) and significant modules (P value ≤0.05) were selected in each dataset. Under control conditions, the strongest positive correlation (r = 0.91) was obtained for the light-green module in the resistant genotype. Functional categorization of this module showed over-representation of ‘defense response to fungus’, ‘response to biotic stimulus’, ‘plant-pathogen interaction’, and ‘MAPK signaling pathways’, signifying the basal pre-formed defense in resistant genotype (Fig. 4B and C).

Download as PowerPoint

Under AB stress, light-cyan1 (r = 0.93) and red (r = 0.96) modules were identified at 24 and 72 hpi, respectively, showing strong positive correlation with the resistance trait. A higher number of genes was identified in the red module (638) compared with light cyan1 (39), which is concordant with our differential gene expression data. The light-cyan1 module over-represents monocarboxylic acid, lipid and fatty acid biosynthetic, fatty acid elongation, cutin, suberin, wax, isoquinoline alkaloid, and secondary metabolites biosynthesis-related terms, suggesting rapid cell-wall remodeling and secondary metabolite production in resistant chickpea to restrict fungal penetration. At the later stage (red module), ‘response to wounding’, ‘flavonoid biosynthetic and metabolic process’, ‘chitin metabolic process’ ‘protein phosphorylation’, ‘biosynthesis of secondary metabolites’, ‘glutathione metabolism’, ‘zeatin biosynthesis’ and ‘plant-pathogen interaction’ terms were over-represented. Moreover, a negatively correlated sky-blue module (r = −0.84) was identified at 72 hpi, harboring a gene-set belonging to ‘callose deposition in phloem sieve plate’, ‘cell wall organization or biogenesis’, and ‘monoterpenoid metabolic and biosynthetic process’ categories in resistant genotype (Fig. 4B and C). The data indicates extensive involvement of post-translational modifications via kinase activity, antioxidant metabolism, and flavonoid biosynthesis pathways in providing resistance to AB disease, which is in agreement with our DEG data.

In the susceptible genotype, a positively correlated (r = 0.77) dark-magenta module was detected at 72 hpi, showing over-representation of ‘polysaccharide biosynthetic process’, ‘plant-type secondary cell wall biogenesis’, ‘alanine, aspartate, and glutamate metabolism’, and ‘phenylpropanoid biosynthesis’ (Fig. 4; Supplementary Dataset S8). Interestingly, the ‘protein phosphorylation’ category was also enriched; however, the number of genes was much higher in the resistant (red; 50) than the susceptible genotype (dark magenta; 32). This signifies the role of PKs in regulating the function of downstream genes mediating AB resistance, which is consistent with DEPK analysis. Noticeably, no significant modules were identified under the control and 24-hpi conditions in the susceptible genotype, which further indicates the delayed or suppressed defense (Fig. 4A). Overall, enrichment analysis presents an overview of the reprogrammed chickpea defense pathways and indicates the presence of genotype-specific hub genes mediating the differential response of both genotypes.

Core regulators orchestrating the host response against AB disease.

To identify the core regulators governing AB resistance and susceptibility, we performed GRN and gene significance and module membership (GS-MM) analysis (Fig. 5). The highly connected hub genes in GRN could be master regulators of the module and, hence, the associated trait. GS-MM analysis identified trait-associated hub genes from light-green (66), light-cyan1 (25), red (356), sky-blue (187), and dark-magenta (53) modules (Supplementary Figs. S11 and S12; Supplementary Dataset S9) and provided crucial insights into the differential response regulators in AB-resistant and susceptible genotypes.

Download as PowerPoint

Basal defense is higher in the AB-resistant genotype.

Basal expression of defense genes has an imperative role in preventing fungal invasion. Notably, no significant module was detected in the susceptible genotype; however, the hub genes identified under a resistant (light-green) module were mainly involved in defense response, immune system process, and plant-pathogen interaction. These genes include 3-ketoacyl-CoA synthase, ERF027 and ABA-8′-hydroxylase. Genotype-level differences under control conditions were also clearly observed in our gene expression analysis, as 2,027 (960 up- and 1,064 downregulated) genes showed differential response between both genotypes (Fig. 1A). Further, GRN revealed calcium-binding protein CML27, ERF9, MYB1, NDR1/HIN113, cysteine-rich, and transmembrane domain–containing protein B as highly connected hub genes, potentially regulating the basal defense mechanisms (Fig. 5A). The relatively high expression of these genes compared with susceptible genotype indicates the pre-activated defense in AB-resistant chickpea genotype.

Cell wall–thickening and secondary metabolite biosynthesis mitigate infection at the early stage.

In the light-cyan1 module, hub genes encoding for fatty alcohol:caffeoyl-CoA acyltransferase, 3-ketoacyl-CoA synthase, alkane hydroxylase MAH1, secoisolariciresinol dehydrogenase, and diacylglycerol O-acyltransferase were associated with cutin, suberin, wax, and lignan biosynthesis. Besides, pleiotropic drug resistance protein 1 (PDR1), tyrosine decarboxylase 1, and PRE6 were also found highly connected in the network (Fig. 5B). PRE6 is an important transcriptional repressor of Aux-response factors (Zheng et al. 2017). Tyrosine decarboxylase is involved in catecholamines biosynthesis (Świędrych et al. 2004). PDR1 participates in both constitutive and JA-induced defense against Botrytis cinerea (Stukkens et al. 2005). Moreover, cytochrome P450 oxidoreductase CYP736A12 and BRO1 domain–containing protein were also present in this module. The BRO1 domain–containing protein is induced in endoplasmic reticulum–stressed rice seeds (Qian et al. 2015). Although, no significant module was identified in the susceptible genotype at 24 hpi, genes related with toxin, phytoalexins, and secondary metabolite biosynthesis were over-represented among DEGs (Supplementary Dataset S2). These data suggest that faster activation of cell-wall thickening and defense-related secondary metabolite biosynthesis is crucial for early host response against A. rabiei.

Signaling, protein phosphorylation, and transcriptional machinery are majorly involved in coordinating AB resistance.

At the later stage (72 hpi), mostly PKs (32) and TFs (44) were over-represented in resistant genotype–associated modules. The modules harbor 28 upregulated RLK-Pelle members including LRR, DLSV, and LRR-like kinases, 12 downregulated receptors that were mainly LRRs, and three differentially expressed MAPKs (MAPK4, MAPK9, and MAPK11) (Supplementary Dataset S9). This indicates extensive involvement of post-translational regulation via protein phosphorylation at current or later stages. The WRKY, NAC, bZIP, and MYB TFs that play an important role in plant immunity were also induced along with other TFs, such as B3, GARP-G2-like, MADS-MIKC, whose functions in defense remain elusive. Downregulation of AP2/ERF, C2C2-Dof, and LIM TFs suggests that they could be negative regulators of AB stress. Further, GRN identified plastid movement-impaired protein, NAD(P)H-dependent 6′-deoxychalcone synthase, MYB14, chitinase, and vestitone reductase hub genes, suggesting that these could be master regulators of overall chickpea response to A. rabiei (Fig. 5C). Intriguingly, OPR11, SAMT, AdoMet synthase, and SDR that are involved in phytohormone pathways were commonly identified between hub and DEGs. PR proteins (PR-1B, PR-4, PR-5a, PR-10), E3 ubiquitin protein ligase DRIP2, thaumatin-like protein, and GSTs, along with calcium signaling–related genes were also activated. In accordance to 24 hpi, coumarate-CoA ligase, feruloyl-CoA orthohydorxylase and isoflavone-2′-hydroxylase, involved in biosynthesis of phenylpropanoids and phytoalexins (Akashi et al. 1998; Kai et al. 2008; Liu et al. 2017; Schmid et al. 2014), were also present among the hub genes at 72 hpi. The enrichment of defense-related genes and over-representation of TFs, based on target abundance in the resistant genotype compared with the susceptible genotype, explains the contrasting behavior of both genotypes during fungal colonization (Supplementary Dataset S10). In the sky-blue network, protein trichome birefringence 36, peroxidase, probable fructokinase, tubulin β-chain, and 2-oxoglutarate-dependent dioxegenase were identified as hub genes (Fig. 5D). In addition, susceptibility-associated genes like pectate lyase, polygalacturanse1, PR1, and Avr9/Cf-9 were identified. Pectate lyase (PMR6) confers powdery mildew susceptibility to Arabidopsis (Vogel et al. 2002), whereas polygalacturanse1 suppresses the programmed cell death (He et al. 2019). Importantly, downregulation of SWEET transporter in the resistant genotype further substantiates its role in disease susceptibility, as observed in an earlier report (Cohn et al. 2014). ENHANCED PSEUDOMONAS SUSCEPTIBILTY1, which provides SA-dependent resistance (Zheng et al. 2009), is down-regulated as SA-mediated defense is suppressed during necrotrophic infections. Intriguingly, development-related hub genes such as FANTASTIC4, B-like cyclin, FT-interacting protein1 were also identified in a downregulated module, implying a growth-immunity trade-off.

Regulators coordinating AB-susceptibility response in chickpea.

The hub genes identified in GRN of the susceptible genotype (dark magenta) were fasciclin-like arabinogalactan protein12, UTP-glucose-1-phosphate uridylyltransferase, cellulose synthase, DUF1336 family protein, and putative aminoacyltransferase (Fig. 5E). Importantly, the role of these genes is underexplored in plant defense and susceptibility. Thus, they could be novel determinants for AB disease. PKs such as CAMK_CAMKL-CHK1/2, RLK-Pelle_LRR-XI-1/2, and bZIP, NAC, NF-YB TFs were identified among the hub genes. In addition, serine/threonine kinases, peroxidase, polyphenol oxidase A1, COBRA, and putative senescence regulator S40 also responded to A. rabiei. The polyphenol oxidase A1 regulates secondary metabolism and cell death in walnut (Araji et al. 2014). Loss-of-function of COBRA invokes JA-regulated defense in Arabidopsis (Ko et al. 2006). However, putative senescence regulator S40 is activated in response to natural senescence, darkness, pathogen, SA, and ABA treatment (Krupinska et al. 2002). Overall, this module contains novel susceptibility factors and also suggests minimal activation of immunity in the susceptible genotype, which was unable to mitigate the disease.

DISCUSSION

To combat invading pathogens, plants modulate their growth, development, and defense via complex and dynamic transcriptional alterations that start with perception and last till the defeat of pathogen or host (i.e., resistance or disease). Hence, the initial stages of infection are critical for the establishment of the pathogen inside the host and disease development (Mine et al. 2018). In chickpea, the early defense signaling and regulatory networks differentially operating in resistant and susceptible cultivars upon AB challenge remain elusive. Here, we investigated the pathogen-induced host transcriptional reprogramming occurring at the early stages of A. rabiei infection in well-characterized AB-resistant (FLIP84-92C) and -susceptible (PI359075) chickpea genotypes, using an RNA-seq approach. Our comprehensive analysis provides an overview of the chickpea defense response to a fungal necrotroph (Fig. 6) and identifies differentially regulated genes and co-expression modules containing receptors, PKs, and TFs, along with other biosynthetic enzymes. In addition, novel hub genes were defined that are potentially responsible for the resistance or susceptibility of chickpea genotypes under AB stress.

Download as PowerPoint

A plethora of proteins participate in basal and induced processes of the multilayered immune system. The pre-formed or constitutively elevated expression of plant defense genes plays a major role in preventing pathogen invasion (Vergne et al. 2010). The wide-spectrum differential gene regulation between AB-resistant and -susceptible chickpea genotypes under control conditions suggests the existence of a functional and efficient pre-formed immune system in the resistant genotype. However, the pre-activated basal defense components in susceptible genotype are not enough to restrict fungal invasion and might serve as susceptibility factors. Extensive cell-wall reinforcement and production of antimicrobial compounds are non-specific plant responses that act as the first barrier against pathogens (Muthappa and Mysore 2013). They prevent or delay fungal penetration, thus buying time for the activation of a downstream defense system. At 24 hpi, mainly lignan, cutin, wax, suberin, and secondary metabolite biosynthesis-related genes were significantly induced, resulting in cell-wall thickening. This data complements the earlier histological observations of physical barriers and cuticle layer in AB-resistant chickpea cultivars (Ilarslan and Dolar 2002). Intriguingly, high-amplitude transcriptional reprogramming in the resistant genotype, especially at 72 hpi, highlighted the importance of this timepoint in deciding the fate of A. rabiei colonization. However, data for later stages would be vital to pinpoint the minimum duration of such a high transcriptional reprogramming. Both genotypes shared a small subset of DEGs at 24 and 72 hpi, which suggests the differential response of these cultivars against AB stress. Enrichment analysis demonstrated an extensive nexus of signaling pathways, post-translational modifications, phytohormone crosstalk, and secondary metabolism in the resistant genotype to overcome pathogen.

The plant surveillance system against pathogens mainly relies on RLKs and NOD-like receptors (Tang et al. 2017; Wang et al. 2020). Compared to the susceptible genotype, high numbers of RLK-Pelle members were identified in the resistant genotype upon AB stress, whereas none of them were induced in the susceptible genotype at 24 hpi. Moreover, rust-resistance kinase Lr10, receptor-like protein EIX1, receptor-like protein kinase HSL1, and wall-associated receptor kinase 20 were identified among hub genes under resistant modules. Lr10 is a coiled coil-nucleotide binding site-LRR receptor involved in immunity against wheat rust disease (Feuillet et al. 2003). MAPK signaling cascades triggered upon pathogen recognition act on WRKY TFs, which are regulators of host defense responses against a variety of pathogens (Eulgem and Somssich 2007). Also, in our datasets, several MAPKs and WRKY TFs, especially WRKY4, WRKY47, WRKY31 and WRKY65 were mainly induced in the FLIP84-92C genotype. Besides, MADS-box, MYB, bHLH, and B3 TFs were also preferentially expressed in resistant genotype. TFs play a crucial role in biotic stress response through complex GRNs (Park et al. 2001; Wang et al. 2009). In our study, MYB14 and PRE6 showed the highest connectivity, hence plausibly modulating other genes of the resistant module. Furthermore, the abundance of SOC1 and BPC1 targets, such as disease resistance protein RPS5, PR genes, and other defense-related transcripts among the DEGs, prospects them as master regulators of chickpea immunity against A. rabiei, apart from their developmental functions. Notably, PR-1B, PR4, PR-5a, and PR-10 were also identified as hub genes in resistant plants. PR-10 provides defense by preventing or disrupting the conidial germination and hyphal growth (Zandvakili et al. 2017). In Arabidopsis enhanced disease susceptibility (eds) mutants, the altered PR1 expression enhances susceptibility to bacterial pathogens (Rogers and Ausubel 1997). Together, results concluded that the resistant cultivar triggered prompt and robust defense via reprogrammed expression of signaling and transcriptional regulatory components in response to A. rabiei. On the contrary, the susceptible genotype could not mount sufficient response to mitigate the disease progression. The minimal activation or suppressed immunity could be attributed to the secreted effector proteins as the Ascochyta genome encodes an armory of predicted effector proteins (Verma et al. 2016b).

SA and JA are classically known to regulate signaling networks involved in induced defense responses against biotrophic and necrotrophic pathogens, respectively (Ryan and Moura 2002; van Loon et al. 2006; Abd El Rahman et al. 2012). In response to A. rabiei, extensive modulation of phytohormone pathways occurred in the host. Our data demonstrate induced expression of JA biosynthesis genes in the resistant genotype; however, SA biosynthesis remained unchanged. Importantly, ET and ABA biosynthesis pathways were also induced, which are not well-known for defense against necrotrophic fungi. This indicates the possibility of their interaction with JA signaling to combat AB stress. The physiologically active hormone levels are controlled through fine-tuning of de novo biosynthesis and catabolism (Saito et al. 2004). Intriguingly, differential regulation of hormone conversion enzymes such as ABA hydroxylase, JMT, SAMT, and methylxanthosine synthase 1 occurred upon A. rabiei challenge. SAMT converts SA into MeSA and methylxanthosine synthase 1 plays a crucial role in alkaloid biosynthesis (Deng et al. 2020; Mizuno et al. 2003). This suggests the active conversion of SA into other forms during infection. Higher expression of JMT, which converts JA into MeJA, mediates jasmonate-regulated plant defense against necrotrophic fungus B. cinerea (Seo et al. 2001). Nowadays, Aux is considered a stress-responsive hormone, besides its role in plant growth and development (Wang and Fu 2011). The downregulation of Aux-related genes suggests suppression of growth response while promoting defense, thus indicating the antagonistic nature of Aux with defense hormones. Together, these hormonal modulations cumulatively orchestrate the early defense response in chickpea and confer AB resistance.

Furthermore, genes involved in plant growth, flowering time, and development were also identified among the DEGs and hub genes. Interestingly, chloroplast and photosynthesis-related genes showed upregulation in the resistant genotype, however, remain unchanged in susceptible. This could possibly be mediated by fungal effectors, as many effector proteins localize to the chloroplast and interfere with photosynthetic processes or ROS levels by interacting with target proteins, thus, promoting susceptibility (Singh et al. 2021b). Moreover, leaf and floral development genes have been located in the QTL regions conferring AB resistance (Daba et al. 2016; Deokar et al. 2019; Kumar et al. 2018). This emphasizes the balanced resource allocation between pathogen combat machinery and growth regulation. Most likely, these genes are employed in the well-known phenomenon of growth-immunity trade-off. However, their functional correlation with chickpea immunity and development needs to be investigated. It will also be interesting to study how pathogen infection affects floral development or vice versa.

In summary, we report the comprehensive transcriptional landscape during Cicer-Ascochyta interaction and divulge the key mechanisms underlying AB resistance and susceptibility in contrasting chickpea genotypes. In addition, we identify early immune signaling components and highlight possible candidate genes that may play important roles in disease resistance. Although the molecular function of identified hub genes requires further experimental investigations, our data strongly indicates their critical role in orchestrating chickpea immunity. Taken together, this study improves our understanding of AB disease development and delivers a handful of candidate genes to be exploited in the genetic manipulation or breeding programs for chickpea improvement.

MATERIALS AND METHODS

Plant material, sample collection, and RNA extraction.

The seeds of FLIP84-92C(2) and PI359075 chickpea genotypes having contrasting phenotypes for AB stress were obtained from the Germplasm Resources Information Network (GRIN). FLIP84-92C is an AB-resistant kabuli chickpea variety with a disease score of 2, whereas PI359075 is a susceptible desi variety with a disease score of 9 (Tekeoglu et al. 2000). Plants were grown under controlled conditions in growth chambers with a 12-h light and dark cycle at 22 ± 2°C and >70% relative humidity. The ArD2 isolate of A. rabiei (Indian Type Culture Collection number 4638) was grown on chickpea-supplemented potato dextrose agar plates. To isolate the fungal spores, A. rabiei–grown plates were flooded with autoclaved distilled water, and the surface was gently rubbed with a sterile loop. The obtained spore suspension was filtered through muslin cloth, and spore concentration was determined using a hemocytometer. For RNA-seq sample collection, the A. rabiei spore suspension containing 10⁶ spores per milliliter was sprayed on 14-day-old chickpea plants. Plants sprayed with sterilized water were used as controls. The aerial tissue of control and infected plants were collected at 24 and 72 hpi, were immediately frozen in liquid nitrogen, and were stored at –80°C. Some of the infected plants from both the genotypes were monitored for two weeks for the progress of AB disease. Total RNA was isolated by utilizing the Plant RNeasy mini kit (Qiagen). The quantity and quality of the RNA samples were checked on 1% denaturing agarose gel, using NanoDrop and Qubit Fluorometer. In total, 24 libraries (CFL, 24FL, 72FL, CPI, 24PI, and 72PI, four replicates each) representing AB-resistant (FLIP84-92C, denoted as FL) and susceptible (PI359075, denoted as PI) genotypes at control (C, uninfected) and AB stress (24 and 72 hpi, A. rabiei–infected) conditions, were sequenced on the Illumina NextSeq500 platform using 2 × 150 bp chemistry.

Data filtering and DEG identification.

The quality of raw reads was analyzed using FastQC. Reads were quality-filtered and trimmed at the ends to remove low-quality bases and adapter sequences, using FASTX-Toolkit. The reads having a minimum PHRED score of 30 were considered for further analysis. The high-quality filtered reads were mapped on the reference genome of chickpea, i.e., CDC Frontier (Varshney et al. 2013), using the HISAT tool with default parameters (Pertea et al. 2016). Further, the mapped reads were assembled via StringTie to generate the reference-guided assembly (Pertea et al. 2016). Read counts were computed using python script prepDE.py3, which created a count file containing the read count matrices for genes. Genes having less than 10 reads on an average per condition were excluded before the differential expression analysis. This filtered gene count matrix was used by a DESeq2 package (Love et al. 2014) in R for the identification of DEGs. DEGs were obtained for all pair-wise combinations (CPIvsCFL, CPIvs24PI, CFLvs24FL, CPIvs72PI, CFLvs72FL, 24PIvs24FL, 72PI vs72FL) with the criteria of adjP value ≤0.05 and fold change ≥1.5.

Functional enrichment and MapMan analysis.

Annotation of the DEGs was performed using Blast2GO application of OmicsBox. In order to identify the over-represented functional categories, GO enrichment analysis was performed using PlantRegMap (Tian et al. 2020). Pathway enrichment analysis was performed by the KOBAS web tool, using hypergeometric model and a significance threshold of P value ≤0.05 (Xie et al. 2011). Venn diagrams, volcano plots, and heat maps were prepared using TB Tool (Chen et al. 2020). Furthermore, MapMan analysis was conducted to visualize the Cicer arietinum gene expression data in the context of biotic stress pathways (Thimm et al. 2004). The Mercator automated annotation pipeline was used to assign Cicer genes to bins. The two-tailed Wilcoxon rank sum test adjusted by the Benjamin-Hochberg method (false discovery rate ≤0.05) was utilized to define the differentially represented MapMan pathways.

Family enrichment, clustering analysis, and TF binding-site prediction.

Among the DEGs, TFs, TRs, and PKs were identified using the iTAK tool (Zheng et al. 2016). Fisher's exact test in R (fisher.test function) was performed to identify significantly over-represented TF families in our datasets. The P value ≤0.05 criterion was considered for significantly enriched families. The normalized count of DETFs, DETRs, and DEPKs were used for hierarchical clustering analysis, using the R hclust function (version 3.6.2). Euclidean distance measure was used to calculate the distance between genes with ward.D as the agglomeration method. The cluster plots were generated with ggplot2 package in R (Wickham 2016). Regulatory sequence analysis tools (RSAT Plants) online server was used to extract the promoter sequences of DEGs. The promoter sequences were subsequently submitted to the binding site prediction tool available on PlantRegMap to identify the TFs having maximum targets in the DEGs. Further, the TF enrichment tool available on the PlantRegMap website was used to find the TFs possessing significantly over-represented targets among the DEGs (P value ≤0.05).

Weighted gene co-expression network analysis.

The co-expression network analysis was performed using the WGCNA package in R (Langfelder and Horvath 2008). Normalized counts of genes in the upper 50% of coefficient of variation for expression across different conditions were used to construct the signed network. The degree distributions in each network followed the power law and satisfied the scale-free topology criterion. The soft threshold power 34 was taken with a minimum module size of 30 genes, mergeCutHeight of 0.15, and deepSplit of 2. The WGCNA modules (co-expression networks) of eigengenes were identified as clusters of highly connected genes. Further, the networks correlated with traits were identified with the criterion of stability correlation P value ≤0.05. The module significance cut-off (Pearson's correlation coefficient) ≥0.5 was considered, and the top-most modules were retained for further analysis. Then, using the GS-MM, hub genes were identified with a cut-off GS and MM ≥0.8 and P value ≤0.05 for a particular condition. Co-expressed genes having weight ≥0.2 were visualized using Cytoscape v3.8.2. The top genes identified through cytoHubba and network analyzer were highlighted in the regulatory network.

qRT-PCR analysis.

To validate the result of RNA-seq, we performed qRT-PCR analysis of 12 randomly selected genes. The gene-specific primers were designed through Primer Express software (Applied Biosystems). The qRT-PCR was carried out on the 7900HT sequence detection system (Applied Biosystems), using SYBR green and 50°C for 2 min, 95°C for 10 min, 40 cycles of 95°C for 15 s and 60°C for 1 min cycling conditions. The data from three biological replicates, each with at least three technical replicates, were used to compute gene expression levels using the 2^−ΔΔCt method. Chickpea β-tubulin (LOC101495306) and EF-1α (LOC105852108) genes were used as internal controls to normalize the data.

Data availability.

All relevant data can be found within the manuscript and its supplementary materials. The raw RNA-seq data are deposited at the National Center for Biotechnology Information Short Read Archive under accessions SAMN29132972, SAMN29132973, SAMN29132974, SAMN29132975, SAMN29132976, and SAMN29132977, which can be retrieved under the BioProject ID PRJNA849761.

AUTHOR-RECOMMENDED INTERNET RESOURCES

FastQC: http://www.bioinformatics.babraham.ac.uk/projects/fastqc

FASTX-Toolkit: http://hannonlab.cshl.edu/fastx_toolkit

ggplot2: https://ggplot2.tidyverse.org

Mercator: http://mapman.gabipd.org/web/guest/mercator

OmicsBox: https://www.biobam.com/omicsbox

PlantRegMap: http://plantregmap.gao-lab.org

prepDE.py3: https://ccb.jhu.edu/software/stringtie/dl/prepDE.py3

RSAT Plants: http://plants.rsat.eu