RED HOT Contributors

 

Divergent and convergent modes of interaction between wheat and Puccinia graminis f. sp. tritici isolates revealed by the comparative gene co-expression network and genome analyses

0

Comparative genomic analysis of Pgt isolates MCCFC and RKQQC

Genomic sequence similarity was assessed between the two North American isolates of

Pgt

59KS19 (henceforth, referred to by its race designation MCCFC), 99KS76A-1 (henceforth, referred to by its race designation RKQQC) [

40

], and the previously sequenced and annotated isolate 75-36-700-3 (henceforth, referred to by its race designation SCCL) [

41

]. Genomic DNA was extracted from the MCCFC and RKQQC urediniospores and sequenced using both short (2 × 100 bp Illumina reads) and long read (Pacific Biosciences) technologies (Additional file

1

: Table S1). In total, 68.3 and 62.8 million quality filtered Illumina paired-end (PE) reads were generated for MCCFC and RKQQC, respectively, and were aligned to the SCCL genomic reference. A total of 50.1 million genomic reads (73.3%) generated from RKQQC, and 25.8 million genomic reads (41.0%) generated from MCCFC were mapped to the SCCL genome. The RKQQC reads covered 84.5% of the SCCL genome at a mean coverage of 81.72×, and the MCCFC genomic reads covered 75.8% of the reference genome at a mean coverage depth of 49.57× (Table 

1

).

Table 1

Summary of the genomic reads from RKQQC and MCCFC aligned to the SCCL reference genome

Total reads

62,781,618

68,290,518

Aligned reads

46,024,317

28,014,613

Proportion reads mapped

73.31%

41.02%

Mean coverage depth

81.72

49.56

Median coverage depth

78

49

% of genome with no coverage

15.52%

24.24%

% of genome with 5× coverage

77.96%

73.80%

% of genome with ≥10× coverage

75.07%

70.82%

To characterize the genomic diversity of the newly sequenced

Pgt

isolates, the paired-end genomic reads mapped to the SCCL genome were processed using the variant calling algorithms implemented in the GATK UnifiedGenotyper [

27

]. A total of 782,717 (read coverage depth ≥ 5 per allele) polymorphic sites including 89,417 indels as well as 693,300 SNPs were identified in the MCCFC and RKQQC genomes. Amongst all the SNPs identified by comparing two isolates, 347,376 (50%) discriminated between MCCFC and RKQQC genotypes, suggesting the high level of inter-isolate genetic divergence (Table 

2

). The potential functional effects of the identified divergent SNPs were assessed using the publically available SCCL gene models (Table 

3

). In total, 6,624 (41.92%) of the gene models from the SCCL genome showed the presence of non-synonymous SNPs.

Table 2

Summary of variant calls generated using the UnifiedGenotyper from the GATK package

Total number of biallelic sites called

782,717

Small insertions/deletions

89,417

SNPs

693,300

Informative SNPs with adequate read coverage in both isolates

557,819

Informative sites with the same genotype call in both isolates

210,443

Discriminatory sites with different genotype call in each isolate

347,376

Discriminatory sites heterozygous in MCCFC

71,264

Discriminatory sites heterozygous in RKQQC

104,945

Table 3

Assessment of the tentative functional impact of SNPs using SNPeffect program

All SNPs in both isolates

 Non-synonymous mutations

3.11

1

106

8620

0.5455

 Synonymous mutations

4.48

1

88

9065

0.5737

 No mutations

NA

NA

NA

5526

0.3497

SNPs in RKQQC alone

 Non-synonymous mutations

2.81

1

106

8251

0.5222

 Synonymous mutations

4.09

1

84

8738

0.553

 No mutations

NA

NA

NA

5850

0.3703

SNPs in MCCFC alone

 Non-synonymous mutations

2.63

0

106

7751

0.4906

 Synonymous mutations

3.91

1

88

8196

0.5187

 No mutations

NA

NA

NA

6432

0.4071

SNPs differentiating RKQQC from MCCFC

 Non-synonymous mutations

1.83

0

96

6624

0.4192

 Synonymous mutations

2.6

0

71

7323

0.4635

 No mutations

NA

NA

NA

7247

0.4587

