Comparative expression analysis of Phytophthora sojae polysaccharide lyase family 3 (pectate lyase) genes during infection of the soybean Glycine max

Abstract

Phytophthora sojae is an important plant pathogen affecting soybean crops worldwide. The specific mechanisms this pathogen uses to penetrate the host cell and initiate infection have not been fully elucidated. However, a greater understanding of the molecular mechanisms that facilitate these processes will be critical to the development of novel strategies to control this devastating agricultural pathogen. It is likely that, to initiate infection, P. sojae uses an appressorium-like structure to break through the plant cell wall. It is also presumed that the pathogen relies on enzymatic activity to weaken the cell wall. The polysaccharide lyase (PL) enzyme superfamily, particularly the PL family 3 (PL3), has been hypothesized to play an important role in the process of host cell penetration by contributing to the degradation of the structural polysaccharides of the cell wall. To investigate this hypothesis, we have scanned the revised version of the annotated P. sojae genome for the presence of putative PL3-coding genes and conducted an extensive sequence analysis of all gene models found. In addition, we have quantified the relative expression of each gene in P. sojae during infection of susceptible (Williams) and resistant (Williams 82) soybean cultivars over a 48 h time period. Twelve PL3-coding gene models were identified and initial results indicate that during infection of the Williams cultivar, several of these genes experience significant up-regulation during the first 48 h of infection. An increase in transcriptional activity is also observed during infection of Williams 82; however, this expression pattern is distinct from that observed in the Williams cultivar. These results provide evidence regarding the potential involvement of PL-coding genes in the early pathogenic processes of P. sojae and suggest that this pathogen expresses genes differentially, depending on whether infection is occurring on a resistant or a susceptible soybean cultivar.

Background

