Plasticity of Phymatotrichopsis omnivora infection strategies is dependent on host and nonhost plant responses

Abstract Necrotrophic fungi constitute the largest group of plant fungal pathogens that cause heavy crop losses worldwide. Phymatotrichopsis omnivora is a broad host, soil‐borne necrotrophic fungal pathogen that infects over 2,000 dicotyledonous plants. The molecular basis of such broad host range is unknown. We conducted cell biology and transcriptomic studies in Medicago truncatula (susceptible), Brachypodium distachyon (resistant/nonhost), and Arabidopsis thaliana (partially resistant) to understand P. omnivora virulence mechanisms. We performed defence gene analysis, gene enrichments, and correlational network studies during key infection stages. We identified that P. omnivora infects the susceptible plant as a traditional necrotroph. However, it infects the partially resistant plant as a hemi‐biotroph triggering salicylic acid‐mediated defence pathways in the plant. Further, the infection strategy in partially resistant plants is determined by the host responses during early infection stages. Mutant analyses in A. thaliana established the role of small peptides PEP1 and PEP2 in defence against P. omnivora. The resistant/nonhost B. distachyon triggered stress responses involving sugars and aromatic acids. Bdwat1 mutant analysis identified the role of cell walls in defence. This is the first report that describes the plasticity in infection strategies of P. omnivora providing insights into broad host range.


| INTRODUCTION
Plant pathogens are categorized based on their mode of acquiring nutrition (Laluk & Mengiste, 2010). Biotrophic pathogens derive nutrition from live plant cells, whereas necrotrophic pathogens derive nutrition from dead or dying plant cells. The hemi-biotrophs have a transient biotrophic phase before becoming necrotrophs. Necrotrophic fungal pathogens are either host specific or broad spectrum. Broad-spectrum necrotrophic fungal pathogens have myriad strategies to infect plant species including production of enzymes, toxins, and effectors, thus causing heavy crop losses annually (Laluk & Mengiste, 2010).
Phymatotrichopsis omnivora (Duggar) Hennebert (G. M. Watkins & Watkins, 1940) is one such broad-spectrum, filamentous, soil-borne, necrotrophic pathogen that causes the destructive Phymatotrichopsis root rot (PRR) disease in Southwest USA and Northern Mexico (Uppalapati et al., 2010). This facultative saprophytic fungus becomes pathogenic during the dry summer months (Rogers, 1942). Penetration may occur through wounds or by mechanical action on the periderm of the roots in the field (Peltier, King, & Sampson, 1926). The typical disease symptoms include chlorosis, rapid wilting, and plant death. P. omnivora infects over 2,000 dicotyledonous plants but cannot infect monocotyledonous plants (Streets, 1937). Several members of Brassicaceae have been reported to escape PRR disease when grown as winter crops (Streets, 1937). P. omnivora causes severe disease in fibre and forage crops like cotton and alfalfa (Medicago sativa), respectively (Lyda, 1978).
Farmers in Southern Oklahoma and Texas are reluctant to grow alfalfa in spite of its high economic and nutritional value due to persistence of PRR disease in this region. Since the first report of this disease in late 1800s (Pammel, 1888), no resistant cultivars have been identified. The molecular mechanisms of pathogen virulence, host susceptibility, and broad host range are not yet understood (Uppalapati et al., 2009).
Extensive attempts to identify P. omnivora-resistant M. truncatula or M. sativa genotypes were unsuccessful. Thus, to further understand the P. omnivora-host or P. omnivora-nonhost interactions and the broad host range of P. omnivora, A. thaliana and B. distachyon were included in this study. We performed comparative cell biology and transcriptional profiling of susceptible (M. truncatula), partially resistant (Arabidopsis thaliana), and resistant (nonhost; Brachypodium distachyon) interactions during P. omnivora infection. We found that P. omnivora exhibits high fungal plasticity and can infect plants either as a necrotrophic pathogen or as a hemi-biotrophic pathogen. This plasticity provides insights into its broad host range. We also report two distinct plant stress responses when challenged with P. omnivora.

| Plant growth, fungal growth, and infection assays
The plants were grown in culture tubes as described previously (Uppalapati et al., 2009) with the following changes. M. truncatula (A17) seeds were scarified with sandpaper. B. distachyon (Bd21-3) seeds were dehusked and surface sterilized by vortexing in 30% bleach solution for 7 min. A. thaliana (col 0) seeds were sterilized in 95% ethanol for 2 min, 30% bleach for 5 min, and followed by four washes in distilled water.
A. thaliana seeds were stratified for 2 days at 4 C. P. omnivora isolate NFPO01 was isolated from infected alfalfa roots at Ardmore, OK in summer 2014 and cultured as described previously (Uppalapati et al., 2009). Four-week-old seedlings were infected with P. omnivora with wheat seed inoculum as described previously (Uppalapati et al., 2009).