Previous diversity studies of the effector encoding genes showed the elevated proportion of non-synonymous SNPs suggesting that this class of genes is likely subject to directional selection [11]. Consistent with these observations, the proportion of non-synonymous and synonymous SNPs in the effector encoding genes (3,567/4,751) compared to the remaining genes (23,405/37,191) in the Pgt genome (Additional file 2: Table S2) in our dataset was significantly different (χ 2 test = 55.5, p-value = 9.4 × 10−14). The mean ratio of non-synonymous to synonymous SNPs in the effector encoding genes (0.98) was 12% higher than that of the remaining genes in the Pgt genome (0.86).

To identify the regions of the Pgt genome missing in the SCCL reference assembly, de novo genome assemblies were produced for each isolate. These assemblies were initially produced with the CLC Bio program (QIAGEN) using the paired-end Illumina reads. The CLC Bio contigs were further extended using the 145,823 and 264,876 PacBio reads generated for MCCFC and RKQQC, respectively. To remove contaminating sequences, the genome assemblies were compared against the NCBI non-redundant nucleotide database, and contigs that showed nucleotide similarity to fungal sequences were retained (see Materials and Methods). The final MCCFC and RKQQC genome assemblies measured 93.3 Mb (N50 7,133 bp) and 107.3 Mb (N50 6,292 bp), respectively (Additional file 1: Table S3), and were comparable with the previously reported Pgt genome sizes [16, 41]. The Pgt genome assemblies and their annotations can be downloaded from the project website http://129.130.90.211/rustgenomics/Download.

Pgt isolates show diversity in gene content

To further assess the level of gene conservation between the two sequenced Pgt isolates and at the same time obtain gene expression data for host and pathogen, RNA-seq analysis was performed on the infected wheat leaf tissues. Total RNA was extracted from three biological replicates at each of 5 time points during infection (12, 24, 48, 72, and 96 h post-inoculation (HPI)) for each isolate separately, as well as three replicates of a mock-inoculated control (0 HPI). RNA-seq libraries were constructed for each of the 33 RNA samples and sequenced using Illumina HiSeq 2500, generating between 36 and 50 million quality-trimmed single 100-bp reads for each RNA-seq library (Additional file 1: Table S1).

To assess the level of host gene expression, all reads were first aligned to the reference genome of wheat [32]. To assure that RNA-seq reads are mapped to the correct copies of genes within the polyploid genome, the parameters of the Tophat aligner [37] were optimized using a subset of reads mapping to a homoeologous set of 100 genes present in single copy in each of the wheat genomes. The alignment parameters were selected to maximize the proportion of reads mapping uniquely to each of the duplicated homoeologs. The final analyses were performed using the program settings that allowed for mapping RNA-seq reads in the training dataset with an error rate of less than 1% (Additional file 1: Figure S1).

RNA-seq reads that did not map to the wheat genes, and therefore enriched for fungal sequences, were combined to perform a de novo transcript assembly using Trinity v2.0.6 [33]. The combined de novo assembled transcripts together with the transcripts predicted in the SCCL genome [41] were used to annotate the MCCFC and RKQQC genomes using the PASA pipeline v2.0.1 [42] (Additional file 1: Table S3). In total 18,166 and 18,777 gene models were annotated within the MCCFC and RKQQC genomes, respectively. The predicted open reading frames from the annotated genes were used to extract 16,716 and 16,253 complete proteins from the MCCFC and RKQQC genomes, respectively (Additional file 1: Table S3).

The proteomes of both isolates were compared with the 15,979 publically available proteins predicted in the SCCL genome and the 15,685 proteins predicted in the

Puccinia triticina

(

Ptt

) genome [

41

,

43

]. These proteomes were used in an all-against-all BLASTP comparison followed by protein clustering using OrthoMCL to identify the groups of orthologous and paralogous genes (Table 

4

) [

34

]. In total, 64,633 proteins from the four genomes were clustered into 13,785 groups containing varying number of orthologous and paralogous proteins (Fig. 

1

). Out of these groups, 1,086 were composed of proteins from only one genome, which with the 7,962 ungrouped proteins results in 11,190 proteins that were unique to one of the four sequenced fungal genomes (Table 

4

). These unique proteins included 1,726 and 1,703 proteins from the MCCFC and RKQQC genomes, respectively (Table 

4

). Additionally, 1,428 groups were found to contain proteins from only the MCCFC and RKQQC genomes, and respectively included 1,633 and 1,619 proteins (Fig. 

1

). In total, 3,359 and 3,322 novel proteins were identified in the MCCFC and RKQQC genomes, respectively. The sequences of these proteins were compared to the proteins predicted in the