Soybeans (Glycine max) are the United States’ second largest cash crop and the largest agricultural export, with a total yearly production value of $41 billion for the 2017 season (https://www.nass.usda.gov/Statistics_by_Subject/result.php?4378EB63-5283-3A7B-89EF-63C62203D62B&sector=CROPS&group=FIELD%20CROPS&comm=SOYBEANS). Glycine max has been described as a hardy crop, well adapted to a variety of soils and conditions (Jin et al. 2003). However, several diseases, including infection with the pathogenic oomycete Phytophthora sojae, result in a devastating reduction in crop health and production (Gally et al. 2007). Total losses worldwide are estimated to reach $1–2 billion annually (Tyler 2007). These figures have grown over the past decade due to widespread use of no-till practices, soybean monoculture, and genetic uniformity of cultivars (Gally et al. 2007), all of which increase the likelihood of infection.

Phytophthora sojae is a soil-borne plant pathogen of the phylum Oomycota. Oomycetes, originally classified as a group of fungi, were later discovered to be phylogenetically closer to brown algae and diatoms (Qutob et al. 2000; Tyler et al. 2006). Currently, oomycetes are considered members of the Kingdom Chromista, which is sometimes referred to as Stramenopila (Cavalier-Smith 1986; Ruggiero et al. 2015; Cavalier-Smith 2018). Although many Phytophthora species have the ability to infect several host species, P. sojae primarily infects soybean plants, resulting in necrotic stem and root rot (Qutob et al. 2000).

Phytophthora sojae produces single-celled, motile zoospores that use positive chemotaxis towards soybean root-derived compounds (Morris et al. 1998). Once near a host, these cells encyst and germinate on the soybean root or hypocotyl surface, resulting in the formation of a germ tube. Some reports indicate that the germ tube may then enlarge to form an appressorium, an extended portion of hyphal tissue that participates in host tissue penetration (Ehrlich and Ehrlich 1966; Morris et al. 1998; Schmitthenner 1999; Li et al. 2010), although others have failed to find conclusive evidence of the formation of such a structure in P. sojae (Enkerli et al. 1997). As the disease progresses, P. sojae continues to spread through the intercellular matrix of the plant, forming extensions termed haustoria (Coffey and Wilson 1983; Enkerli et al. 1997), which absorb nutrients from the host plant cells via intimate contact with their plasma membrane (Kamoun 2003; Panstruga and Dodds 2009). The association of haustoria with the host cell wall, and the extrahaustorial membrane (which is derived from the host plasma membrane) provide physical separation between the pathogen and the host (Panstruga and Dodds 2009). This three-membrane assembly serves as a conduit for the free flow of pathogenic proteins and degrading enzymes from the pathogen to the host (Ellis et al. 2006).

Current research on P. sojae has focused on elucidating the specific pathogenic determinants that facilitate penetration of the host. Research from our laboratory has centered on genes encoding enzymes involved in carbohydrate metabolism such as carbohydrate esterases (CE), glycoside hydrolases (GH), glycosyl transferases (GT), and polysaccharide lyases (PL), collectively known as carbohydrate-active enzymes or CAZymes (Ospina-Giraldo et al. 2010a; Ospina-Giraldo et al. 2010b; Mingora et al. 2014; Horowitz and Ospina-Giraldo 2015). Initial annotation analysis of the genome assembly (version 1.1) indicated that the genome of P. sojae contains a total of 379 CAZy homologs. With the exception of the genes encoding glycosyl transferases, these major gene superfamilies encode enzymes referred to as cell wall degrading enzymes, or CWDE (Ospina-Giraldo et al. 2010a). CAZyme orthologs have been found in the genomes of all oomycete species studied thus far, and in many fungal genera, including Fusarium, Magnaporthe, and Neurospora, among others (Ospina-Giraldo et al. 2010a; Horowitz and Ospina-Giraldo 2015; Hinkel and Ospina-Giraldo 2017).

The cell wall, a structure composed largely of carbohydrate-derived compounds including pectin, cellulose, and various other polysaccharides, constitutes the primary barrier that must be breached by a pathogen in order to successfully colonize the plant. However, evidence suggests that, in Phytophthora spp., the force generated by the appressorium alone is insufficient to breach the rigidity and strength of the plant cell wall (Money et al. 2004). Therefore, it is likely that the enzymatic degradation of this barrier may significantly contribute to successful pathogenesis, and occurs by the direct action of CAZymes encoded by the numerous CAZy homologs contained within the P. sojae genome.

Pectin is one of the most abundant components of the cell wall and can account for up to 35% of cell wall dry mass (Rogers and Perkins 1968). It consists of as many as 17 distinct monosaccharides, with a significant percentage of galacturonic acid (Pelloux et al. 2007), which can form linear polymers joined by α-(1,4)-glycosidic bonds, often modified at the C-6 carboxyl end by the addition of a methyl group via an ester bond. A major role in the degradation of pectin is played by enzymes of the PL superfamily, in particular PL family 3 (PL3) members (also known as pectate lyases). The primary function of this enzyme family is to hydrolyze uronic acid-containing polysaccharide chains via a β-elimination mechanism, leading to the release of an unsaturated hexenuronic acid residue and a new reducing end (Yip and Withers 2006; Cantarel et al. 2009). PL3 enzymes thus randomly cleave (1 → 4)-α-D-galacturonan into α,β-unsaturated reaction products (Edstrom and Phaff 1964). Genes coding for PL3 enzymes are among the most represented CWDE genes in the P. sojae genome (Ospina-Giraldo et al. 2010a). Further analyses of their expression and secretion may serve to expand current understanding of P. sojae pathogenicity and may provide information critical to the development of novel infection control methods.

Results

The putative PL family 3 gene complement in Phytophthora

A previously compiled list of P. sojae PL3-coding genes (Ospina-Giraldo et al. 2010a), based on the 1.1 version of the annotated genome released by the JGI-DOE, was reevaluated to create an updated list based on the more recent version (3.0) of the P. sojae genome assembly. These sequences were further compared with NCBI’s GenBank database using BLAST, and accession numbers were obtained for each gene. A new grouping of PL3-coding genes was established, resulting in a set of 12 putative gene models (Table 1), all of which were confirmed to contain the PFAM domain PF03211, characteristic of the pectate lyase family (Fig. 1). Further sequence comparison analysis revealed that two of the gene models had 100% sequence similarity with each other, which prevented their distinct evaluation via qPCR techniques. HMMER analysis of the final PL3-coding gene set revealed homology to pectate lyases from P. infestans, P. sojae, P. nicotianae, P. parasitica, P. palmivora, and P. kernoviae, in addition to several Pythium and Colletotrichum spp. All hits had E-values of e^-100 or lower.

Table 1 The updated PL3-coding gene complement in P. sojae genome

Fig. 1

HMM logo of Pfam03211 consensus sequence. Hidden Markov Model analysis revealed the presence of Pfam03211 domain in all PL3 sequences. The Pfam03211 Pectate Lyase Family 3 domain is 96 amino acids long and is typically located closer to the amino terminus of pectate lyase proteins. Analysis was conducted using the EMBL-EBI Pfam server (http://pfam.xfam.org/), using HMMER (Schuster-Böckler et al. 2004)

Eight out of the 12 PL3-coding gene models identified on the P. sojae genome assembly are distributed in two clusters of three and five sequences, respectively, located on scaffold 3. Two more sequences are closely located on scaffold 13 and the remaining two are found on scaffolds 2 and 6. One of the gene models (XM_009526560.1) is fully enclosed within a much larger ORF in the complementary strand, whose hypothetical function has not been described.

Nine gene models consisted of a continuous exon. Conversely, three gene models (those encoding transcripts XM_0095239513.1, XM_009526562.1, and XM_009526401.1) appeared to contain introns of 22, 213, and 531 bp, respectively. The putative 22 bp intron in gene model encoding transcript XM_0095239513.1 is predicted to be located between positions 280 and 301. Since introns (particularly those larger than 90 bp) are rare in Phytophthora spp., we decided to further investigate these three sequences. Sequence analysis of the gene models encoding transcripts XM_009526562.1, and XM_009526401.1 revealed that both of them contained one single, complete, and uninterrupted open reading frame. Furthermore, PCR and RT-PCR amplification using primers targeting regions either adjacent to the putative introns or within the intron itself produced amplicons of the same size (Additional file 1: Figure S1). We sequenced these amplicons and confirmed that, for each of the two genes, the sequences of the fragments obtained by PCR and RT-PCR are identical to each other, and are also identical to the corresponding gene model sequences available at JGI. Analysis of the third gene model, XM_0095239513.1, indicated that the only stop codon present in this sequence is actually located 41 bp downstream from the end (position 301) of the originally predicted putative intron. Therefore, we reanalyzed this sequence using different intron prediction algorithms (Additional file 2: Figure S2) and found that signals for a potential intron exist at a location spanning the only stop codon found in the sequence. This new putative intron is 37 bp long and is located between positions 310 and 346 of the gene model sequence. PCR and RT-PCR amplification using primers targeting sequences immediately adjacent to the putative intron produced fragments differing by ~ approximately 40 bp (Additional file 1: Figure S1). Sequence analysis confirmed that the fragment amplified by RT-PCR specifically lacks the 37 bp intron sequence.

Evolutionary analysis

Several evolutionary trees were constructed based on a CLUSTAL W alignment of the compiled gene sequences using three different phylogeny reconstruction methods. All trees have a relatively similar topology, with all gene sequences contained within three clades, regardless of the inference method used (Fig. 2). The largest of these clades can be subdivided in two subclades, one of them containing gene models encoding transcripts XM_009523714.1, XM_009523712.1, XM_009526401.1, and XM_009526391.1 while the other subclade contains sequences from XM_009526562.1, XM_009526560.1, and XM_009523709.1. A second, smaller clade contains gene models encoding transcripts XM_009539509.1, XM_009539513.1 and XM_009526324.1. The remaining third clade contains sequences from XM_009532098.1, and XM_009522292.1. Bootstrap analysis with 1000 replications yielded seven branch points with 100% confidence. Genes located on the same scaffold were positioned together in the same clade. Furthermore, an estimate of the evolutionary divergence between sequences conducted using the Maximum likelihood model (Tamura et al. 2004) and the bootstrap analyses conducted on all inferred trees support the hypothesis that at least four of the PL3-coding gene models found in P. sojae can be grouped in two pairs of true paralogous genes (Fig. 2). An additional evolutionary analysis involving representative sequences from P. sojae and other Phytophthora spp. confirmed the evolutionary relatedness of PL3-coding genes within the genus and each of the two major clades identified (based on bootstrap values greater than 80%) after the phylogenetic analysis included sequences from all Phytophthora spp. (Fig. 3).

Fig. 2

Maximum Parsimony analysis of P. sojae PL3-coding gene sequences. The evolutionary history was inferred using the Maximum Parsimony method and the consensus tree inferred from two most parsimonious trees is shown. The percentage of parsimonious trees in which the associated taxa clustered together are shown next to the branches. The analysis involved 12 nucleotide sequences. Codon positions included in the analysis were 1st + 2nd + 3rd + Noncoding. There was a total of 674 positions in the final dataset and all positions containing gaps and missing data were eliminated. Tree shown is representative of the topology generated by the Neighbor Joining and Maximum Likelihood methods using the same dataset. Evolutionary analyses were conducted in MEGA7

Fig. 3

Maximum Parsimony analysis of PL3-coding gene sequences from multiple Phytophthora spp. The evolutionary history of 18 Phytophthora spp. PL3-coding gene sequences was inferred using the Maximum Parsimony method and the consensus tree inferred from two most parsimonious trees is shown. The percentage (cut-off value = 80%) of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) is shown next to the branches. There was a total of 656 positions in the final dataset and all positions containing gaps and missing data were eliminated. The scale bar represents the number of nucleotide changes over the whole sequence. Taxa included in this analysis are P. sojae (magenta square), P. infestans (yellow circle), P. ramorum (black circle), P. capsici (cyan diamond), P. cinnamomi (green triangle), and P. parasitica (red circle). No NCBI accession identifiers are available for P. ramorum or P. cinnamomi; therefore, identifiers from each genome project, as shown in FungiDB, were used instead. Tree shown is representative of the topology generated by the Neighbor Joining and Maximum Likelihood methods using the same dataset. Evolutionary analyses were conducted in MEGA7

Promoter and signal peptide analysis

Analysis of the 500 bp sequence upstream of the translational start point using the MEME suite revealed three motifs within the 12 PL3-coding gene promoter regions. Two of the three motifs were scattered and appeared to have no conservation regarding their position in the sequences. However, one motif present in five of the 12 sequences analyzed by MEME appeared to be in a very similar location across the sequences. This motif is 40 nucleotides in length and is found consistently between positions 22–105 upstream of the transcription start site (Fig. 4).

Fig. 4

Shared motif in the promoter sequence region of several P. sojae PL3-coding genes. The motif was identified by MEME-Suite analysis of aligned nucleotide sequences, up to 500 bp upstream of the transcription start site for each PL3-coding gene. Five, out of eleven genes (XM_009532098.1, XM_009522992.1, XM_009523712.1, XM_009526562.1, and XM_009526401.1), contained this conserved motif, and its presence correlated with observed increases in transcriptional activity during pathogenicity. Size of base indicates relative level of conservation in the examined sequences

SignalP 4.1 analysis indicated that eight out of 12 predicted amino acid sequences possess a signal peptide that allows secretion into the extracellular spaces, with a consistent cleavage site between residues 17–22. SecretomeP 2.0 analysis equally revealed that all 12 predicted amino acid sequences have the potential for non-canonical secretion (Additional file 3: Table S1).

Polysaccharide lyase gene expression during infection

Custom primers were designed to be used in the subsequent qPCR analysis of the transcriptional levels of the PL3-coding genes (Additional file 4: Table S2). (Because genes encoding transcripts XM_009523714.1 and XM_009523712.1 are identical, it should be understood that expression results regarding transcript XM_009523712.1 are actually combined values of XM_009523712.1 and XM_009523714.1 transcriptional levels.) These primers were tested in a standard PCR with P. sojae DNA to confirm amplification, and with extracted G. max DNA to confirm the primers did not amplify any plant material.

Expression in cv. Williams

Results of the qPCR experiments conducted with samples from infected Williams plants showed an increase in expression of most PL3-coding genes (except XM_009539509.1, whose expression did not appear to change) at all time points during infection when compared to mycelial expression. However, the degree to which transcript levels increased showed considerable variation. After the initial expression increase observed at 12 h post inoculation (hpi), there was a relative decrease in expression for three of the transcripts (between 12 and 24 hpi). By 36 hpi, all genes, with the exception of XM_009539513.1 and XM_009526324.1, showed a decrease in expression compared to the previous time points. By 48 hpi, when visible lesions on infected soybean leaves could be observed, there was a surge in expression of PL3-coding genes, with most (with the exception of XM_009523709.1) showing a considerable increase in transcriptional levels (Fig. 5).

Fig. 5

qRT-PCR analysis of P. sojae PL3-coding genes at several time points after inoculation of soybean plants (cv Williams) with P. sojae. Fold expression was calculated by the Livak (ΔΔCt) method using mycelia expression as the calibrating sample and actin as the reference gene. Bars represent standard error

Through 36 hpi, the transcripts that consistently experienced the highest expression were XM_009532098.1, XM_009526324.1, and XM_009526560.1, whose levels showed a several-hundred-fold increase compared to those of the mycelial stage (up to this point, the highest expression was detected in XM_009532098.1, which had an approximately 500-fold increase at 12 hpi). Moreover, most of the remaining transcripts experienced varying levels of up-regulation, ranging from a 10-fold to approximately a 400-fold increase compared to mycelial expression levels. Following a general trend, eight out of 11 expressed genes typically experienced a slight decrease in expression after 24 hpi, with the exception of transcript XM_009526560.1, whose expression continued increasing until peaking at 48 hpi (Fig. 5). All genes experienced further up-regulation at 48 hpi, with the exception of XM_009523709.1, whose expression was lower than at 36 hpi. Furthermore, four transcripts (XM_009522292.1, XM_009526560.1, XM_009526324.1, and XM_009532098.1) had the largest increase in expression (between 12,000- and 2200-fold) compared to mycelial levels.

Expression in cv. Williams 82

Results for Williams 82 showed a different expression pattern from that seen in Williams (Fig. 6). All genes increased their expression at 12 hpi and continued so at 24 hpi, with transcripts XM_009526401.1, XM_009523712.1, XM_009522292.1, XM_009526391.1, and XM_009539509.1 showing a greater than 2000-fold increase in expression compared to mycelial levels. In fact, transcripts XM_009526401.1 and XM_009523712.1 reached over a 6000-fold and 4000-fold increase in expression, respectively. By 36 hpi, five of the 11 genes showed increased transcriptional levels; in contrast, the expression in the remaining six decreased slightly. At 48 hpi the degree of expression of eight of the transcripts decreased considerably. However, two transcripts (XM_009526401.1 and XM_009526391.1), showed greater than 15,000-fold and 8000-fold increases in expression, respectively, while XM_009623712.1 expression remained approximately constant.

Fig. 6

qRT-PCR analysis of P. sojae PL3-coding genes at several time points after inoculation of soybean plants (cv Williams 82) with P. sojae. Fold expression was calculated by the Livak (ΔΔCt) method using mycelia expression as the calibrating sample and actin as the reference gene. Bars represent standard error

Discussion

CAZymes, enzymes involved in the metabolism of carbohydrates, have been presumed to play a large role in the initiation of infection by Phytophthora spp. due to their potential cell-wall degrading capabilities. In general, several of these enzymes, including cutinase, esterase, polygalacturonase, cellulase, and generalized lyases have previously been found to be expressed at higher levels throughout infection (Torto-Alalibo et al. 2007; Mingora et al. 2014; Horowitz and Ospina-Giraldo 2015). However, relatively little is known about the variation in gene expression during pathogenesis, or the specific role that many of these CAZyme gene families may serve. The main goal of this study was to characterize the sequences of the PL3-coding genes found in P. sojae and gather quantitative data concerning their relative expression levels during infection of soybean cultivars susceptible (Williams) and resistant (Williams 82) to P. sojae.

The large complement of CAZy homologs in the genome of P. sojae was first described comprehensively in 2010 (Ospina-Giraldo et al. 2010a). The catalog was based on the 1.1 version of the annotated genome released by the JGI-DOE and included at least 54 members of the polysaccharide lyase superfamily, 35 of which were classified as PL3-coding genes. Although new bioinformatic analyses (mostly, homology searches using hidden Markov models) conducted on the more recent 3.0 version of the P. sojae genome annotation have reduced the number of potential PL3-coding genes to 12 homologs, it is clear that P. sojae retains significant gene redundancy that may enable rapid and successful colonization of the plant host.

Phylogenetic analyses of the PL3-coding gene sequences revealed an active evolutionary history, exemplified by multiple gene duplications reminiscent of the plasticity found in the effector gene families (Tyler et al. 2006; Haas et al. 2009). Sequences located on the same scaffold are part of the same clade in the phylogenetic tree (Fig. 2) and several of them appear to have originated as a result of gene duplication events in their most recent ancestors. In fact, an estimate of the evolutionary divergence between sequences obtained using the Maximum likelihood model (Tamura et al. 2004) and the bootstrap values on all inferred trees support the hypothesis that at least two pairs of true paralogous genes exist in the P. sojae genome (Additional file 5: Table S3). Most of the PL3-coding genes are organized in groups, with their ORFs frequently positioned in opposite orientations; indeed, eight genes found on scaffold 8 are distributed in two clusters and the genes on each cluster are consistently separated by relatively short distances. The first of these clusters, which contains three genes, spans a sequence of approximately 23 kb while the second one, which contains the remaining five genes, spans a sequence of approximately 74 kb (Table 1). When the phylogenetic analysis was expanded to include sequences from five additional members of the genus Phytophthora, the evolutionary relatedness of PL3-coding genes within the genus was confirmed. Moreover, the distinctiveness of each P. sojae paralog was further supported (Fig. 3).

Although genome assembly and annotation algorithms have substantially improved, accuracy limitations still persist. As a result, some expressed sequences are incorrectly annotated as introns. There are limited experimental reports on P. sojae introns (Costanzo et al. 2006; Shen et al. 2011; Horowitz and Ospina-Giraldo 2015), which is partly explained by the relatively low density of 0.62 introns per gene (Shen et al. 2011). Currently, it is commonly accepted, based mostly on computational analyses, that the average intron size in P. sojae is 88 bp, with most introns between 60 and 90 bp in length (Shen et al. 2011). Although introns with sizes outside of that range have been computationally predicted (Ospina-Giraldo et al. 2010a), introns greater than 90 bp have rarely been experimentally confirmed. One exception is a 236 bp intron within a pectin methylesterase-coding gene reported by Horowitz et al. (Horowitz and Ospina-Giraldo 2015). Out of the 12 genes analyzed in our study, only one (XM_009539513.1) was confirmed to be interrupted by an intronic sequence. Furthermore, our PCR, RT-PCR, sequencing and bioinformatic analyses appear to validate our hypothesis that introns predicted in the gene models encoding transcripts XM_009526562.1, and XM_009526401.1 are likely an artifact of the gene prediction algorithm used at the time the P. sojae genome was annotated.

The very low percentage (9%) of intron-containing PL3-coding genes and the lack of multiple introns in one gene considerably diverge from what we previously reported for pectin methylesterase-coding genes in P. infestans (Mingora et al. 2014) but agree with the general consensus that P. sojae is an intron-poor organism (Shen et al. 2011). The gene density (0.09 introns per gene), in particular, is lower than the density values reported for fungi such as Aspergillus nidulans, Cryptococcus neoformans, Neurospora crassa, and Saccharomyces pombe, although slightly greater than that reported for S. cerevisiae (Kupfer et al. 2004). The reduced intron density we have found here and elsewhere (Horowitz and Ospina-Giraldo 2015) is in stark contrast with the hypothesis that Phytophthora has acquired a considerable number of introns in its early evolutionary history while experiencing only a modest decrease in intron density in comparison to its chromist ancestor (Csűrös et al. 2008). This discrepancy will only be resolved when additional evidence from genome-wide experimental studies can be obtained.

Glycine max cv. Williams is a standard soybean cultivar that has no Rsp (Resistance to Phytophthora sojae) genes within its genome (Gao et al. 2005) and, as such, is susceptible to P. sojae infection. However, cv. Williams is the recurrent parent to cv. Williams 82 (Haun et al. 2011), which is a resistant hybrid cultivar combined by crossing four resistant back-crossed lines of Williams and Kingwa (Bernard and Cremeens 1988). Phytophthora resistance is in large part due to the persistence of large Kingwa introgressions in the Williams 82 genome. These Kingwa introgressions contain several different Rsp genes, most notably Rsp1, which are actively transcribed during a P. sojae infection (Gao et al. 2005). Williams 82 is ideal for juxtaposition with Williams because the two cultivars are identical in all aspects except the Phytophthora resistance found in Williams 82 (Bernard and Cremeens 1988). Due to the presence of these Rsp genes in Williams 82, we expected P. sojae to display a differential gene expression profile when infecting susceptible and resistant soybean cultivars.

Although most evidence suggests that P. sojae (like P. infestans and many pathogenic fungi) produces a specialized hypha called the appressorium, which enables host penetration by the fungi that produce it, there is also evidence that, in Phytophthora spp., the force generated by this structure alone is insufficient to breach the rigidity and strength of the plant cell wall (Money et al. 2004). Therefore, it has been suggested that polysaccharide lyases may aid in the penetration process by playing a role in host cell-wall breakdown via polysaccharide degradation. However, there had previously been little evidence to support this hypothesis. In this study, the mRNA levels of all P. sojae PL3-coding genes were quantitatively determined in infected and control samples, providing a measure of their relative transcriptional activity.

Analysis of the qPCR data using T-tests revealed a statistically significant difference in gene expression at 12, 24, 36, and 48 hpi, in both cultivars, between the test (infecting) and calibrator (non-infecting) samples (the only exception being XM_009539509.1; see below). Furthermore, the qPCR results also revealed an up-regulation of most PL3-coding genes in Williams soybean by 48 hpi, a trend observed in other CAZy gene studies (Mingora et al. 2014; Horowitz and Ospina-Giraldo 2015), except for transcript XM_009539509.1, whose expression did not change compared to mycelial levels. However, the up-regulation of these genes was not linear, but was instead punctuated by relatively decreased transcriptional activity of most genes between 24 and 36 hpi. When compared to mycelium expression levels, 11 of the 12 genes showed heightened levels of expression at each time-point during the first 48 hpi but the patterns of expression change were distinct for each gene. Similarly, the quantified fold-difference values ranged widely, from just a few hundred to over 2000 times greater than corresponding mycelial expression levels.

Expression analysis of transcript XM_009539509.1 revealed an unexpected pattern. This transcript could not be detected in RT- and qPCR experiments using any of the Williams samples; however, it was highly expressed in Williams 82 (furthermore, successful amplification was obtained when P. sojae genomic DNA was used as template in PCR experiments). These results support the notion that, while XM_009539509.1 is not expressed by P. sojae during infection of Williams, it may play an unknown role in the interaction with a resistant host.

As in the experimental infection of Williams, the specific expression patterns of P. sojae PL3-coding genes in Williams 82 considerably varied for each gene, with transcriptional levels ranging from a minimal increase (XM_009539513.1 at time points 12 and 48 hpi) to more than 15,000-fold increase (XM_009526401.1 at 48 hpi) in comparison to mycelial levels. These high values are in agreement with those reported in other studies (Mingora et al. 2014; Horowitz and Ospina-Giraldo 2015). During the transition from 24 to 36 hpi, some genes showed an increased transcriptional activity while others did just the opposite. The most intriguing change occurred at 48 hpi, when only two transcripts (XM_009526401.1 and XM_009526391.1), showed increased expression (> 5000-fold and > 8000-fold, respectively). Interestingly, however, the degree of expression change of eight of the transcripts decreased considerably, while the expression of the remaining one (XM_009623712.1) remained approximately constant. This may perhaps indicate the end of an unsuccessful infection attempt.

Most importantly, in some cases, the transcript levels in Williams 82 were much higher than in Williams. For example, at 36 hpi, all genes exhibited considerably higher expression in Williams 82 compared to Williams. These results supported our hypothesis that, during interaction with Williams 82, P. sojae displays altered PL3-coding gene expression patterns in a (fruitless) effort to successfully colonize this resistant cultivar. These differences can be seen in terms of the specific transcripts that were up-regulated and in the degree of increase in transcriptional levels.

Taken together, these results suggest that P. sojae produces considerably lower levels of PL3-coding transcripts to successfully colonize susceptible Williams plants. Furthermore, in Williams, a specific subset of P. sojae transcripts (XM_009526324.1, XM_009532098.1, and XM_009526560.1, and at 48 hpi, XM_009522292.1) appears to play a leading role during infection. No similar statement (perhaps with the potential exceptions of XM_009526401.1 and XM_009539513.1, particularly at 48 hpi) can be made about the interaction with Williams 82 (Fig. 5 and Fig. 6).

Using the MEME-suite, three motifs within the 12 PL3-coding gene promoter regions were discovered, the largest one of which is shown (Fig. 4). This motif, in particular, was present in five of the 12 sequences (42%) and was located in a very similar location (between positions 22–105 upstream of the transcription start site) in all the sequences containing it. Noticeably, this motif bears a striking resemblance to the Eukaryotic Initiator and the Flanking Promoter Region (INR–FPR) expression-increasing promoter motifs, which may occur either separately, or combined (albeit with a two-bp spacer) to form a “supramotif” (YCAYTYY-NN-MWTTTNC), as described previously for Phytophthora genes (McLeod et al. 2004; Roy et al. 2013b). The motif found in the PL3-coding sequences fits within the predicted sequence of the “supramotif” and is only lacking the two-bp spacer that has previously been described (Roy et al. 2013b).

The genes that possess this motif within the promoter region also appear to be expressed in higher amounts in both cultivars of soybean compared to genes that do not contain the motif. This is in agreement with previous reports (Roy et al. 2013a; Roy et al. 2013b), which describe that the INR – FPR motif increased expression in 10%–15% of Phytophthora genes. Similarly, SignalP 4.1 analysis indicated that eight out of 12 predicted amino acid sequences possess a canonical signal peptide. As was the case with genes containing a common motif in their promoter sequences, these genes also appear to be expressed in higher amounts, on average, than those that do not contain a secretion signal, in both cultivars of soybean. These data are also in strong agreement with recent data obtained in our laboratory for the P. sojae CE Family 10 genes (unpublished).

All 12 predicted amino acid sequences also have the potential for non-canonical secretion, according to SecretomeP 2.0 (including the three sequences that were not predicted to contain a canonical signal by SignalP 4.1). Thus, all proteins encoded by the genes analyzed in this study have the potential to be secreted into the extracellular space, potentially facilitating direct targeting of the plant cell wall. Taken together, several highly expressed PL3-coding genes contain both the conserved motif in the promoter region and a signal-peptide coding sequence. These genes (XM_009532098.1, XM_009523712.1, and XM_009522992.1, in Williams, and XM_009523712.1 and XM_009522992.1, in Williams 82) were highly up-regulated during the infection process in comparison to other PL3-coding genes. Two of these genes, XM_009532098.1 and XM_009522992.1 appear within the same clade in our evolutionary analyses. Interestingly, transcripts XM_009526324.1 and XM_009526560.1, whose respective genes do not contain the conserved motif in their promoter region, were also up-regulated in Williams, although not as highly as the three previously mentioned. Other genes that were slightly up-regulated only possessed either the signal peptide or the motif, but not both. It is tantalizing, therefore, to hypothesize that possession of both the signal peptide and motif leads to a larger fold increase in expression during infection, whereas possession of only one of these elements may lead to a more moderate fold increase in expression. Furthermore, the presence of a sequence motif in the promoter of those genes with increased transcriptional activity may suggest the existence of a concerted regulatory mechanism critical to the infection process.

Indirect support for the potential role of CAZymes in P. sojae pathogenesis is observed when the activity data are compared to those of similar enzymes found in other organisms. For example, the fungus Uromyces fabae infects species of beans in a similar manner to that observed in oomycete pathogenesis. During interactions between U. fabae and its host, cell wall-degrading enzymes with PL activity firmly adhere to the host cell wall and participate in localized degradation via the alteration of carbohydrate structure and composition (Münch et al. 2008). After penetration of the host cell wall, hemibiotrophic organisms such as P. sojae, initially display a biotrophic activity. Shortly after infection, they transition into a necrotrophic stage, using the decaying plant tissue as an energy source. Differences in expression seen during this transition period may arise due to the dramatic shift in host-pathogen interactions as well as the changes in pathogen nutrition. In the fungal pathogen Colletotrichum graminicola, it is known that metabolism during biotrophy depends entirely on obtaining nutrients from the living host cell (Münch et al. 2008). However, during the necrotic phase of growth, C. graminicola resumes active infection processes, deploying secondary hyphae. These hyphae use a vast array of carbohydrate-active enzymes and depolymerases to cause extensive degradation of the host cell wall (Münch et al. 2008).

The dramatic up-regulation and nearly synchronous peak in expression of PL3-coding genes in P. sojae during the first 24 hpi is similar to that reported by previous studies documenting the patterns of pathogenesis-related gene expression over the course of infection. In fact, the first 24 h following inoculation, when initiation of infection and biotrophy occur, has been continuously noted as the most dramatic period of widespread gene up-regulation (Fig. 5 and Fig. 6) (Moy et al. 2004; Horowitz and Ospina-Giraldo 2015; Hinkel and Ospina-Giraldo 2017). The up-regulation (compared to mycelial levels) observed at 12 and 24 hpi is still observed at 36 hpi for each of the eleven genes studied. Furthermore, in Williams and Williams 82, four and two of the analyzed genes, respectively, experienced their highest levels of expression at 48 hpi. Interestingly, two of the genes consistently experiencing the highest levels of expression in Williams, XM_009522292.1, XM_009532098.1, are grouped within a single clade. However, this type of close evolutionary relationship is not observed for the genes that experience the highest upregulation in Williams 82, which were placed in very different clades by our evolutionary analyses.

The resulting parabolic nature observed in the majority of these expression curves is further supported by previous findings (Moy et al. 2004) and by what is currently known about the generalized oomycete life cycle. Following initiation of infection, P. sojae exhibits a rapid progression from biotrophy to necrotrophy, a change that characterizes the life cycle of a hemibiotrophic pathogen, and which generally occurs after 24 h of active infection (Moy et al. 2004). This dramatic physiological transition fundamentally alters the nature of host-pathogen interactions as well as many internal characteristics of the pathogen itself, such as growth and metabolism (Lee and Rose 2010). Depending on the specific environmental conditions and organismal strains, associations between the host and pathogen may still be at various stages of development at 24 hpi, whereas by 48 hpi, the infection has confidently progressed well into necrotrophy (Moy et al. 2004). The expression profile observed for the PL3-coding genes is very similar to that previously reported for P. sojae genes encoding pectin methylesterases (PMEs), enzymes belonging to the Carbohydrate Esterase family 8 (CE8), which also experience a rapid increase of expression within the first 36 hpi (Horowitz and Ospina-Giraldo 2015). Other previous studies show that in P. infestans PMEs are also clearly upregulated within the first 24 hpi (Mingora et al. 2014). However, less of a clear temporal association was observed in that case, and it appeared that the expression of several P. infestans PMEs varied over a large range during the first 48 hpi, while still remaining up-regulated (Mingora et al. 2014). Thus, it is likely that PL3-coding genes are important to infection and are working in conjunction with other virulence proteins to establish an infection.

Further research is needed to advance our understanding of the role PL and other CAZyme-encoding genes may have in Phytophthora spp. and other important agricultural pathogens. Biochemical analysis and physical modeling would provide key insights into the structure and true function of PL3 enzymes, further elucidating their role in growth and pathogenesis. In addition, new insights into this process will increase the feasibility of developing effective control mechanisms in the future.

Conclusions

We have conducted an extensive sequence analysis of all putative PL3-coding gene models found in the revised version of the annotated P. sojae genome. In addition, we have quantified the relative expression of each gene in P. sojae during infection of susceptible (Williams) and resistant (Williams 82) soybean cultivars over a 48 h time period. Results indicate that during infection of the Williams cultivar, several of the PL3-coding genes experience significant up-regulation during the first 48 h of infection. An increase in transcriptional activity is also observed during infection of Williams 82; however, this expression pattern is distinct from that observed in the Williams cultivar. These results provide evidence regarding the potential involvement of PL-coding genes in the early pathogenic processes of P. sojae and suggest that this pathogen expresses genes differentially, depending on whether infection is occurring on a resistant or a susceptible soybean cultivar.

Methods

Living material

Cultures of Phytophthora sojae strain R4 were obtained from the United States Department of Agriculture (USDA) in Beltsville, Maryland. The isolate was subcultured and maintained on Rye A agar plates for long-term laboratory storage at 18 °C.

Glycine max seeds from cultivars Williams and Williams 82 were planted in terracotta pots containing Miracle Grow® potting soil. Plants were watered daily and kept at 24 °C in a growth chamber with a 16 h light/8 h dark photoperiod (Arriagada et al. 2004).

Growth media

Rye A agar was prepared by first soaking 60 g rye seeds in distilled water for approximately 48 h until germination had occurred. Grain particles were then strained out of the liquid and placed in a beaker with enough distilled water to completely cover the grain. This mixture was homogenized using a blender on high setting for 2 min. Following homogenization, the mixture was incubated for 1 h at 68 °C and then passed through 4 layers of gauze to remove any solid particles remaining in the mixture. The filtrate was then combined with the original germination supernatant along with 20 g sucrose and 15 g Difco® granulated agar. The final volume was brought to 1 L with distilled water and the resulting solution was autoclaved for 20 min at 15 psi before being poured into plates (Caten and Jinks 1968).

Unclarified V8/Lima bean agar media was prepared by first autoclaving 80 g of frozen baby Lima beans in 200 mL distilled water for 10 min at 15 psi. The solid Lima beans were then filtered out of the mixture after it had been allowed to cool sufficiently. Next, 1.4 g of calcium carbonate and 200 mL of V8 juice were added to the liquid. The volume was then brought to 1 L with distilled water. The resultant pH was measured and adjusted to a value of 6.0 by adding the necessary amount of KOH. Following the addition of 30 g Difco agar, the mixture was autoclaved for 15 min at 15 psi and poured into plates upon cooling (Miller 1955).

Infection assays

Infection assays were conducted as previously described (Schmitthenner and Bhat 1994; Ryley et al. 1998; Dorrance et al. 2008). Cultures of P. sojae were grown in unclarified V8/Lima bean agar plates at 25 °C in the dark for two weeks, until the mycelium covered the plates. Mycelial inoculum was prepared by slicing the V8/Lima bean agar into long strips, which were then placed in a sterile syringe. Subsequently, the culture was macerated by forcing the agar through the syringe, and collected in a sterile 50 mL conical tube. At this time, the macerated culture was placed back into the syringe, and a sterile needle was attached to the syringe. Three-week-old plants were infected by producing a 2 cm lesion into the soybean hypocotyl using the sterile needle and injecting between 0.2 and 0.4 mL of macerated P. sojae culture into the lesion. Plants were covered with a clear plastic bag for the first 12 hpi to maintain high humidity. Infected plants were then incubated at 25 °C until tissue collection (obtained from the first true leaves), which occurred at 12, 24, 36, and 48 hpi. Samples were then flash-frozen in liquid nitrogen and stored at − 80 °C. In addition to these inoculated samples, a control sample was obtained in the same manner from each plant just before inoculation.

Sequence analysis

Database mining

From our database of predicted P. sojae CAZy homologs, which had been identified using the 1.1 version of the annotated P. sojae genome assembly released by the Department of Energy Joint Genome Institute (JGI-DOE), gene models listed as members of PL3 were further investigated. These original gene models (a number of which were incomplete or wrongly predicted in the genome assembly) were used for BLAST-based comparisons against the currently available version (3.0) of the annotated P. sojae genome database and a new set of predicted gene models, more accurate than the one based on version 1.1, was generated. The database for automated carbohydrate-active enzyme annotation (dbCAN), which provides automated CAZyme annotation for any genome or set of nucleotide/protein sequences using improved Hidden Markov Models (HMM), was used to confirm that the selected P. sojae CAZy gene models did, in fact, belong to the PL3-coding gene complement. Furthermore, in order to validate the dbCAN results and confirm the presence of a pectate lyase PFAM domain, PL3 protein sequences were analyzed using the HMMER tool (Schuster-Böckler et al. 2004) at the European Molecular Biology Laboratory (EMBL-EBI. https://www.ebi.ac.uk/Tools/hmmer/).

Genes to be included in subsequent studies were required to meet several criteria including an E-value (confidence measure) less than 1 × 10^− 25, inclusion in the current JGI-DOE complete P. sojae genome annotation, and presence of no more than three introns, each one no longer than 90 bp (because experimental evidence has shown that Phytophthora genes usually have no more than three introns, each not exceeding 90 bp in length, the number of introns – and their respective size – within a predicted transcript is often indicative of the likelihood that the proposed annotation is accurate).

Promoter and signal peptide analysis

A 500 bp sequence upstream of the specified transcription start site of each PL3-coding gene was downloaded from GenBank. Sequences were scanned for the presence of common motifs using the Multiple Em for Motif Elicitation Suite (MEME), version 5.0.2 (http://meme-suite.org/tools/meme). Protein sequences were analyzed for the presence of canonical and non-canonical signal peptides using SignalP 4.1 and SecretomeP 2.0, respectively, found at http://www.cbs.dtu.dk/services/SignalP/ and at http://www.cbs.dtu.dk/services/SecretomeP/.

Evolutionary analysis

PL3-coding gene sequences were aligned using Clustal W (Thompson et al. 1994), available within the MEGA (Molecular Evolution Genetic Analysis) suite, version 7.0.26 (Kumar et al. 2016). Data provided by the Clustal W alignment were used to construct phylogenetic trees based on DNA sequence similarity. Three separate trees were constructed using the Neighbor Joining, Maximum Parsimony, and Maximum Likelihood algorithms. Robustness of the inferred trees was tested by bootstrap analysis consisting of 1000 replications (Felsenstein 1985). While the Neighbor-joining method uses estimations of phylogenetic diversity to construct hypothesized evolutionary relationships (Saitou and Nei 1987), the Maximum Parsimony one prioritizes minimizing the total number of character state changes, thereby minimizing the total amount of homoplasy (Fitch 1977). Maximum Likelihood, in turn, evaluates the statistical likelihood of a proposed evolutionary history, and then creates a tree with the highest probability (Felsenstein 1981).

Primer design

For all valid PL3-coding gene models, specific primers were designed using the NCBI Primer-BLAST algorithm (https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_LOC=BlastHome) following specifications provided by the Bio-Rad iQ5 Multicolor Real-Time PCR handbook, which ensured that the designed primers were appropriate for use in reverse transcription (RT) and quantitative (q) PCR experiments. Primer sequences can be found in Additional file 4: Table S2.

RNA extraction

Phytophthora. sojae RNA was extracted from the leaves of infected soybean plants. To obtain samples from a variety of infection stages, primary leaf samples were collected at 12, 24, 36, and 48 hpi. As a control, an RNA sample from non-infecting mycelium grown in pea broth (Hinkel and Ospina-Giraldo 2017) was also collected. Extraction of RNA from these samples was completed using the Qiagen RNeasy Mini RNA extraction kit, following the suggested procedure for fungi. This procedure included the optional DNase digestion to ensure elimination of genomic DNA. Successful completion of RNA isolation was then confirmed with agarose gel electrophoresis and spectrophotometry.

cDNA synthesis, RT-PCR, and qPCR

An RT-PCR reaction was performed with Qiagen’s Quanti-Tect Reverse Transcription Kit according to the suggested protocol. Extracted P. sojae RNA (mycelial and infection time-point samples) and the designed primers were used to amplify specific regions of the desired PL3-coding mRNA molecules, thus creating a PL3 cDNA library. qPCR was conducted using Qiagen Quantitect SYBR Green Supermix (Qiagen, Germantown, MD). Each cDNA sample produced in the previous RT-PCR reaction was analyzed in triplicate. Actin, which is constitutively expressed in P. sojae, served as a reference gene while non-infecting mycelial P. sojae RNA was used as a calibrator. Analysis of qPCR results used the ΔΔC_t (also known as Livak) method, which provides a quantitative value for relative differences in expression level of each gene (Livak and Schmittgen 2001). Quantified expression levels were then normalized to the expression of the reference gene for reliable multi- sample comparison and analyzed according to the expression levels detected in the mycelial calibrator. Statistical significance for qRT-PCR results was assessed by a paired Student’s t-test between non-infecting mycelial and infected sample gene expression for each of the PL3-coding genes studied.