M. truncatula A17 and F83005.5 plants were grown in hydroponic culture on Farhaeus medium as described by Ben et al. (2013), with 16 h of light (170 μmol m⁻² s⁻¹) at 25°C and 8 h of darkness at 23°C.

V. alfalfae V31-2 (formerly V. albo-atrum), its GFP-expressing transformant Va-A1b2 (Ben et al., 2013), V. albo-atrum LPP0323 (provided by A. v. Tiedemann, Germany), and V. dahliae JR2 (provided by B. Thomma, Netherlands) were grown on PDA medium at 24°C in the dark. Spore suspensions were obtained as described in Ben et al. (2013).

Ten day old plants were inoculated after cutting the root apex by root-dipping in a conidial suspension (10⁶ sp/ml) for 60 min. Then roots were kept in sterile water for 10 min, transferred back to nutritive solution and incubated in a growth chamber at 20°C with 16 h photoperiod. Symptoms were scored on a scale from 0 (healthy) to 4 (dead plant) (Supplementary Figure S1B). Area Under the Disease Progress Curves (AUDPC, Shaner and Finney, 1977) were computed using the “agricolae” package of the R system for statistical computing and graphing (R Core Team, 2012). All data were analyzed with ANOVA using the appropriate model depending on the experimental designs used. When required, data transformations were applied to achieve normality and homoscedasticity of ANOVA residuals. Pairwise treatment differences were determined by a Tukey's test or Newman-Keuls test using the “agricolae” package.

M. truncatula roots were rinsed briefly in water and 1 cm fragments were embedded in 5% low melting point agarose. Longitudinal sections of 45 μm were prepared on a vibratom (Leica VT 1000S) and mounted in distilled water. Confocal images were acquired with a spectral confocal laser scanning system (SP2 SE, Leica) equipped with an upright microscope (DM 6000, Leica, Germany). Phenolic compounds were observed using an inverted microscope (Leica DM IRBE, Leica) under UV light. Images were acquired with a CCD camera (Color Coolview; Photonic Science) using an objective with 40 × magnification.

Total RNA was extracted using TRIzol reagent (Invitrogen) following the manufacturer's instructions. After treatment with RQ1 RNase-Free DNase (Promega) cDNA was synthesized using the ImProm-II™Reverse Transcription System (Promega). Quantitative PCR (qPCR) was carried out with standard protocol (see protocol details in Supplementary Figure S4) using 3 μl of cDNA (3 ng/μl).

Data analysis was performed with the SDS 2.3 software (Applied Biosystems). C_T-values were normalized against the harmonic mean of four reference genes Medtr2g099090.1, Medtr2g033910.1, Medtr3g085850.1, and TC117750 (encoding Glyoxalase 2, Phosphatase 2C, G3PDH, and H3L, respectively). Primers are shown in Supplementary Table S1.

Fold change in gene expression was determined with the ΔΔC_T method, using the mock-inoculated condition as the reference (Livak and Schmittgen, 2001).

Six roots of Va- and mock-inoculated A17 and F83005.5 plants were harvested at t0, 4, 8, and 24 hpi, and shock-frozen in liquid nitrogen. Total RNA was extracted with TRIzol reagent (Invitrogen). Samples from 4, 8, and 24 hpi were equimolarily pooled for mock-inoculated or Va-inoculated plants. Twenty μg RNA were processed for Massive Analysis of cDNA Ends.

MACE was conducted by GenXPro GmbH as described in Behringer et al. (2015). Tag assembly and annotation were conducted by GenXPro GmbH. After sequencing, each 93 nt tag was mapped allowing only 1 mismatch to sequences in databases including M. truncatula genome version 3.5v5 and M. truncatula ESTs of MtGI Release 11.0.

Differential expression analysis was conducted using “edgeR” package of the R/Bioconductor statistical language (Robinson et al., 2010). Reads non-uniquely aligned and in incorrect orientation with respect to genome annotation were discarded. Genes with ≥5 reads were retained.

Normalized read counts were obtained by Trimmed Mean of M-values (TMM) normalization and were modeled by a negative binomial distribution, which allows over-dispersion of counts, and fitted using generalized linear models to test for differential abundances (i) between A17 and F83005.5 “t0” libraries to characterize putative genotype effect on gene expression before inoculation (ii) between mock- and Va-inoculated libraries to detect genes responding to inoculation in each genotype. Biological replicate was fitted as an effect. Differential expression was tested in “tag-wise” mode, where the dispersion is set to the observed dispersion for the transcripts. Logistic regression of the transcript proportions was also used (Collett, 2002). P-values were adjusted to control the false discovery rate using the Benjamini-Hochberg method.