Pgt

PANaus pan-genome assembly compiled from five Australian

Pgt

isolates  [

16

]. In total, 2,594 (78%) of the novel RKQQC proteins and 2,617 (77.9%) of the novel MCCFC proteins showed sequence similarity to the proteins in the PANaus dataset [

16

]. The identification of homologous sequences in these independently sequenced

Pgt

isolates suggests that the isolate-specific genes identified in the RKQQC and MCCFC genomes are unlikely to be the result of contamination from other DNA sources, but rather represent true presence-absence polymorphisms among the

Pgt

isolates. Additionally, 728 proteins in the RKQQC genome and 742 proteins in the MCCFC genome showed no significant similarity to any

Pgt

proteins and likely represent unique genes within these newly sequenced

Pgt

isolates. Out of these novel proteins, 626 from the RKQQC genome and 625 from the MCCFC genome showed similarity to the known PFAM domains (Additional file

3

: Table S4).

Table 4

Summary of OrthoMCL protein clustering results

Number of proteins used for clustering

64,633

15,979

16,716

16,253

15,685

Number of ungrouped proteins (unique proteins)

7,962

2,158

1,365

1,388

3,051

Number of proteins clustered into ortholog/paralog groups

56,671

13,821

15,351

14,865

12,634

Number of ortholog/paralog groups

13,785

11,100

11,067

10,736

7,696

Number of unique (single genome) groups

1,086

141

166

141

638

Number of proteins in unique (single genome) groups

3,228

333

361

315

2,219

Maximum number of proteins/group

260

107

61

69

106

Total number of unique proteins (unique groups + ungrouped proteins)

11,190

2,491

1,726

1,703

5,270

Fig. 1

Venn diagram displays the results of OrthoMCL clustering. The analysis included 15,979, 16,716, and 16,253 proteins predicted in the SCCL, MCCFC, and RKQQC genomes, respectively. In addition, OrthoMCL clustering included 15,685 proteins predicted in the Puccinia triticina f.sp. tritici (Ptt) genome. Numbers represent the counts of orthologous/paralogous protein groups. Groups composed exclusively of proteins from the RKQQC (a) and MCCFC genome (b) included 315 and 361 proteins, respectively

An alternative approach to assessing the gene content divergence between the newly sequenced Pgt isolates was based on the combined comparative analysis of genomes and RNA-seq data. For this purpose, the RNA-seq data from both isolates were separately aligned against the SCCL reference transcripts. RNA-seq reads from the RKQQC and MCCFC infected leaf tissues were successfully aligned to 10,154 and 10,437 SCCL genome gene models, respectively, and individual transcript abundances were calculated for each isolate. A conservative approach was taken to qualify genes as present or absent in the sequenced genomes. A gene annotated in the SCCL genome was considered present in a newly sequenced isolate if it showed expression (FPKM > 0) within that isolate and/or if it resided within an orthologous group with a protein from that isolate. For a gene annotated in the SCCL genome to be considered absent from the genomes of RKQQC or MCCFC it would have to lack any detectable expression from that isolate, and have no detected orthologs from that isolate within the OrthoMCL dataset. In total, 810 genes (5.1%) were found to be present in the RKQQC genome and absent in the MCCFC genome. Conversely, 825 SCCL genes (5.2%) were present in the MCCFC genome and absent in the RKQQC genome. A total of 1,509 genes (9.4%) annotated in the SCCL genome were not detected in the transcript data from either isolate. The high degree of divergence in the gene contents of RKQQC, MCCFC, and SCCL genomes is consistent with the degree of divergence observed in other Pgt isolates as well as between isolates of other fungal pathogens [14, 16, 44]. Expectedly, the highest level of divergence in the gene content was found between the Ptt genome and the genomes of Pgt isolates.

MCCFC and RKQQC Pgt races have distinct effector complements

Biotrophic fungal pathogens interact with and manipulate their host plants through the use of secreted effector proteins. For this study, several criteria were used to identify and partition the likely effector candidates present in the genomes of the two

Pgt