| Light, confocal, and scanning electron microscopy
Infected roots were stained with 10-μg/ml wheat germ agglutinin coupled to green fluorescent dye Alexa Fluor 488 (WGA Alexa Fluor 488; Invitrogen Corp., Carlsbad, CA, USA), which was dissolved in distilled water for 30 min. The roots' cell walls were counterstained with 10-μg/ml propidium iodide dissolved in distilled water for 15 min. Imaging was done as described previously (Ray, Guo, Kolape, & Craven, 2017

| Transcriptional profiling and RT-qPCR
Infected roots at 0, 1, 3, 5, 7, and 10 days post infection (dpi) were frozen in liquid nitrogen followed by homogenization with a mortar and pestle. RNA extraction was done with the Spectrum™ plant total RNA kit (Sigma-Aldrich, cat # STRN50 pooled in equimolar ratio. The pooled libraries were sequenced on a NextSeq 500 Sequencing system (Illumina, CA, USA). RT-qPCRs were done as described previously (Gill et al., 2018) with the primer sets listed in supporting information Table S1.
Count data was then obtained using Kallisto (kmer length = 31; Bray, Pimentel, Melsted, & Pachter, 2016). To improve the confidence in the read quantification, we performed 1,000 bootstrapped samples and took the average as our transcript quantification. Expression count for each respective plant species were analysed after trimmed mean of M values normalization (Robinson, McCarthy, & Smyth, 2010). This generated a matrix with samples as the columns and transcripts as the rows. Alignment and read statistics are provided in Table S2. The estimated count data was then processed using VOOM (Law, Chen, Shi, & Smyth, 2014) for the purposes of differential analysis using the R package LIMMA (Smyth, 2005). Multiple hypotheses correction was performed using the Benjamini-Hochberg procedure, with a false discovery rate cut-off of 0.01 (Benjamini & Hochberg, 1995). For differential comparisons, we compared the transcript counts of the time points post infection with the time 0 control. In addition, differential comparisons were also performed between consecutive time points post infection.
Gene ontology enrichments were done using AgriGO (Tian et al., 2017) for the three plant interactions. Pathway enrichments and protein domain enrichments in M. truncatula were performed using MedicMine (Krishnakumar et al., 2015) and in A. thaliana were performed using ThaleMine . The co-expression analysis was performed on the genes that are statistically significantly different compared with genes in control sample (0 dpi). The expression values were grouped into early stage (1 and 3 dpi) and late stage (7 and 10 dpi), respectively. For each of the respective early-stage/ late-stage data sets, an all-pairs Pearson correlation coefficient calculation of the differential transcripts was performed. Gene pairs with positive correlation above 0.80 or a negative correlation below −0.75 were further analysed as putatively co-expressed genes. Gene enrichment analysis was performed on the early-stage/late-stage categorized normalized gene counts using the right-tailed Fisher exact test, accounting for multiple hypotheses bias (false discovery rate < 0.05). Further enrichment analyses were performed based on culture tube and cell biology observations overlapping specific time points post infection, as given in Table 1. Mercator (Lohse et al., 2014) was used to annotate the reference transcripts and obtain respective Mapman descriptive terms. The co-expression clusters, differential expression results, along with the Mapman annotations were all visualized in Cytoscape 3.6 (Shannon et al., 2003). A further combined analysis of all the expression count matrices of the plant species was performed, after collectively normalizing the raw expression count data from the plant species using trimmed mean of M values normalization (Robinson et al., 2010). An all-pairs Pearson correlation coefficient calculation was then performed. Gene pairs that had an absolute correlation value greater than .98 were then clustered using Markov clustering (Enright, Van Dongen, & Ouzounis, 2002), with an inflation value of 7.

| RESULTS
3.1 | M. truncatula, B. distachyon, and A. thaliana respond differently to P. omnivora infection in culture tube assays We used a previously established culture tube-based disease assay system (Uppalapati et al., 2009) with three model plant species, M. truncatula (susceptible), B. distachyon (nonhost/resistant), and A. thaliana (partially resistant). In the susceptible M. truncatula, the fungal mycelia grew into the culture media tubes from 1 dpi onwards leading to wilting by 10 dpi and chlorosis and death by 18 dpi ( Figure S1). The fungal growth features were similar in the resistant B. distachyon culture tubes except that the plants did not wilt and die. Although B. distachyon plants were stressed by 18 dpi with yellow leaves, they survived for 8 weeks ( Figure S1). Interestingly, P. omnivora growth was inhibited in A. thaliana at 1 dpi. The mycelia grew slowly by 3 dpi into the agar, avoiding roots when possible. The fungal growth was active at 5 dpi, still avoiding the roots when possible. By 10 dpi, the A. thaliana plants exhibited leaf yellowing and died by 18 dpi ( Figure S1). Dark pigmentation was seen in the agar at the root-shoot junctions in all the three plant species between 5 and 7 dpi when the mycelia grew to the bottom of the tubes, potentially increasing the nutrient competition between plant and fungus.
3.2 | P. omnivora grows as a necrotrophic pathogen in susceptible M. truncatula  (Table S3). These results indicate that A. thaliana efficiently blocks P. omnivora growth initially, but the fungus eventually penetrates the cells and grows intracellularly. Taken together, these data indicate that P. omnivora employs two different infection strategies in the susceptible and partially resistant plants.
Further, the difference in infection biology in the resistant and partially resistant plants ( Figure S2) also suggests differential host perception and defence responses.

| Transcriptional analysis indicated the infection process as a gradual step-up phenomenon
To determine plant molecular responses of all three species tested to P. omnivora infection, RNA sequencing was conducted in the  (Table S4). In the resistant/nonhost B. distachyon, cellular carbohydrate metabolism processes and genes involved in cellular detoxification like peroxidases and monoxygenases were upregulated from 1 dpi onwards and sustained throughout the infection process. At 3 dpi, protein kinase SD-2b family genes, LRK10L-2 subfamily genes, and genes involved in lignan synthesis, which play a role in plant defence, were upregulated. At 7 and 10 dpi, genes encoding aldolase, chorismates, steroids, and other aromatic amino acids were significantly upregulated (Table S6a-b).
GO enrichments in partially resistant A. thaliana indicated that the plant launches an active defence response from 1 dpi onwards (Table S7a- Table S8).

| Gene enrichments at early and late infection time points indicated distinct molecular processes were induced in the three plant species
We did gene enrichments corresponding to early (1 and 3 dpi) and late (7 and 10 dpi) infection time points to identify key plant processes.  (Tables 1 and S9).
F I G U R E 6 Differentially expressed gene (DEG) analysis across three plant species and defence-related gene analysis. (a) Comparative gene ontology enrichment of upregulated and downregulated DEGs across three plant species. Scale bar indicates increasing significance levels from 1 to 9. The details of these gene ontology classes are listed in Table S8a-b. (b) Defence-related gene analysis in the upregulated and downregulated DEGs in three plant species. Further details of these genes are listed in Table S11 In the partially resistant A. thaliana, unique disease-responsive genes were induced at 1, 3, 7, and 10 dpi (Tables 1, and S9) 3.8 | Nonhost and partially resistant plants induce two distinct stress response pathways against P. omnivora We examined the expression of genes that encode proteins involved in plant immunity summarized in Table S10 in (Table S11). The B. distachyon F I G U R E 7 Correlation network data analysis. Graphs (a,b) indicate the number of hub genes in functional gene categories in twofold and above upregulated differentially expressed genes in the three plant species at early (1 and 3 days post infection [dpi] compared with control) and late infection stages (7 and 10 dpi compared with control), respectively. (c) Interspecies correlation networks clusters of differentially upregulated genes at 5 dpi. Orange edges indicate Arabidopsis thaliana genes, purple edges indicate Brachypodium distachyon genes, and green edges indicate Medicago truncatula genes. The three plant species formed separate clusters indicating absence of interspecies correlations. Datasets S12-S14 were used for correlation network analysis interaction was characterized with the consistent upregulation of a BAK1 homolog (Bradi4g18280) from 3 to 10 dpi (Figure 6b).  Figure S4). The signalling gene At5g64890 was induced from 3 dpi onwards and encodes for the PROPEP2 small peptide. PROPEP2 is processed into PEP2 and binds to PepR2 DAMP receptor that triggers defence responses ( Figure S4). This analysis further confirmed that the two plants induced different defences.
3.9 | Correlation networks and hub genes determined diverse functional pathways in A. thaliana and B. distachyon To characterize the functional pathways during plant responses to P. omnivora, we generated correlation networks for twofold and above (c) Wild-type BD21-3 roots infected with P. omnivora stained with wheat germ agglutinin Alexa 488 (green) and propidium iodide (magenta). The fungal hyphae grow intercellularly, indicated by the green colour. (d) Collapsed image of 17 z-stack optical sections imaged at 1.5-μm thickness of Bdwat1 mutant root infected with P. omnivora. Stars indicate penetration point into epidermal cells. Arrow indicates intracellular hyphae growing from cell to cell. The arrow heads point to the hyphae growing on the surface of the root that is stained with wheat germ agglutinin Alexa 488 in green coloured hub genes that are hypothesized to play significant roles in the infection process (Albert, Jeong, & Barabási, 2000;Langfelder, Mischel, & Horvath, 2013). Mapman annotation was used to identify functional classes of these hub genes. M. truncatula had both abiotic and biotic stress response genes, secondary metabolism genes, signalling receptor kinases, pectin methyl esterases (PMEs), and cell wall degradation genes upregulated at the early infection stage (1-3 dpi). By the late infection stage (7-10 dpi), the genes encoding signalling receptor kinases, secondary metabolism, and PMEs were downregulated (Figure 7a-b). The resistant B. distachyon plants had enhanced abiotic stress response genes.