Weighted gene co-expression networks were constructed using the “WGCNA” R package (Langfelder and Horvath, 2008). Transcripts with a coefficient of variation (CV) ≥ 0.35 for normalized read counts across the six libraries from the A17 resistant line were selected. This CV threshold was chosen to include all A17 differentially expressed genes (DEGs) responsive to inoculation in the resulting trimmed dataset. WGCNA was performed using similar settings as in Formey et al. (2014) except a beta power of 20. Eigengene value was calculated for each identified module and used to test for association with different contrasts of biological interest (i.e., “mock-inoculated vs. Va-inoculated roots”) and contrasts accounting for the effect of the biological replicates. The five modules most correlated to inoculation response were visualized using Cytoscape v3.0.1 (Smoot et al., 2011).

Genes from the Greenyellow module possessing Affymetrix probes were analyzed using the legumeGRN gene regulatory network prediction server (http://legumegrn.noble.org/, Wang M. et al., 2013). Unweighted Gene Co-expression Network was predicted using compiled data from the M. truncatula Gene Expression Atlas (MtGea) regarding different pathogenic root interactions. A Pearson's correlation threshold of |0.7| was applied.

A hydroponic culture which allows easy access to roots as described by Ben et al. (2013) was used in this work. After inoculation with V. alfalfae (Va) strain V31-2, the susceptible line F83005.5 developed disease symptoms and aerial fresh biomass was reduced at the end of the experiment, whereas the resistant line A17 was not affected (Supplementary Figure S1). To assess if this difference was correlated to root colonization, both lines were inoculated with Va-A1b2, a GFP-expressing strain, and root sections were observed at different times after inoculation using confocal laser scanning microscopy.

During initial stages of root colonization, no difference between the two genotypes was observed (Figure 1). Conidia were sucked into the xylem vessels at the cut root ends by the transpiration stream, in resistant and susceptible line, as observed 2 h after inoculation (Figures 1A,B). Germination of conidia occurred in both lines, within 24 h (Figures 1C,D). Hyphae developed afterwards and colonized the vessels in both lines similarly (data not shown), as already described for the susceptible line F83005.5 (Ben et al., 2013).

Differences became apparent at intermediate stages [4–7 days post-inoculation (dpi)] (Figure 2). Indeed, at 7 dpi, A17 roots were devoid of Va (Figure 2A) while in roots of the susceptible line, hyphae were abundant in the central cylinder (Figure 2B). Typical for vascular diseases, colonization at this stage was strictly limited to the xylem vessels. At latest stages, the roots of susceptible plants were highly colonized and hyphae were also growing in the cortex (Figure 2D). In contrast, cortical cells of the resistant line showed intensive auto-fluorescence suggesting the accumulation of soluble phenolic compounds (Figures 2A,C; Supplementary Figure S2). Staining with the chitin-binding lectin WGA-FITC did not reveal any hyphae in roots of line A17 at 10 dpi, demonstrating that absence of GFP fluorescence was not due to metabolic inactivity of Va (Supplementary Figure S3).

The disappearance of Va in line A17 as observed by microscopy, was confirmed by quantification of fungal DNA in M. truncatula roots and aerial parts by qPCR using Va-specific primers (Supplementary Figure S4). It can be concluded that the fungus was eliminated from roots of the resistant line A17 plants after initial colonization of the vessels.

Strain V31-2, together with strains V. albo-atrum LPP0323 and V. dahliae JR2 used in previous studies (Negahi et al., 2013), was characterized by PCR with species-specific primers (Inderbitzin et al., 2013; Supplementary Figure S5). The results confirmed that the alfalfa isolate V31-2 belongs to the newly formed species V. alfalfae, and incidentally showed that strain LPP0323 belongs to V. nonalfalfae. PCR with Ave1-specific primers (de Jonge et al., 2012) conducted on these 3 strains showed that neither V31-2 nor LPP0323 contained an Ave1 homolog, in contrast to the positive control Vd JR2. In Ve1 tomato, resistance to race 1 V. dahliae and V. albo-atrum is mediated through recognition of the fungal Ave1 gene. The absence of Ave1 in V31-2 suggests that resistance in M. truncatula toward this strain is different from the Ve1-Ave1 interaction described in tomato.

To get more insight into resistance mechanisms which eliminate Va from roots in line A17, a transcriptomic study was undertaken. Based on root colonization pattern, the period between 0 and 24 hpi was considered as decisive, preceding the accumulation of putative defense compounds, and elimination of the pathogen. Two independent experiments with lines A17 and F83005.5 (mock-inoculated, Va-inoculated) exhibiting similar time course of infection were conducted (Supplementary Figure S6). cDNA libraries prepared from RNA extracted from roots at t0 (before inoculation) and from roots at 4, 8, and 24 hpi pooled thereafter (named “Early” pools) were constructed. Between 13 and 29 million reads were obtained for each of the 12 libraries by MACE sequencing (Supplementary Table S2).

Sequence mapping on the M. truncatula genome gave a total number of 58,186 putative genes. This was reduced to less than half when only genes with ≥5 reads in at least two libraries were considered, with 17,978 sequences for A17 and 18,089 for F83005.5. Among them, 1261 genes were differentially expressed between the two lines at t0, with 735 (58%) over-expressed in A17, establishing a strong genotype effect on basal gene expression (Supplementary Table S3). Interestingly, analysis of GO term enrichment showed that functions associated to “response to stimulus” and to secondary metabolism were strongly enriched for the transcripts with higher expression level in resistant line A17. More strikingly, the category “defense response” was only found in this set of differentially expressed genes (DEGs; Table 1). Line F83005.5 exhibited higher expression of genes involved in primary metabolism, and in categories “biological regulation and signaling,” “localization,” “cellular component organization and process,” and “developmental process.” These strong and large differences in gene expression before inoculation suggest that the resistant line might possess preformed defenses or might be better prepared to induce a fast and efficient response.

Due to the strong genotype effect, independent data analyses were performed for each M. truncatula line. Three groups of DEGs were distinguished, corresponding to genes regulated by time/development (Early-Mock vs. t0), by inoculation (Early-inoculated vs. Early-Mock), or by experimental batch factors (Experiment 1 vs. Experiment 2; Figures 3A,B). In total, 1055 and 4053 DEGs were identified in A17 and F83005.5, respectively, evidencing a stronger transcriptional response in the susceptible line, although colonization at this stage was still very low and not much different from that in the resistant line (Figure 1). Among the 45 DEGs responding to inoculation in the resistant line A17 (0.3% of the 17,978 genes), 13 responded specifically to Va, 29 responded also to development; most genes (72%) were induced (Supplementary Table S4). Among the 302 DEGs responsive to inoculation in the susceptible line F83005.5 1.7% of the 18,089 genes), 114 (37.7%) responded specifically to Va and 111 (36.7%) responded also to development; here, up- and down-regulated genes accounted each for roughly half (45 and 55%) of the DEGs (Supplementary Table S5). The gene expression profiles revealed by MACE were validated by qRT-PCR for 24 DEGs which appeared up- or down-regulated in inoculated A17 and/or F83005.5 (Supplementary Table S6). Gene expression levels were highly similar with both methods (Pearson's correlation coefficient, r = 0.86) supporting the reliability of MACE results.

Twenty DEGs showed similar expression patterns in both lines. Nineteen were functionally annotated, including genes involved in secondary metabolism, oxidation-reduction process, response to stress, and regulation of cellular process (Table 2). Both A17 and F83005.5 showed increased expression of three chalcone synthase genes after inoculation with Va, but induction was stronger in line A17.

To get further insights into the biological pathways involved in the response to Va infection, GO term enrichment analyses were performed for the DEGs responding to fungal inoculation in each line. Regarding line A17, all significantly enriched functional categories contained only up-regulated DEGs and concerned a limited number of biological processes (Table 3). The category “response to stimulus” included notably “response to biotic stimulus,” “defense response,” and “response to organic substance,” and the category “metabolic process” included “aromatic compound,” “amino acid derivative,” and “secondary metabolism.” This picture suggests a targeted and organized biological response of the resistant line early after pathogen inoculation. In contrast, the response to Va in the susceptible line F83005.5 appeared more complex and disordered. Processes related to “response to stimulus” as well as “metabolism process” including very diverse biological pathways were highly affected (Table 4). Among those, single functional categories such as “cellular catabolic process” or “response to biotic stimulus” were both significantly induced and repressed. A majority of cellular pathways such as “response to oxidative stress” and “response to hormone stimuli” appeared to be down-regulated, and all genes of the category “biological regulation,” among them notably “ion homeostasis” were consistently repressed. Strikingly, the GO term associated with “defense response” was not significantly enriched in the susceptible line.

To get a more focused vision into specific cellular processes, we overlaid DEGs responding to Va inoculation onto metabolic networks related to biotic stress response using MapMan (Usadel et al., 2005; Supplementary Figure S7). The response of the resistant line A17 was characterized by induction of a few metabolic pathways putatively involved in pathogenic interactions. In particular “plant secondary metabolism” accounted for eight of the significantly induced genes in A17 (Supplementary Figure S7A), with genes involved in metabolism of lignins and lignans, phenylpropanoids, and flavonoids (including chalcones and isoflavonoids; Supplementary Figure S8A). Genes involved in secondary metabolic process were induced in the susceptible line as well, but with lower expression levels (Supplementary Figure S8B). In addition, many cellular functions were modified (induced and/or repressed) in line F83005.5 in response to Verticillium, including genes participating to cell wall and proteolysis, to defense (mostly PR proteins), to oxidative stress, signaling, and transcription factors. Notably, expression of genes involved in phytohormone signaling was strongly affected in this line: genes related to ethylene (ET), jasmonate (JA), abscisic acid (ABA), and auxin signaling were induced whereas genes of salicylic acid (SA) signaling were repressed (Supplementary Figure S7B).

Experiments where plants were pre-treated with ET, MeJA, ABA, auxin, and SA showed that ABA, SA, and auxin protected against symptoms and fungal colonization and ET delayed disease symptoms whereas MeJA had no effect (Supplementary Figures S9, S10). None of these hormones had a direct effect on the fungus (Supplementary Table S7). Hence, the induction of these genes in line F83005.5 does not seem sufficient to fully trigger defenses regulated by the respective hormone signaling pathways.

Incidentally “Response to abiotic stress” showed a similar contrasted pattern in the two lines: induction of genes in A17 and simultaneous induction and repression of genes in F83005.5.

On the whole, GO enrichment and MapMan analyses suggest a structured molecular response in the resistant line A17 leading to the onset of defense mechanisms early after Va inoculation. In contrast, the susceptible line F83005.5 suffered a complete transcriptional dysregulation of plant metabolism and defense pathways which probably favors disease development.

We further investigated the transcriptome organization in resistant line A17 by applying weighted gene co-expression network analysis (WGCNA) to 2083 sequences exhibiting a coefficient of variation of at least 0.35 across the six libraries (t0, Early-Mock, Early-inoculated, two repetitions). This dataset comprised all 45 inoculation-responsive DEGs identified above. Complementary to the analyses described above which considered only DEGs with a FDR of <0.05, this approach will further reveal subtle effects and show the degree of relationship between genes expression patterns by nodes distance in the graphical presentation. We identified five modules associated to inoculation response which gathered 794 genes (35 to 369 genes per module; Figure 4A, Table 5). This structure points out a coordinated transcriptional response in the incompatible interaction between M. truncatula and V. alfalfae. In contrast, WGCNA on data from line F83005.5 provided only a blurred picture without any distinct co-expression module, suggesting the absence of a structured transcriptional regulation in the susceptible line after inoculation (data not shown).

Statistical analysis of GO term enrichment revealed distinct functional assignments for the different co-expression modules (Table 5). Biological processes involved in response to various biotic and abiotic stimuli were significantly over-represented in all modules, and three of them (Greenyellow, Black, and Darkgray) appeared to be related to specific regulatory pathways such as “negative regulation of biological process” or to different metabolic activities, with some involved in defense against pathogens (e.g., “secondary metabolic process”).

The Greenyellow module is highly enriched in biological functions related to plant defense and encompasses genes exhibiting expression profiles the most correlated to inoculation response. Notably, it contains 29 (64%) of the previously identified DEGs responding to Va inoculation in line A17, together with non-differentially expressed genes (Figure 4B). Hence it appeared to be of major importance for resistance against Va. Detailed examination of this resistance-associated module led to the identification of genes known to play key roles in plant-pathogen interactions, such as genes involved in plant defense (e.g., encoding chitinases, disease resistance RPM1-like proteins, pathogenesis-related protein PR10 or PR5-like receptor kinase (PRK5), LRR protein, …), in ROS metabolism (e.g., peroxidases…), and in secondary metabolism (e.g., PAL, chalcone synthases,…; Supplementary Table S8). Phytohormone-responsive genes, most of them related to ET, ABA, and auxin signaling, are also present in this module, together with genes encoding receptor-like kinases and protein kinases and 21 putative transcription factors (TFs) from different gene families (namely, AP2/ERF and B3, ethylene responsive, MYB and MYB/SANT, NAC, bHLH, WRKY, TIFY, LOB, Trihelix family protein, and C2H2-type zinc finger protein families). Homologs of several genes involved in PAMP-triggered immunity (PTI) were identified, such as BAK1 (FM886833.1), MKK2 (Medtr1g113960.1), EF-Tu Receptor (EFR, Medtr5g026090.1), NDR1 (contig_99115_1.1), and Wat1 (walls are thin 1; Medtr3g072500). A KEGG analysis revealed the presence of additional homologs of genes related to PTI, namely CNGC10 (cyclic nucleotide gated channel 10), HSP90 (heat shock protein 90), and CERK1 (chitin elicitor receptor kinase 1). These genes may play a crucial role in the perception of Va and in signal transduction in M. truncatula roots as well as in the subsequent transcriptional regulation leading to resistance.

Finally, to address the question if the Greenyellow co-regulatory module was specific for M. truncatula resistance toward Verticillium or potentially involved in the response to other root pathogens, we used the LegumeGRN prediction server (Wang M. et al., 2013). An in silico analysis of co-expression patterns of the 333 (90%) genes from the Greenyellow module which possess Affymetrix probes on transcriptional data of M. truncatula roots inoculated with Aphanomyces euteiches, Macrophomina phaseolina, Phymatotrichopsis omnivore, and Ralstonia solanacearum led us to identify a stringent co-expression network gathering 186 genes (Figure 5). This highlights a substantial level of conservation in the transcriptional response of M. truncatula toward various root pathogens. Indeed, a high proportion (56%) of the genes with Affymetrix probes from the Greenyellow module are also co-regulated in response to other root pathogens (either fungi, bacteria, or oomycetes). A core of 15 TFs holds central positions within this co-regulatory module. Hence, not only defense mechanisms but also important regulators of gene expression are conserved among the responses of M. truncatula to different root pathogens.

The response against V. alfalfae in resistant M. truncatula is different from the well-characterized Ve1 dependent resistance against V. dahliae and V. albo-atrum in host plants such as tomato, cotton, or Solanum torvum (Fradin et al., 2009; Liu et al., 2012; Zhang et al., 2012). Homologs of Ve1 were detected in the M. truncatula genome but were not localized under the resistance QTLs against Va V31-2 (Ben et al., 2013) and are not expressed at important early stages of the interaction, since they were not detected in our MACE libraries. Ve1 recognizes Ave1, the Verticillium effector of race1 strains (Fradin et al., 2011). Ave1 homologs were found in three fungal and many plant species; the Fusarium and Cercospora homologs were both recognized by tomato Ve1, supporting the hypothesis that Ve1—Ave1 mediated resistance had traits of PAMP-triggered immunity (de Jonge et al., 2012). Va strain V31-2 does not possess Ave1 homologs although they were reported in the alfalfa strain VaMS102 (de Jonge et al., 2012).

The study of root colonization showed that after initial colonization of the stele, the fungus was eliminated from roots of the resistant line A17. In other plants, resistance takes place in the shoot. V. dahliae and V. longisporum were shown to infect roots of resistant and susceptible hosts but resistant hosts inhibited pathogen development in the shoot (Heinz et al., 1998; Eynck et al., 2009; Zhang W. W. et al., 2013). The destruction of the fungus by the plant, and the absence of both Ave1 in Va strain V31-2 and of Ve1 in line A17 suggest that resistance to Verticillium wilt in M. truncatula A17 is novel and different from that described in other plants.

A common transcriptional response involving a few genes related to flavonoid biosynthesis and stress response was found in both M. truncatula lines. However, cytological studies showed that only line A17 accumulated autofluorescent soluble compounds, probably flavonoids, in inoculated roots. Such autofluorescent phenolic compounds were also induced by A. euteiches in roots of the resistant line A17 and to a lesser degree in the susceptible line F83005.5 (Djébali et al., 2009). Induction of genes involved in phenylpropanoid, terpenoid, and flavonoid biosynthesis is the most common response to Verticillium infection in resistant and in susceptible plants (Pegg and Brady, 2002; Fradin and Thomma, 2006). Accumulation of medicarpin, the major isoflavonoid phytoalexin in Medicago spp. (Naoumkina et al., 2010), has been observed in alfalfa callus inoculated with V. albo-atrum (Latunde-Dada and Lucas, 1986). The antifungal properties of phenolics known to be involved in defense against Verticillium wilt (Beckman, 2000; Báidez et al., 2007; Eynck et al., 2009). The phenolics-storing parenchyma cells may release these compounds into the xylem vessels to inhibit pathogen growth.

The resistant line A17 exhibited a higher level of gene expression for the functional category “secondary metabolism” at t0, compared to line F83005.5. This might indicate a state of “preparedness” and allow for a more efficient induction of phenolics after inoculation with Va (this study) or A. euteiches (Djébali et al., 2009). It has been reported that Va and other alfalfa pathogenic fungi are able to metabolize medicarpin (Soby et al., 1996) which might explain that phenolic compounds were not observed in the heavily colonized susceptible line.

Analysis of the response to Va inoculation by GO term enrichment showed that resistant line A17, in addition to genes of secondary metabolism, induced genes related to stress response, and notably to defense such as genes for PR4 protein and BetV1 proteins belonging to the PR10 family (Breiteneder et al., 1989).

In contrast, the susceptible line F83005.5 was severely affected in response to stress, defense mechanisms and in several metabolic pathways. Such a strong effect at early infection stages when only a very small part of the root is colonized, suggests that V31-2 is able to repress defense mechanisms and dysregulate plant metabolism by secreted effectors or toxins. Verticillium is described to produce several toxins which induce wilting symptoms but may also act as elicitors of defense, depending on their concentration and/or the plant genotype (Meyer et al., 1994; Wang et al., 2004; Palmer et al., 2005). Preliminary experiments showed that culture filtrate of strain Va V31-2 induces wilting symptoms and necrosis in M. truncatula (data not shown). SA has been reported to protect cotton callus against Verticillium toxin (Zhen and Li, 2004). The SA pathway was down-regulated in the susceptible M. truncatula line, notably genes encoding the SA receptor SABP2 (Kumar and Klessig, 2003) and the SA-regulated glutathione S-transferase (Uquillas et al., 2004), which might lead to enhanced susceptibility to toxins putatively produced by Va during infection.

By applying WGCNA using the MACE data of resistant line A17 we identified one module (called Greenyellow) of the co-expression network which was highly related to inoculation response and contained genes encoding LRR-RLKs, NB-LRRs or regulatory elements including transcription factors.

Various gene homologs encoding proteins involved in PAMP perception were detected, including BAK1, EFR, the chitin receptor CERK1, and PRK5 (Miya et al., 2007). Chitinases which release chitin fragments from fungal cell walls (Adams, 2004) were also present in the module. PRK5 encodes a protein kinase receptor with a PR5 family extracellular ligand binding domain (Wang et al., 1996). PR5 specifically binds and hydrolyzes β-1,3-glucans, components of fungi, and oomycete cell walls (Osmond et al., 2001). PRK5 probably fixes the same ligands as PR5. The receptor kinase BAK1 is required for the function of certain Pattern Recognition Receptors (PRR) lacking cytoplasmic signaling domain such as Ve1 (Heese et al., 2007; Wang et al., 2010; Fradin et al., 2014). BAK1 has been shown to be involved in Arabidopsis resistance to V. longisporum (Roos et al., 2014), together with the brassinolide, JA and ABA signaling pathways. BAK1 interacts with the receptor EFR (Heese et al., 2007), which recognizes the bacterial PAMP EF-Tu (Zipfel et al., 2006). Surprisingly, BAK1 was also shown to be involved in the susceptibility of Arabidopsis Col-0 to V. dahliae (Gkizi et al., 2016).

Several homologs of genes involved in PTI signaling, such as NDR1, EDR1, and MKK2 (or MEK2) were also present in the resistance-associated module Greenyellow. NDR1 which is required for the activity of several CC-NB-LRR R proteins such as RPM1 (Aarts et al., 1998), participates in ETI and in PTI (Knepper et al., 2011), while EDR1 is a MAP kinase which exerts negative feedback on the activation of MAPK cascades and plays a key role in regulating defense mechanisms (Zhao et al., 2014).

The whole picture suggests that the M. truncatula response may resemble innate immunity and PTI rather than the more specific ETI, similar to the situation in resistant tomato and cotton. Indeed, NDR1, BAK1, and MEK2 were found to be required in tomato (Fradin et al., 2009), and FLS2, EFR, CERK1, BAK1, RPM1, PRS2, and MKK2 were found to be induced in resistant cotton G. raimondii after inoculation with V. dahliae (Chen et al., 2015). However, in M. truncatula the response is triggered without implication of a Ve1 homolog, suggesting a different mechanism of recognition of Verticillium and possibly of signaling.

In addition to genes involved in PAMP perception and signal transduction, the Greenyellow module contains genes related to hormone signaling. Although MapMan and GO term enrichment did not show significant induction of these genes in response to inoculation or over-expression at t0 in the resistant line, their co-regulation suggests their involvement in the resistance response. We identified the nodulin MtN21 (a homolog of Wat1), and mlo2 (mildew resistance locus 2 O) in this co-expression module. MLO proteins are involved in resistance against downy mildew in barley but also in dicot plants and prevent pathogen penetration (Consonni et al., 2006; Acevedo-Garcia et al., 2014). MtN21/WAT1 involved in secondary cell wall formation The Arabidopsis Wat1 mutant is more resistant to several vascular pathogens, probably due to constitutive activation of the SA pathway and inhibition of the auxin pathway (Denancé et al., 2013). Such effects on SA and auxin pathways were also observed in the mlo2 mutant (Consonni et al., 2006). The expression of these homologs in M. truncatula indirectly suggests activation of the auxin pathway in response to V. alfalfae.

WGCNA showed that the resistance-related response of line A17 was associated with induction of 21 putative TFs, most of them (71%) being also co-expressed in response to other pathogens of M. truncatula. Four among them are involved in ET signaling (Medtr1g069960.1, Medtr1g087920.1, Medtr1g093600.1, Medtr6g037610.1). In addition, genes encoding for ACC oxidases and genes with GO annotations “response to ethylene stimulus” like EBF1 (EIN3-binding F-box protein 1; Medtr4g114640.1) were identified. This indicates that ET signaling is activated in response to V. alfalfae. However, since ET signaling genes were also induced in the susceptible line and ACC pretreatment did only delay disease symptoms, ET does probably not contribute significantly to resistance.

Some TFs of the Greenyellow module belong to the MYB, NAC, bHLH, WRKY, and ERF TFs families that are also the most represented in the cotton transcriptome in response to V. dahliae (Zhang Y. et al., 2013). WRKY TFs play a crucial role in defense mechanisms (Pandey and Somssich, 2009), NAC TFs are involved in the response to abiotic stress and pathogens (Nuruzzaman et al., 2013) and bHLH93 is implicated in flowering in A. thaliana. Interestingly, flowering time and symptom severity were linked in V. longisporum—inoculated Arabidopsis (Häffner et al., 2010). These particular TFs were also found to participate in some ABA responses (Cutler et al., 2010; Wang Y. et al., 2013), and among the five genes encoding for PR10 proteins in this module, three encode for ABA-responsive protein ABR17 with a BetV1 domain. These results together with efficient protection by ABA pretreatment give some indication that the ABA pathway might be involved in M. truncatula resistance against V. alfalfae.

A large number of genes identified in the resistance-associated Greenyellow module was also found to be co-expressed during the interaction of M. truncatula line A17 with four other root pathogens, among them 15 TFs (Figure 5). Because line A17 is susceptible toward three of these other pathogens, co-expression of these genes does not suggest, at first glance, a straightforward and decisive role uniquely in disease resistance. However, the experimental conditions were differing between the pathosystems and the data analyzed usually represent only a snapshot of the whole infection process. As a consequence, differences in intensities and/or timing of gene expression in this composite dataset might still not be sufficient to explain resistance to one given pathogen and susceptibility to others. We advocate that there is probably not one ultimate master switch responsible for the phenotype in QDR, but rather a collaboration of several regulatory genes which contribute to the final result. The core pathogen-responsive genes identified in our study might thus still participate in basal defense of M. truncatula against root pathogens.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.