isolates (Fig. 

2

). All 15,979 proteins predicted in the SCCL genome were screened for the presence of N-terminal signal peptide [

45

] and the absence of transmembrane domains outside the signal peptide region [

35

]. As described above, the OrthoMCL protein family information along with the RNA-seq data were used to identify those effector candidates that were present, and those that were absent from each isolate. Of the 1,799 secreted proteins identified in the SCCL genome, 1,586 were found present in both RKQQC and MCCFC, while 68 and 80 were found exclusively in either RKQQC or MCCFC, respectively, and 65 candidates were not detected in either isolate (Additional file

2

: Table S2).

Fig. 2

Diagram showing the pipeline for candidate effector classification into three groups. The SCCL gene models were compared with the genomic and transcriptomic sequence data generated for MCCFC and RKQQC. Effector candidates were categorized as unique, conserved polymorphic, or perfectly conserved between the two isolates

Using variant calling data, at least one non-synonymous mutation was found in 787 candidate effector-encoding genes. Based on these data the effector candidates were categorized into three groups (Additional file 2: Table S2). The first was composed of 148 candidates that were exclusively found in either the MCCFC or RKQQC genomes. Although it appears that these proteins are ostensibly dispensable for infection of the wheat cultivar Morocco, they are still high priority candidates because this type of presence/absence variation is the most likely to have an effect on the Pgt virulence. The second group of candidate effectors is composed of 728 genes that are conserved in both the MCCFC and RKQQC genomes yet show non-synonymous coding sequence variation between the isolates. These polymorphisms have the potential to alter the function of specific proteins and affect the interaction between the pathogen and host. The remaining 858 genes were conserved between the MCCFC, RKQQC, and SCCL genomes, containing only synonymous SNP variation in the coding sequences. These conserved genes likely represent a core group of effector candidates enriched for effectors that may be essential for successful fungal infection.

Co-regulation of host and pathogen transcriptomes during the course of infection

To better understand the effects of effector complement divergence and conservation between the Pgt isolates on host gene co-regulation in the rust-wheat pathosystem, the time-resolved RNA-seq expression profiles of infected leaf tissues were generated. The RNA-seq reads mapped to the publically available wheat [32] and Pgt [41] reference genomes were used to obtain the FPKM values for their respective gene models at each of the six infection time points 0, 12, 24, 48, 72, and 96 HPI (NCBI GEO accession number GSE93015). The FPKM values were used to compare the joint expression of fungal (except 0 HPI) and wheat genes between different isolate treatments at the same time points, as well as across the time points of the same isolate treatment (Additional file 1: Figure S2). Based on the relative proportion of fungal reads mapped to the reference genome, both Pgt isolates demonstrated quite similar temporal patters of gene expression with the RKQQC race showing the slightly reduced proportion of mapped reads at 72 and 96 HPI (Additional file 1: Figure S3). However, the difference in the fraction of mapped reads at the 96 HPI time-point was not statistically significant.

In total 1,238

Pgt

genes and 6,750 

T. aestivum

genes were identified as differentially expressed between the two isolate treatments and/or across the time course of infection for at least one isolate (FDR ≤ 0.05, log

2

-fold expression change ≥ 2) (Table 

5

). For the detailed analyses of expression profiles, the data was further filtered to retain only those genes that have data available for at least six biological replicates in the entire dataset, resulting in 1,054

Pgt

genes and 3,877 wheat genes (Additional file

4

: Tables S5 and S6). The pair-wise correlation analysis between the expression profiles of each fungal gene in the MCCFC and RKQQC datasets showed that the majority of genes (57%) have similar expression patterns during the course of infection with the Pearson correlation coefficient (PCC) above 0.5 (Additional file

1

: Figure S4 and Additional file

5

: Table S7). Comparison of the wheat expression profiles between the MCCFC and RKQQC datasets revealed that the substantial fraction of genes (71%) show similar patterns with the PCC > 0.5 (Additional file

1

: Figure S4 and Additional file

5

: Table S8).

Table 5

Wheat and Pgt genes differentially expressed between the Pgt isolates or across different time-points

P. graminis

317

308

257

533

571

153

707

1,238

T. aestivum

167

182

177

419

703

0

6562

6750