| DISCUSSION
We conducted a comprehensive cell biology and transcriptional analysis of P. omnivora interactions in three different plant species: M. truncatula, B. distachyon, and A. thaliana. The stronger virulence of the pathogen in M. truncatula was supported by the fact that P. omnivora caused typical disease symptoms like mycelial mantle formation, wilting, and chlorosis as described previously (Uppalapati et al., 2009;G. M. Watkins, 1938;G. M. a. W. Watkins, M.O., 1940). In previous studies, M. truncatula defence responses were induced by 3 dpi but returned to basal levels by 5 dpi (Uppalapati et al., 2009) (Vidhyasekaran, 1974). In a more recent study, cold acclimatized grasses accumulated sugars that enhanced their resistance to fungal pathogens (Rapacz, Pla_ zek, & Niemczyk, 2000). The study of vascular wilt pathogen Fusarium oxysporum on Lupinus luteus grown with and without sucrose indicated that sugars may be involved in resistance mechanism (Morkunas, Marczak, Stachowiak, & Stobiecki, 2005 were also upregulated in the P. omnivora-A. thaliana interaction as part of ETI response. The induction of these NBS-LRR genes coincided with the induction of SA-dependent defence responsive genes and kinase-encoding genes. SA defence pathway gene NDR1 and several mitogen-activated protein kinase pathway genes that were induced in P. omnivora-A. thaliana interaction were implicated in resistance mechanisms to other fungal pathogens (Century et al., 1995;Genenncher et al., 2016;Song et al., 1995;Wang, Song, Ruan, Sideris, & Ronald, 1996).
One of the initial responses of A. thaliana to P. omnivora infection was the induction of genes encoding sHSPs and regulators of sHSPs.
These included AtHSP90/70/83/81/40. The roles of AtHSP90/70/40 in plant immunity were previously described (Hubert et al., 2003;Park & Seo, 2015). The defence mechanism in A. thaliana was also characterized by the induction of WRKY genes and glucosinolate production genes, both of which were uniquely upregulated in this inter- susceptibility to the pathogen indicating that PEPR-mediated defences play an important role in the basal immunity response against P. omnivora as early as 1 dpi. The PEPR pathway induced basal resistance when microbe-associated molecular pattern-triggered defences were compromised during Colletotrichum higginsianum infection in A. thaliana (Yamada et al., 2016). In our study, SA-mediated defence pathways were upregulated along with genes encoding JAZ proteins, indicating suppression of JA-mediated biotic defence pathway from 3 dpi onwards. This corresponded to slow growth of the fungus into the culture tubes as well as intracellular root epidermal colonization. The invasive hyphae did not stain with WGA Alexa 488 unlike the intercellular hyphae. Alexa 488 binds to chitin in cell walls. Plants recognize chitin and induce primary immune responses.
In order to evade these responses, pathogenic fungi are known to alter the localization of chitin in their cell walls. Studies in Magnaporthe grisea showed that the bulbous invasive hyphae localize chitin further into the cell wall, whereas α-1,3, glucans and chitosan were localized into the accessible regions of the cell walls as a strategy to circumvent host recognition (Fujikawa et al., 2009). We hypothesize that P. omnivora alters chitin localization in a similar manner to evade detection during root epidermal cell penetration in A. thaliana.
In the susceptible M. truncatula, P. omnivora adopts a necrotrophic infection strategy, effectively damaging cells ahead of its infection indicated by apoplastic and ER markers. However, a similar attempt in A. thaliana launched an efficient defence response against the P. omnivora. In order to evade these defences, P. omnivora appeared to alter its infection strategy, characterized by intracellular invasive hyphae that had potentially altered chitin localization in the outer cell This research indicates exploring novel avenues to accomplish PRR disease control methods in the field. These could involve RNA interference silencing-based approaches or small peptide-based approaches.