Further, we have functionally classified genes shared between the MCCFC and RKQQC datasets using the GO terms and assessed the average PCC for genes within each GO group (Additional file 5: Table S9 and S10). In Pgt, the most similar expression profiles were obtained for genes involved protein biosynthesis (GO:0043043, GO:0005840) and various metabolic processes (GO:1901135, GO:0019637), while the genes involved in the regulation of transcription (GO:1903506, GO:0003677) and RNA biosynthesis (GO:2001141) showed the lowest level of gene expression correlation. Among the wheat genes, those that were involved in protein biosynthesis (GO:0032544, GO:0008135, GO:0006412), stress response (GO:0009409, GO:0048583), transcription factor activity (GO:0001071) showed the high level of gene expression correlation between the MCCFC and RKQQC datasets. The lowest correlation values were found for genes involved in lipid transport (GO:0006869) and cellulose metabolism (GO:0030243).

To obtain more detailed picture of the complex transcriptional events that occur in both plant and pathogen over the course of infection, we performed k-means gene clustering (Fig. 

3

) for each isolate-specific dataset (Additional file

6

: Table S11). A number of gene clusters containing both the

Pgt

and wheat genes have been identified suggesting the connectedness of host’s and pathogen’s regulatory processes. The biological significance of each gene cluster was assessed by performing the GO term enrichment analysis (Additional file

6

: Tables S12-S17). We have selected top 10 and 9 clusters for the wheat genes in the RKQQC and MCCFC datasets, respectively (Fig. 

3

). The host’s gene clusters in plants infected with different

Pgt

isolates were enriched for similar GO terms likely indicating the significance of associated pathways for plant-pathogen interaction. Among the identified clusters there were genes known for their involvement in response to the pathogen infection such as salicylic acid-dependent systemic acquired resistance [

46

], response to chitin [

47

], defense response to fungus, and the regulation of immune response (Fig. 

3

). The GO-enriched clusters also included genes involved in photorespiration and chloroplast organization that play critical role in the outcome of fungal infection [

48

]. Similar to the results of the transcriptome profiling of infected wheat tissues [

24

], the stress response pathways associated with salicylic acid- and jasmonic acid-induced metabolic processes were also over-represented in the wheat clusters.

Fig. 3

Clusters of co-expressed Pgt and wheat genes. Co-expressed gene clusters were obtained by k-means clustering of the expression data generated by the RNA-seq profiling of the MCCFC- and RKQQC-infected leaves. The GO term enrichment levels for each gene cluster are shown on the heat maps as the log2 fold changes compared to the genome-wide estimates

The pathogen clusters generated for both the MCCFC and RKQQC datasets were enriched for genes with the oxidoreductase and antioxidant activities (Additional file 6: Tables S12-S17). Similar to previous study of the stripe rust infected wheat leaves [24], we also found over-representation of the genes involved in nucleic acid, protein and polysaccharide metabolism. The development of fungi was also associated with the increase in the transmembrane transport (Fig. 3) that was also demonstrated for Fusarium oxysporum during the colonization of Medicago truncatula host [49]. We have identified 7 clusters in the RKQQC dataset and 8 clusters in the MCCFC dataset enriched for the fungal genes encoding candidate effectors with the secretion signal peptides (Fig. 3, Additional file 6: Tables S14 and S17), consistent with the role of effectors in manipulating host’s responses to establish compatible interaction.

Comparative analysis of GCNs identifies conserved and divergent sub-networks around specific gene ontologies

GCNs have been shown to be a powerful tool to characterize system level interactions between a host and a pathogen [50]. As with the K-means clustering, the Pgt isolate-specific GCNs were constructed using the FPKM expression values obtained for plant and pathogen genes that showed statistically significant changes over the experiment. To account for the non-scale free nature of the time-course RNA-seq data, a Graphical Gaussian model implemented in the R package GeneNet was used to calculate partial correlations between the expression profiles of each gene within the respective isolate treatment [38, 51]. Significant partial correlations (network edges) between genes (network nodes) were selected using the FDR cutoff value of 0.001 (Additional file 7: Tables S18 and S19).

The RKQQC and MCCFC race treatment networks contained 2,811 and 2,843 significantly connected nodes from both host and pathogen, linked respectively by 87,975 and 141,343 edges. The direct comparison of all significant edges between the isolate-specific networks revealed that 18,023 edges are conserved representing 20.5% of the RKQQC and 12.8% of the MCCFC treatment edges. Interestingly, the edge conservation was not equally distributed across the networks. The edges connecting certain groups of nodes show a much higher degree of network conservation than others. This variation in edge conservation is most apparent when dividing edges by the species-of-origin of nodes they connect. Between the two isolate-specific networks, edges linking wheat nodes to other wheat nodes were 13.8% conserved, edges connecting

Pgt

to

Pgt

nodes were 2.9% conserved, and those linking

Pgt

to wheat nodes were only 1.7% conserved (Table 

6

).

Table 6

The level of edge and node conservation between the isolate-specific GCNs

Conserved

18,023

957

1,131

15,935

1,471

1,062

409

5.3%

6.3%

88.4%

72.2%

27.8%

Unique to RKQQC

69,952

17,215

19,491

33,246

2,788

1,875

913

24.6%

27.9%

47.5%

67.3%

32.7%

Unique to MCCFC

123,320

38,763

17,939

66,618

2,842

2,031

811

31.4%

14.6%

54.0%

71.5%

28.5%

Percent conservation

8.5%

1.7%

2.9%

13.8%

20.7%

21.4%

19.2%

To understand which biological processes are affected by each

Pgt

isolate in the susceptible host, all wheat genes within each isolate-specific GCN were functionally annotated using the Gene Ontology (GO) terms [

52

] (Fig. 

4a

). Significantly enriched GO terms were identified for three types of nodes in the GCNs: nodes connected by the RKQQC-specific edges, nodes connected by the MCCFC-specific edges, and the nodes connected by the conserved edges found in both GCNs. In total 82 GO terms were significantly enriched (Fisher exact test,

p

-value 

1

: Table S20). After removing redundancy inherent in the GO term hierarchy, this list was winnowed down to 51 groups of GO annotated wheat genes. To identify network modules that are associated with each enriched GO term, the conserved GO annotated genes were used to extract sub-networks from each isolate (henceforth, GO sub-networks). These sub-networks contained first order-connected nodes originating from wheat and

Pgt

including among other fungal genes the

Pgt

effector candidates. To identify the nodes and edges that were conserved or unique between the isolate-specific networks, combined sub-networks containing nodes from both isolate-specific GCNs were created by merging the GO sub-networks from each isolate based on the conserved nodes and edges (Fig. 

4b

).

Fig. 4

Comparison of Pgt-wheat gene co-expression networks. a A diagram of the analysis pipeline used to generate the combined GO sub-networks from each isolate-specific GCN. b. The example of combined network composed of three nodes the are conserved between both networks, and two nodes that are unique to the RKQQC- and MCCFC-specific networks. In this example, edges a, b, and c are first order edges for node CON1, while d, e, and f are second order edges for the same node. In this example, only edge a is conserved between the two isolate-specific networks. A sub-network for node CON1 contains nodes CON2, RKQQC, and MCCFC, but excludes CON3 because it doesn’t have a first order edge connecting it to CON1. The edges of the CON1 sub-network include edges a, b, c, and d, but exclude e and f because they connect with a node outside of the sub-network

Comparisons of the combined sub-networks revealed that they vary widely in terms of the proportion of nodes and edges that are conserved (Fig. 

5

) (Additional file

1

: Table S21). Some combined GO sub-networks showed a high degree of node and edge conservation, while others had mostly unique sets of nodes and edges (Fig. 

5a

). The proportion of conserved edges ranged from 0 to 26%, and the proportion of node conservation in the GO sub-networks ranged from 1.3 to 55%. Randomized networks, created by randomly regenerating all the edges of each isolate network yet preserving the same number of nodes and the total number of edges, showed significantly lower level of node/edge conservation pattern (Fig. 

5b

). The level of node and edge conservation was assessed for each of one thousand sub-networks that were created by randomly choosing sets of 5-nodes from both the randomized and the real-world isolate networks (Fig. 

5b and c

). The results of the permutations revealed that many of the GO sub-networks within the experimental dataset show a significantly higher degree of edge conservation than is expected between randomly generated networks of equal size and degree of connectedness. Furthermore, this edge conservation varied widely across the real-world sub-networks, with some sub-networks showing a high level of conservation and others showing little or no conservation (Fig. 

5a

). The observed heterogeneity in node/edge conservation across the real-world GCNs suggests that while certain host plant pathways can be modulated in a similar manner by both

Pgt

isolates, likely utilizing the virulence factors (effectors) shared by both

Pgt

isolates, other pathways might be affected by the virulence factors differentiating one

Pgt

isolate from another (Additional file

1

: Table S21).

Fig. 5

Wheat GO sub-networks vary in the degree of network conservation between the Pgt isolate treatments. a Examples of four different GO sub-networks, each including five wheat GO annotated genes, show a high degree of variability in both edge and node conservation between the two Pgt isolate treatments. b Each GO sub-network was randomly regenerated using the same number of nodes and edges observed in our experimentally reconstructed GCNs. For this purpose, sub-networks including five wheat nodes were randomly sampled 1,000 times, and produced a normal distribution of edge conservation values with a mean of 0.18%. c Five node sub-networks were randomly sampled from the real world data, and produced a multi-modal distribution of edge conservation values ranging from 0.05 to 24.81%

The Pgt effector candidates within the networks showed conservation in terms of the GO sub-networks they were associated with. Although the overall edge conservation between network conserved fungal effector candidates and wheat genes was relatively low (2.3%), the rate of network conserved effectors being associated with the wheat genes in the same GO term, on average, was 13.7%, and was as high as 32.7% for some candidates. Indeed, of the 126 Pgt effector candidates present in both isolate-specific networks, 100 (79.4%) were associated with the same GO sub-network in each isolate.

To investigate whether the degree of GO sub-network conservation is related to the sequence level conservation of effector candidates between the isolates, we compared the proportion conserved edges within each sub-network with the degree of DNA sequence conservation (as described in Fig. 

2

) for all effector candidates associated with each sub-network. We found that there is a significant positive correlation between the proportion of sequence conserved effector candidates associated with a specific GO sub-network and the level of edge conservation in that sub-network (Fig. 

6

) (R

2

 = 0.2). These results suggest that GO sub-networks with the higher levels of edge and node conservation between networks also show tendency to be associated with a higher proportion of effectors showing the high levels of sequence conservation. This trend has lead us to hypothesize that these effector candidates may be directly involved in modulating these biological pathways within the host plant, and that the degree of sequence divergence in the effector complements likely influences the regulation of host pathways.

Fig. 6

The relationship between GO sub-network edge conservation and proportion of perfectly conserved effectors. The percent edge conservation for each GO sub-network plotted against the percentage of perfectly conserved at the sequence level (as determined by genomic and transcriptomic sequence comparisons) effector candidates associated with that sub-network. Geometric point size reflects the total number of effector candidates associated with each sub-network. The best-fit regression line (blue) and standard error intervals (dark grey) are shown. Those sub-networks with a higher percentage of conserved effectors, particularly those with few effector candidates overall tend to have a higher degree of sub-network conservation

RKQQC and MCCFC isolate treatments produce distinct co-expression networks around salicylic acid response genes

The observed variation in the proportion of conserved host-pathogen edges among different GO sub-networks suggest the existence of convergent and divergent modes of evolution among the host-pathogen regulatory modules, where some modules are regulated by a conserved set of effectors and some are regulated by a divergent set of effectors. The latter examples include quite intriguing cases of the GO sub-networks, which are affected by distinct sets of effectors from different Pgt isolates.

Here, we have performed a more detailed analysis of one set of genes that produce a distinct sub-network structure within each isolate and include the salicylic acid (SA) responsive genes (GO:0009751) (Fig. 

7a

). Despite the fact that all of the five SA-responsive genes themselves showed similar expression profiles between the isolate treatments (Fig. 

7c

), relatively few of the other nodes and edges, including the effector-encoding genes, within the SA sub-network were conserved between the isolates (15.8% node conservation and 2.7% edge conservation (Fig. 

7a

, Additional file

1

: Table S21). In other words, it seems that very different systematic changes caused by each isolate treatment can produce very similar results with regard to these SA-responsive genes. The low level of edge/node conservation suggests that the SA response pathway is modulated differently by each

Pgt

isolate, perhaps by deploying different sets of effectors (Fig. 

7b

). A total of 27

Pgt

effector candidates were associated with the SA co-expression sub-network, with only 2 candidates (PGTG_18238 and PGTG_02185) consistently associated with the sub-network in both

Pgt

isolates.

Fig. 7

A comparison of the GO:0009751 co-expression sub-networks from each isolate. a The conserved (grey) nodes and edges, versus the nodes and edges that are unique to the RKQQC treatment (red) and MCCFC treatment (blue) are shown. The GO annotated genes themselves are highlighted as yellow nodes. b The candidate effector nodes within both sub-networks (orange squares), and their connections (orange edges) with the host plant nodes (green) are shown. c The log2 FPKM expression values of the 5 genes from GO:0009751, which do not show significant difference in expression under each isolate treatment. d Relative transcript abundance T. aestivum genes annotated as (GO:009751) responsive to salicylic acid. Transcript specific primers were used to perform qRT-PCR on mRNA extracted from wheat lines 18 h after salicylic acid (SA), abscisic acid (ABA), or mock treatments

The five GO:0009751 annotated genes are known for their association with the plant-fungal interaction. Three of the genes are homologous to the MYB transcription factors (Traes_5AS_7D519210E.5, Traes_5BS_D1C03C165.1, Traes_5DS_1EF547639.5) and thus have been designated MYBa, MYBb, and MYBd, respectively. All three of these genes showed protein sequence similarity with MYB59 (AT5G59780) and MYB48 (AT3G46130) from Arabidopsis. In Arabidopsis, these two transcription factors were shown to be specifically up regulated in response to SA treatment [53]. Additionally, the expression of AT5G59780 has also been shown to be responsive to chitin treatment [54] indicating that this gene is likely involved in the recognition of or response to pathogenic fungi.

The presence of these MYB transcription factors in the SA-responsive sub-network raised the possibility that they are involved in the co-regulation of gene expression within this sub-network. To investigate this possibility further we used Nsite program [55] to predict putative transcription factor binding sites within the promoter regions 3 Kb upstream of each wheat gene from the SA-responsive sub-network. A total of 50 MYB-responsive elements have been identified in the promoter regions of genes within the sub-network, a significant enrichment relative to the genes in the GCNs as a whole (Fisher’s exact test, P-value = 0.017) (Additional file 8: Table S22).

The two other GO:0009751 annotated genes Traes_2AL_6A8D574C4.1 and Traes_2AL_52C4B6996.2 showed amino acid sequence similarity to caleosins and peroxygenases from other plant species and were designate POX1 and POX2, respectively (Fig. 7a). Peroxygenases are part of the oxylipin metabolic pathway, which is known to generate anti-fungal compounds [56, 57]. RD20, a caleosin/peroxygenase from Arabidopsis, has been shown to respond to salicylic acid and increase plant resistance to pathogenic fungi via the oxylipin pathway [58, 59]. Furthermore, oxylipin metabolism has also been shown to play an important role in the fungal pathogenesis of wheat. Lipoxygenases, another group of enzymes involved in the oxylipin pathway, were recently shown to function in the wheat defense response against Fusarium graminearum [60]. Therefore, the finding of 7 transcripts (Ta2alMIPSv2Loc183543.1, Traes_2BL_77148B8D8.1, Traes_2DL_CE85DC5C0.1, Traes_2DL_B5B62EE11.2, Traes_5BS_060785740.2, Traes_5DS_E8892706A.2, and Traes_6DS_E66547E66.2) with strong similarity to lipoxygenases from wheat and other grasses in the combined SA-responsive sub-network provides further support in the validity of the constructed GCN, and suggests that the oxylipin metabolic pathway possibly plays important roles in the interaction between wheat and Pgt.

To confirm responsiveness, either positive or negative, of these genes to SA in wheat, their expression levels were assessed in the wheat seedlings of cultivar Morocco treated with the solutions of SA. Using qRT-PCR we compared the expression for four of the GO:0009751 annotated genes after SA treatment against mock and abscisic acid-treated controls (Fig. 7d). All four genes showed a significant change in expression, either positive or negative, in response to SA treatment (P-value < 0.05).

These data are a further indication that the homology-based GO annotations for these genes likely reflect true functional roles in the wheat SA response pathway. The SA mediated repression of the three wheat MYB transcription factors is consistent with the transcriptional changes observed in a minority of Arabidopsis MYBs after treatment with SA [53]. Conversely, the consistent and slightly increased expression of these MYBs over the course of both compatible interactions seems to support their potential for having roles in the plant defense response against Pgt.

About

We Support OUR Contributors

Get Our Newsletter

 Receive podcast updates
Exclusive insights
Patient Engagement Tips from industry experts
We hate SPAM as much as you do and promise to keep your email address safe.
  • Subscribe to the Podcast