Skip to main content

Interaction between Rag genes results in a unique synergistic transcriptional response that enhances soybean resistance to soybean aphids



Pyramiding different resistance genes into one plant genotype confers enhanced resistance at the phenotypic level, but the molecular mechanisms underlying this effect are not well-understood. In soybean, aphid resistance is conferred by Rag genes. We compared the transcriptional response of four soybean genotypes to aphid feeding to assess how the combination of Rag genes enhanced the soybean resistance to aphid infestation.


A strong synergistic interaction between Rag1 and Rag2, defined as genes differentially expressed only in the pyramid genotype, was identified. This synergistic effect in the Rag1/2 phenotype was very evident early (6 h after infestation) and involved unique biological processes. However, the response of susceptible and resistant genotypes had a large overlap 12 h after aphid infestation. Transcription factor (TF) analyses identified a network of interacting TF that potentially integrates signaling from Rag1 and Rag2 to produce the unique Rag1/2 response. Pyramiding resulted in rapid induction of phytochemicals production and deposition of lignin to strengthen the secondary cell wall, while repressing photosynthesis. We also identified Glyma.07G063700 as a novel, strong candidate for the Rag1 gene.


The synergistic interaction between Rag1 and Rag2 in the Rag1/2 genotype can explain its enhanced resistance phenotype. Understanding molecular mechanisms that support enhanced resistance in pyramid genotypes could facilitate more directed approaches for crop improvement.


Plants are in a constant battle against pathogens and pests. A common outcome of the selective pressure imposed on plants by pests is the evolution of plant resistance genes (R) that can recognize pest challenges [1]. The detection triggers a physiological response that reduces the invader’s ability to colonize or feed on the plant. In turn, the selective pressure thus imposed on the pest can result in the evolution of different strains or biotypes with the ability to evade recognition by specific R gene products. This well-defined evolutionary arms race has been identified as one of the processes that shape plant and pest evolution [2]. Many types of R genes have been identified [1], although the most abundant and well-characterized belong to the intracellular nucleotide-binding/leucine-rich-repeat (NLR) class [3]. These receptors detect pathogen or herbivore effectors and evoke a defense response, effector-triggered immunity (ETI) that shares common signaling events and outputs with the basal defense response triggered by recognition of pathogen or herbivore associated molecular patterns (PTI); although the ETI response is faster and seems to overcome negative regulators that normally temper the strength of the PTI response [3].

Plant response to stress needs to display both robustness and plasticity simultaneously [4]. In the case of plant immunity, robustness is necessary to mount a strong defense in the presence of a large variety of pests, while plasticity is needed to display a different flavor of the defense arsenal as specific pathogens change [3, 5]. Activation of ETI can be triggered by direct interaction between an effector and a NLR, indirectly through NLRs that can sense host targets modified by effectors, or by sensing effector-decoy interactions [6]. It has been proposed that, upon activation, sensor NLRs transmit the signal to helper NLRs that do not interact directly with effectors [5, 6]. Moreover, while a large number of sensor NLRs exist in many angiosperm species, helper NLRs are fewer; thus R gene signaling converge on helper NLR nodes to create a regulatory network with built-in redundancy [5, 7]. These regulatory networks also converge downstream in a series of phytohormone crosstalk events, mainly synergistic or compensatory interactions between the jasmonate (JA) and salicylate (SA) pathways, that contribute to the connectivity of the network [8, 9]. Finally, these signals trigger a transcriptional reprogramming that is used by the plant to mount different chemical, structural, and physiological changes depending on the nature of the attacking pest [3]. Signal integration that leads to a specific transcriptional response is controlled by transcription factor (TF) regulatory networks that rely on a selected group of TFs and interacting proteins to provide specificity [10, 11]. This network architecture with highly connected nodes, including helper NLRs, phytohormone regulatory factors, and TF hubs, imparts plasticity and robustness to the system.

R genes are commonly used as part of management strategies to control pests in crop production, as a sustainable alternative to chemical treatments. However, the selective pressure exerted on the pathogens or insects favors the evolution of pests, negatively impacting the durability of individual R traits [12, 13]. One strategy used by plant breeders to overcome this limitation is the use of gene pyramids: the incorporation of multiple different R genes able to control a specific pathogen or insect in an individual crop variety [14]. Gene pyramids have been successfully used to manage different bacterial, fungal, and insect pests in a variety of species including rice, wheat, and soybean, among others [14, 15]. In many cases, these crop varieties with multigenic resistance have enhanced pest control when compared with monogenic lines, indicating an additive interaction between R genes [14,15,16,17]. In other cases, however, this interaction may result in incompatibility and less favorable outcomes [18]. These R gene interactions suggest that each R gene triggers a subset of the defense arsenal. The combination of multiple R genes results in complementary, or antagonistic, outcomes that can produce plants with a more robust or compromised defense response depending on the nature of this interaction [15]. While these interactions are well-described from a phenotypic perspective and pyramids are used frequently in breeding programs, the molecular events underlying the observed traits are normally not known.

The soybean (Glycine max)-soybean aphid (Aphis glycines) system is an attractive model to study the effects of gene pyramiding. The soybean aphid, an invasive insect species, is one of the most economically damaging pests of soybeans in the United States, causing up to 40% yield reduction if left unmanaged [19, 20]. Resistance to soybean aphids is provided by Resistance to Aphis glycines (Rag) genes. Twelve Rag genes and four aphid resistance quantitative trait loci have been identified to date [21], although none have been cloned so far. The resistance provided by Rag1 and Rag2 has been well-characterized, and transcriptome studies describing the response of plants carrying either gene to soybean aphids are available [22,23,24]. The Rag1 gene was the first resistance gene to be identified in the soybean cultivars Dowling and Jackson [25,26,27]. The gene is located in a 115 kb interval on chromosome 7 [28]. Rag2 was identified in two plant introductions, PI 243540 and PI 200538 [29, 30], and it is located in a 54 kb interval on chromosome 13 [31]. Both these regions contain NLR genes which have been proposed as candidates for the resistance genes [28, 31]. Both the Rag1 and Rag2 genes have been shown to primarily confer an antibiosis type of aphid resistance [32]. The Rag1 transcriptional response to aphid feeding involved repression of transcripts related to cell wall modifications, induction of SA regulated plant defense, induction of the phenylpropanoid pathway, and upregulation of cuticle production [22, 23]. The Rag2 transcriptional response to aphid feeding involved upregulation of genes involved in cell wall modifications, secondary metabolism, hormone metabolism and stress signaling while transcripts involved in carbon metabolism and photosynthesis were downregulated [24].

While different Rag genes are available for breeding programs, the discovery of soybean aphid biotypes [33,34,35] that can successfully colonize resistant soybean genotypes led to the development of soybean lines carrying gene pyramids with two or three resistance genes [17, 36,37,38]. Greenhouse and field phenotypic studies showed that pyramid lines carrying Rag1 and Rag2 can control aphid populations more effectively than monogenic resistance lines [17, 36, 39], and no adverse effects have been reported for this interaction. Given that each Rag gene seems to trigger a somewhat different transcriptional response, we hypothesized that the enhanced resistance observed at the phenotype level would be based on an additive effect at the transcriptional level, that is, the transcriptional response of the Rag1/2 pyramid should overlap with the transcriptional responses of Rag1 plus Rag2. Alternatively, the interaction of signaling events triggered by Rag1 and Rag2 in the pyramid could lead to an output different from the addition of individual Rag gene responses. In this case, the interaction can be considered synergistic and is measured at the molecular level by the number of genes whose expression changes significantly only in the pyramid. This type of transcriptional synergistic response has been previously observed in the cellular response to drug combinations when compared to the response to monotherapies [40]. The enhanced transcriptional response, due to either additive or synergistic interaction, would then translate into a more robust resistant phenotype.

To test these hypotheses, we conducted a transcriptome analysis using RNA sequencing (RNA-seq) to determine the molecular response of four soybean genotypes to soybean aphids 6 or 12 h after infestation. The soybean genotypes used were: IA3027RA12 which contains both the Rag1 and Rag2 genes, IA3027RA1 which has Rag1 alone, IA3027RA2 which has Rag2 alone and IA3027, the susceptible control. These soybean genotypes are referred to as Rag1/2, Rag1, Rag2 and susceptible, respectively. Our results showed that, unexpectedly, the Rag1/2 response presents a unique transcriptome signature that indicates a synergistic interaction between the Rag1 and Rag2 genes.


Confirmation of enhanced resistance in the pyramid genotype

A no-choice experiment was used to determine whether the resistance phenotype of the monogenic and pyramid genotypes could be observed using our experimental conditions. As expected, the Rag1/2 soybean genotype had significantly reduced aphid populations compared to genotypes with the Rag1 or Rag2 gene alone, both of which had significantly fewer aphids than the susceptible control (Fig. 1A). There were no significant differences in aphid populations between soybean genotypes that contained the Rag1 or Rag2 gene alone.

Fig. 1
figure 1

A. aphid resistance phenotypes of the genotypes used in this study. Phenotyping results for the four soybean genotypes obtained using no-choice experiments conducted in growth chambers. Data was normalized to the susceptible control, IA3027. Statistical data analysis for each comparison was done using Student’s t-tests (Two-Sample Assuming Unequal Variances). Error bars represent the standard error. Different letters indicate statistically significant differences (p < 0.05). B. Global changes in the soybean transcriptional response to aphids. Total number of DE genes for aphid versus mock comparisons for each soybean genotype at each time point is shown as induced genes (positive fold change) and repressed genes (negative fold change). Total number of induced and repressed genes for each comparison are indicated next to each of the bars

Transcriptional response to soybean aphids

The transcriptional response to aphids of each genotype was assessed to understand the molecular effect of pyramiding aphid resistance genes in soybean. The response was analyzed at 6 and 12 h after infestation because previous experiments have shown that the most significant transcriptome changes in resistant plants are observed earlier than 24 h after aphid infestation [22,23,24, 41].

In total, 1,042,480,636 reads were generated from 48 leaf samples, with an overall mapping rate of > 95% with reference to version 2 of the Williams 82 soybean genome [42]. Of those, 4646 differentially expressed (DE) genes responded to soybean aphids in one or more of the four soybean genotypes for the eight (aphid versus mock) comparisons made for both time points (Additional File 1). The number of DE genes for each of the three aphid-resistant soybean genotypes was larger than in the susceptible response 6 h after infestation. In contrast, at the 12 h time point, the number of DE genes was similar among the four soybean genotypes. The majority of the DE genes (> 70%) were repressed for all four soybean genotypes (Fig. 1B) at the 12 h time point.

We also examined the expression differences between soybean genotypes in the absence of aphids by comparing mock samples across genotypes and time points (Additional File 2). Comparing the three resistant genotypes to the susceptible line, and the Rag1/2 pyramid resistance line to the single resistance gene genotypes, at six and twelve hours, identified between 75 and 218 DE genes per comparison. Among these genes, of specific interest were those located within the Rag1 [28] and Rag2 loci [24]. In version 2 of the soybean genome, the Rag1 locus corresponds to Gm07: 5,531,331 to 5,769,789 and contains the range of genes from Glyma.07G062300 to Glyma.07G064700 (Fig. 2). Six genes from this region were DE between genotypes when comparing Rag1 to susceptible, Rag1/2 to susceptible or Rag1/2 to Rag2 (Glyma.07G062400, Glyma.07G063100, Glyma.07G063300, Glyma.07G063600, Glyma.07G063700 and Glyma.07G064400, Table 1). Four of the DE genes in the region have BLASTP homology to soybean NLRs (E = 0, Glyma.07G063100, Glyma.07G063300, Glyma.07G063600 and Glyma.07G063700). Of these, Glyma.07G063100 was expressed at somewhat higher level in Rag1-containing genotypes in three comparisons (Table 1), while Glyma.07G063700 was expressed at much higher level in all Rag1-containing genotypes across all six comparisons (FC between 5.6–10.4; Table 1). This suggests that the orthologs of Glyma.07G063700, or less likely Glyma.07G063100, correspond to Rag1. These genes were not identified as candidates for Rag1 by Kim et al. [28], who used the first release of the soybean genome (G. max 1.0) in their study. Glyma07g06910 (corresponding to Glyma.07G063100) was predicted as a low confidence gene model with no NLR homology in the first genome release. The G.max 1.0 gene corresponding to Glyma.07G063700 (Glyma7G07145) was located outside of the genomic interval corresponding to Rag1. However, between the G. max 1.0 and G. max 2.0 versions of the soybean genome, additional sequencing and markers were used to fill gaps and improve sequence assembly. The genomic interval between the markers used to map Rag1 in genome version 2 has increased by 123 kb and now includes Glyma.07G063700, making it a high priority candidate for Rag1 (Fig. 2).

Fig. 2
figure 2

Comparison of the Rag1 locus in versions 1 and 2 of the Williams 82 genome. The top line represents the sequence of the Rag1 locus in version 1 of the soybean genome as reported by Kim et al. (2010). The bottom line represents the Rag1 locus in version 2 of the Williams 82 genome assembly. Conserved regions are represented by red shading. Candidate NLRs are represented as arrows and color-coded as follows: not differentially expressed (white), expressed more in the susceptible genotype (grey) or expressed more in the resistant genotypes (green)

Table 1 Differentially expressed genes among genotypes, in the absence of aphid feeding near the Rag1 and Rag2 loci

The same approach was used to identify candidate genes for Rag2. In G. max 2.0, the Rag2 locus corresponds to Gm13: 30,412,470 to 30,466,615 and contains the range of genes from Glyma.13G190400 to Glyma.13G191200. Three genes from this region were DE between genotypes when comparing Rag2 to susceptible, Rag1/2 to susceptible or Rag1/2 to Rag1 (Glyma.13G190400, Glyma.13G190800 and Glyma.13G190900, Table 1). Of these three, two have BLASTP homology to soybean NLRs (E = 0, Glyma.13G190400 and Glyma.13G190800), however, both were repressed across all six genotype comparisons. Unexpectedly, both genes were repressed in response to aphids in Rag1 at 12 h (Additional File 1) suggesting the Rag1 and Rag2 loci could somehow interact. Glyma.13G190400 (Glyma13g25970 in G. max 1.0) was proposed as a candidate for Rag2 gene based on gene expression [24], while Glyma.13G190800 (Glyma13g26000 in G. max 1.0) has been proposed as the candidate Rag2 gene based on fine mapping [31]. Repression of these genes in our experiment could suggest they are not candidates for Rag2, but it is important to note that the time points used in this study are just a small snapshot of the resistance response. It has previously been proposed that the presence of Rag1 and Rag2 genes causes constitutive expression of some defense genes even in the absence of herbivores [23, 43]; however, this effect seems to be minor for the genotypes analyzed here.

Transcriptional effect of pyramiding Rag1 and Rag2 genes at the global level

We hypothesized that the Rag1/2 pyramid genotype response to aphids would be the results of an additive effect of the monogenic genotypes, i.e. the DE transcriptome of the Rag1/2 plants would be the sum of DE genes in the Rag1 response and the Rag2 response. In this case, we expected a nearly complete overlap such that there would be few DE genes unique to the Rag1/2 transcriptome response to soybean aphids. Alternatively, a synergistic response would be indicated by a set of genes only DE in the pyramid and absent from either individual Rag genotype response. When the lists of DE genes for each genotype were compared, there were unique and common gene sets at each time point (Fig. 3A). Rag1/2 showed large unique DE sets at both time points (1000 genes at 6 h and 259 genes at 12 h), indicating that pyramiding the Rag1 and Rag2 genes in one soybean genotype results in a synergistic effect on aphid resistance at the transcriptome level. The DE genes unique to the Rag1/2 response are referred to as “synergistic genes.” The Rag1/2 synergistic response to soybean aphids comprised 65% of the DE genes at 6 h and 19% of the DE genes at 12 h. There was a fivefold increase in the number of DE genes that were common to all four soybean genotypes at 12 h when compared to the 6 h time point. The number of other unique and common DE genes varied among respective comparisons (Fig. 3).

Fig. 3
figure 3

The pyramid genotype has a unique transcriptome response. A. DE genes in each genotype after aphid infestation for 6 h (left) and 12 h (right) unique to each genotype or shared among different genotypes. The DE genes unique to the Rag1/2 response were designated as the “synergistic genes” set. B. Clustering heatmap showing expression patterns for the 4646 genes differentially expressed after aphid infestation in any genotype after 6 h or 12 h. Each column represents the aphid versus mock comparison for each soybean genotype and time point combination. Each row represents a gene. Genes included in this figure are differentially expressed for at least one comparison. Z-scores were calculated for each gene in all comparisons. Soybean genotypes are designated as: Sus: Susceptible, R1: Rag1, R2: Rag2, and R1R2: Rag1/2. C. Heatmap corresponding to the 1000 genes differentially expressed exclusively in the pyramid genotype at 6 h after aphid infestation, referred in the text as the synergistic response

In total, there were 4646 aphid-responsive genes in all aphid versus mock comparisons at both time points. Cluster analysis was conducted on these genes to visualize how aphids impacted gene expression on a global level. The DE genes modified by aphids clustered by time point (Fig. 3B). Interestingly, at 6 h, the Rag1/2 gene expression pattern was more similar to that of Rag2, while at the 12 h time point the Rag1/2 response was closer to the Rag1 response. Based on similarity in gene expression pattern, DE genes grouped into six clusters (Fig. 3B). Expression patterns for DE genes on the clustering heatmap showed a characteristic set of genes that were induced by aphids at 6 h in all soybean genotypes but repressed at 12 h (clusters 1 and 4). Other sets of genes were mostly repressed at 6 h but induced 12 h after aphid infestation (clusters 2 and 5, Fig. 3B). Subsets located in clusters 2 and 3 with distinctive expression patterns in the 6 h Rag1/2 response contained some of the synergistic genes. The distinctive region in cluster 2 had DE genes that were strongly repressed in the 6 h Rag1/2 response only but not modified by aphids in the other three soybean genotypes. In cluster 3, the distinctive region contained DE genes that were strongly induced in the Rag1/2 response only. The rest of the synergistic genes were distributed among the clusters with differences in expression patterns.

Given the size of the early synergistic gene set, this 6 h set (1000 DE genes, Additional File 1) was reanalyzed to evaluate the possibility that these genes were the result of a quantitative additive effect, that is, genes that were marginally regulated in the monogenic genotypes and only became DE in the pyramid as a result of additive effects that placed their expression above our arbitrary FDR cutoff. A cluster analysis of the synergistic genes for the 6 h time point (Fig. 3C) identified two clusters of similar size that were defined by the Rag1/2 genotype response (either repressed or induced). Of the 1000 Rag1/2 DE genes depicted, 664 (66%) were significant in Rag1/2 (FDR < 0.05), but not significant in Rag1 or in Rag2 using a relaxed cutoff (FDR > 0.25). A similar approach has been used to differentiate genes affected by combined stresses versus the effect of individual stress conditions [44]. We also compared the observed log2FC values for the 6 h synergistic gene dataset to the expected additive values (log2FCRag1 + log2FCRag2) (Additional File 3, Fig. S1) in aggregate using a χ2 test. This analysis showed that the observed values for the full 6 h synergistic dataset are significantly different from the expected additive values (χ2 = 2.96− 28). Analysis of the induced or repressed gene subsets, however, showed that observed values for repressed synergistic genes are not significantly different from expected additive changes (χ2 = 1), while observed values for induced synergistic genes are significantly different from expected additive values (χ2 = 4.08− 92). These results indicate that the synergistic gene set is, at least in part, uniquely regulated by aphids when the two resistance genes are present in the same genotype.

Dynamics of the response to aphids for each genotype

To determine the time progression of the response of each soybean genotype to aphids, gene lists (aphid versus mock) for each soybean genotype at 6 h and 12 h were compared (Additional File 3, Fig. S2). Overall, there was very little overlap of DE genes at both time points (< 6%) within each soybean genotype. The number of DE genes that were common at both time points for each soybean genotype was: 38 for the susceptible, 108 for Rag1, 114 for Rag2, and 108 for Rag1/2. Within each soybean genotype, some DE genes changed in the same direction at both time points, but others showed opposite regulation (induced at 6 h but repressed at the 12 h time point or vice versa). The very small overlap of DE genes between the two time points for each soybean genotype indicates a dynamic response to aphids for all genotypes. Moreover, given the limited overlap in the transcriptional responses at 6 h versus elevated similarity at 12 h, the results suggest that earlier events determine the differences in the resistance phenotype observed for each genotype.

Importantly, most 6 h synergistic genes (910 genes out of the 1000 gene set) are not DE in the 12 h response in Rag1/2 or any other genotype, suggesting that the 6 h synergistic response in the Rag1/2 line is not the result of a more rapid induction of the defenses triggered by individual Rag genes.

Biological processes modified by aphids in each soybean genotype

Based on the unique global response of the Rag1/2 genotype and its enhanced resistance compared with each of the monogenic genotypes, we hypothesized two scenarios that could be occurring in Rag1/2 plants. First, it is possible that the synergistic gene set is simply a larger number of genes that function within the same or similar biological pathways as those regulated in either of the monogenic responses. In this case, we would expect to find a large overlap in the type of response (biological processes) in Rag1/2 when compared to each monogenic genotype. Alternatively, each response, and particularly the synergistic response observed for Rag1/2, could regulate unique biological processes that would result in distinct metabolic outcomes.

To determine the similarities and differences in the biological processes occurring within each resistance response, an analysis of gene ontology (biological processes) terms that were significantly overrepresented at corrected P < 0.05 was conducted for each genotype. Many biological processes were significantly overrepresented in only one genotype, and few were shared at 6 h of aphid infestation (Fig. 4A; Additional File 4). The Rag1/2 response showed a number of overrepresented biological processes related to plant immunity (including phytohormone (JA and SA) mediated signaling, SA biosynthesis, secondary cell wall biogenesis, regulation of plant-type hypersensitive response, regulation of hydrogen peroxide metabolism, defense response - incompatible interaction, systemic acquired resistance, detection of biotic stimulus, MAPK cascade and regulation of multi-organism process, and others). The Rag1 response uniquely involved biological processes related to chloroplasts (i.e., chloroplast organization, chloroplast RNA processing and protein targeting to chloroplast, and transcription from plastid promoter) and several primary metabolism processes. Chalcone biosynthesis and response to gravity were both uniquely overrepresented in the Rag2 response. For each of the three aphid-resistant soybean genotypes, the majority of the defense-related transcripts were induced (Additional File 4). In the susceptible response, no defense-related GO term was significantly overrepresented at 6 h. However, the unique susceptible response involved starch and maltose biosynthesis with the DE genes being mostly induced.

Fig. 4
figure 4

Distinct biological processes are regulated early in the pyramid response to aphids. Comparison of GO terms for biological processes significantly modified after 6 h (A) or 12 h (B) of infestation by soybean aphids in the susceptible (blue bars), Rag1 (red bars), Rag2 (green bars) and Rag1/2 (yellow bars) soybean genotypes. On the Y axis, soybean genotypes for which specific GO terms (biological process) were significantly overrepresented (p ≤ 0.05) are indicated in parenthesis. S: Susceptible, R1: Rag1, R2: Rag2, and R1R2: Rag1/2

Contrary to the findings for the 6 h response, biological processes modified by aphids after 12 h were mostly similar among the four soybean genotypes (Fig. 4B). Eleven out of 24 overrepresented GO terms at the 12 h time point were common to all four soybean genotypes. Of these, five were likely related to plant defense responses, including calcium ion transport, regulation of hormone levels, response to hormone stimulus, cell wall organization and plant-type cell wall biogenesis (Fig. 4B). For most of the common GO terms, the DE genes were repressed in all soybean genotypes (Additional File 4).

Overall, our results showed that the transcriptional responses of the four soybean genotypes at 6 h were mostly unique to each genotype and affected distinct biological processes. Specifically, GO analysis indicated that in the 6 h Rag1/2 response, the biological functions that were differentially regulated by aphid presence were not significantly overrepresented in the response of the individual Rag1 or Rag2 genotypes. On the other hand, at 12 h the transcriptional responses and the biological processes were more similar among the four soybean genotypes.

Biological processes modified by soybean aphids in the Rag1/2 synergistic response

To examine the biological processes regulated by the synergistic genes in the Rag1/2 soybean genotype, GO analysis was conducted for DE genes that were unique to this genotype at each time point. There were 1000 transcripts unique to Rag1/2, and 13 biological process terms were significantly overrepresented (Table 2) 6 h after aphid infestation. No significantly overrepresented GO terms were identified in the 12 h response. A number of significantly overrepresented GO terms for the 6 h synergistic response were related to plant defense (secondary cell wall biogenesis, regulation of multi-organism process, detection of biotic stimulus, systemic acquired resistance, SA mediated signaling, and MAPK cascade), and most of the DE genes in these categories were induced by aphids (Table 2). Another group of overrepresented terms were involved in photosynthesis-related biological processes (photosynthesis-light harvesting, light reaction, photosystem II assembly, and chlorophyll biosynthesis), but the majority of these transcripts were repressed in the Rag1/2 response (Table 2). The majority of the synergistic transcripts included in the defense-related GO categories were distributed in clusters 1, 2 and a few were in cluster 3 in Fig. 3B, and in cluster 2 in Fig. 3C. Most of the photosynthesis-related transcripts, strongly repressed only in Rag1/2, were located in cluster 2 in Fig. 3B, and in cluster 1 in Fig. 3C. This further confirmed the uniqueness of the Rag1/2 response to aphids when compared to the other three soybean genotypes.

Table 2 Biological processes differentially represented in the Rag1/2 synergistic response to soybean aphids after 6 h

Hormone signaling in the plant response to soybean aphids

Given the strong representation of phytohormone signaling in our gene ontology results and that the crosstalk among hormone signals is an important component of the R gene-mediated response, we identified phytohormone signatures using the Hormonometer tool [45] that compares a given dataset to several Arabidopsis datasets of hormone-regulated gene expression. Since the tool has been developed for Arabidopsis, it was adapted by first identifying the closest Arabidopsis homolog for each DE gene, as has been previously described [46, 47]. Hormonometer results (Additional File 3, Fig. S3), which determine the correlation of a given dataset to a hormone response dataset, indicated that genes associated with the two main defense hormones, JA (MJ, Methyl Jasmonate, for the Arabidopsis datasets) and SA, are upregulated in response to aphids in all four genotypes. The JA response seems to be higher in the resistant genotypes, more specifically in Rag1/2 and Rag2 at 6 h, and it tapers down at 12 h. This result agrees with a recent experiment showing accumulation of JA-Ile in Rag1/2 plants after 6 h of aphid infestation [48]. The SA response seems more sustained through the two time points in all genotypes, although weaker overall in the susceptible line. In contrast, another hormone that is frequently related to defense responses, ethylene (ET), is generally suppressed throughout the period analyzed in all genotypes, with a stronger suppression in the susceptible line. Auxin signaling was observed at 6 h but was less prevalent at 12 h. ABA response was observed at 12 h for all genotypes. Both gibberellin and brassinosteroid signaling were generally repressed, while cytokinin signaling was not prevalent in any genotype. Although it provided an important insight into phytohormone activity during the response to aphids, this analysis failed to identify a phytohormone signaling pattern that could explain the large difference in gene expression observed for the Rag1/2 line at 6 h after aphid infestation.

Identification of transcription factors that participate in the Rag1/2 synergistic response to aphids

Another important component of a plastic defense response is the presence of transcriptional networks that add regulatory nodes to signaling pathways. We used the SoyDB transcription factor database [49] to identify DE TF in order to determine their participation in the Rag1/2 synergistic response to soybean aphids. Given the importance of JAZ transcriptional repressors in control of the JA pathway [50,51,52], JAZ homologs were also added to our analysis even though these proteins do not bind DNA directly. Like the full transcriptome set, TF expression patterns (Fig. 5) exhibited distinct responses among genotypes at 6 h but similar responses at 12 h. Across all treatments, there were 522 TF that responded to infestation by soybean aphids, belonging to 44 TF families.

Fig. 5
figure 5

The pyramid response includes a large number of transcription factors (TF). Expression patterns of TF that were differentially expressed between aphid-treated and mock-treated samples for each soybean genotype at both time points were plotted (with and absolute fold change larger than 1.5). The x axis indicates absolute fold change in gene expression, and the y axis shows TF families. Several differentially expressed transcription factors per family are plotted for each comparison

We identified DE TFs belonging to 33 different families in the 6 h response, with 107, 89, 149 and 43 differentially expressed TF in Rag1, Rag2, Rag1/2 and susceptible genotypes, respectively. We identified individual DE TFs that were strongly differentially expressed in response to aphids and TF families where multiple members respond to aphid feeding (Fig. 5). Ten TF families of interest include: ZIM, WRKY, TPR, NAC, MYB/HD-like, JUMONJI, Homeodomain/HOMEOBOX, C2C2 (Zn) CO-like, bHLH and AP2-EREBP. For most of these families, the magnitude of gene expression is greatest in Rag1/2. While many of the TF families had mixed expression patterns, genes in the WRKY and ZIM families were highly induced in the three resistant genotypes. TFs belonging to 39 families were present in the 12 h response. We identified 192, 169, 201 and 152 differentially expressed TF in Rag1, Rag2, Rag1/2 and susceptible genotypes, respectively. Ten TF families responded to aphids at 12 h: WRKY, TPR, NAC, MYB/HD-like, C2H2 (Zn), Homeodomain/HOMEOBOX, C2C2 (Zn) CO-like, BZIP, bHLH and AP2-EREBP (Fig. 5). Unlike the 6 h timepoint, by 12 h after aphid treatment, most of these families respond similarly across all four genotypes.

The difference in Rag1, Rag2 and Rag1/2 TF expression patterns at 6 h again suggests a synergistic response. To evaluate this further, we identified all TF families significantly over represented within each genotype and time point and within the 6 h and 12 h synergy genes. C2C2(Zn)CO-like TFs were significantly over represented at 6 h in Rag1 (Corrected P < 3.63E-12), WRKY TFs were significantly overrepresented in Rag2 at 6 h (Corrected P < 4.52E-2), and both C2C2(Zn)CO-like (Corrected P < 2.72E-4) and WRKY (Corrected P < 5.72E-7) TFs were significantly overrepresented in Rag1/2 at 6H. The Rag1/2 genotype had the highest number of DE WRKY TF (23 WRKYs) when compared with Rag2 (11 DE WRKYs) or Rag1 (4 DE WRKYs), while in the susceptible control only 1 DE WRKY transcription factor was observed, and this WRKY was repressed (Additional File 5). Remarkably, 13 WRKY TFs (57% of Rag1/2 TFs) were differentially expressed only in Rag1/2 and are included in the 1000 synergy genes identified at 6 h. WRKY TFs were significantly overrepresented among the 91 TFs found within the 1000 synergy genes (Corrected P < 2.53E-3). In contrast, when we examined the C2C2(Zn)CO-like TFs, only two of the 11 (18%) differentially expressed genes in Rag1/2 appeared synergistic. At 12 h, BZIP TFs were significantly over represented in Rag1 (Corrected P < 1.53E-2), while TCP TFs were significantly overrepresented in Rag1, Rag2 and Rag1/2 (Corrected P < 3.30E-4, 6.97E-3, and 4.06E-3, respectively). Of the 10 TCP TFs expressed in Rag1/2 at 12 h, three corresponded to 12 h synergistic genes.

With these examples in mind, we then reexamined the data to identify all TF families at 6 and 12 h where 50% or more of the TFs expressed in Rag1/2 also corresponded to synergistic genes. This approach identified 19 TF families at 6H (Additional File 6), including WRYKY (57%), BHLH (75%), AP2-EREBP (71%) and NAC (86%) TFs, the major plant TF families involved in regulation of plant defense [11]. Similarly, we identified ten TF families at 12 h with both NAC and WRKY TF familied identified at both 6 and 12 h. Remarkably, this approach also identified ZIM transcription factors as important. While the soybean genome contains 24 ZIM TFs, two (orthologs of AtJAZ1 and AtJAZ2) were unique to Rag1/2 at 6 h (Fig. 5). JASMONATE-ZIM-DOMAIN (JAZ) TF act as repressors of the JA pathway [53]. Overall, the observed expression patterns for the WRKY and JAZ TF suggests that the Rag1/2 response is distinct from that of the monogenic genotypes and that for a subset of genes the magnitude of change in expression is larger in the pyramid than in either of the single Rag genotypes.

The large number of TF differentially regulated in the Rag1/2 response associates well with the observed transcriptional changes in this genotype in response to aphids. However, on its own it does not explain why the pyramid response is so different from the response observed in the Rag1 and Rag2 monogenic genotypes. One possibility to explain this observation is that when Rag1 and Rag2 are combined, the suite of TF regulated by Rag1 signaling interact and cooperate with the repertoire of TFs in the Rag2-mediated response and combine in unique ways to induce changes in the expression of the synergistic genes. Thus, we could expect to see TF that are shared between Rag1/2 and each of the individual Rag genotypes to act as nodes in a TF network to regulate the synergistic response. To test this idea, we performed a STRING analysis to identify potential interactions among TF differentially expressed in the Rag1/2 response. Since there is scarcity of soybean-specific experiments incorporated in available databases, we performed this analysis using the Arabidopsis homolog more closely related to each DE soybean TF. After removing TF that did not participate in any interaction, each remaining TF was cataloged to indicate whether it was DE only in the Rag1/2 response or if it was also DE in other genotypes. Given the ancient hexaploidy nature of the soybean genome [42], several soybean genes can have close homology to a unique Arabidopsis gene. For this annotated network (Fig. 6), if more than one soybean TF corresponding to the same Arabidopsis homolog was present in the Rag1/2 DE dataset, it is indicated as a split circle. Differences in regulation are indicated with different colors. The analysis identified a large Rag1/2-specific TF network with high connectivity and several smaller networks with few genes. Remarkably, of the 82 TFs depicted in Figs. 6, 52 correspond to synergistic TFs.

Fig. 6
figure 6

Transcription factor networks in the 6 h pyramid response are regulated by Rag1 and Rag2. Protein-protein interaction network for TF differentially regulated in the Rag1/2 genotype after 6 h of aphid infestation. The network was built using the closest Arabidopsis homolog for each soybean TF. Colors indicate whether each TF was DE in the pyramid only or in the pyramid and other genotypes in the 6 h response. Split circles indicate more than one soybean homolog corresponding to the same Arabidopsis TF

The large network showed sectors with different levels of influence from individual Rag genes (Fig. 6). The expected cooperative contribution from Rag1 and Rag2 signaling, represented with blue (genes DE in Rag1 and Rag1/2) and red (genes DE in Rag2 and Rag1/2), is apparent in sector 1, centered around MYB15. Two soybean TFs homologous to the AtMYB15 gene were differentially upregulated, one in Rag1 and Rag1/2 (Glyma.02G005600, FC = 2 fold and FC = 2.8 respectively) and the other in Rag2 and Rag1/2 (Glyma.20G209700, FC = 3.6 fold and FC = 5.2 respectively). These TF are connected to other TF DE in Rag2 and the pyramid (homologs of AtATAF1, AtWRKY33 and AtWRKY40) or DE in all resistant genotypes (AtWRKY46), and this cluster of TF is then connected to a series of synergistic (yellow) TF. Sector 1 contains TF mainly associated with phytohormone signaling, defense and response to stress, and genes in this sector are mainly induced in the response to aphids. Sector 2 has less connected nodes, but also shows TF that are regulated by Rag1 or Rag2 signaling interacting with each other. Particularly interesting are the homologs of AFO (two homologs, Glyma.05G056000 DE in Rag1 and Rag1/2, and Glyma.17G138200 only in Rag1/2) and PHB (two homologs, Glyma.15G129700 DE in Rag2 and Rag1/2, and Glyma.09G023600 only in Rag1/2) that interact with each other and several synergistic TF. The TF in this sector are mainly associated with development of axial symmetry, epidermal cell patterning, and chlorophyll biosynthesis, and expression patterns suggest that these processes are repressed in the pyramid genotype. Sector 3 contains TF involved in circadian regulation and flowering, with many of the genes showing DE in response to aphids in all genotypes, including susceptible plants.


Gene pyramiding has become an important plant breeding and resistance management tool to develop more durable resistant crop cultivars against insect pests and pathogens. Thus, the effect of gene pyramiding has been widely investigated at the phenotypic level. For instance, resistant cultivars containing multiple R genes have been developed to effectively protect plants against soybean aphids [17, 36,37,38], Soybean mosaic virus [54], bacterial blight and blast in rice [55,56,57], late blight of potato [58, 59], wheat leaf rust [60], and wheat powdery mildew [61]. In many cases, the combination of R genes in an individual genotype results in enhanced resistance when compared to genotypes with one R gene (for a summary see [14]). However, few reports have studied the molecular bases for the increase in resistance observed in R pyramids.

Recently, Kamphuis et al. [16] analyzed the effects of pyramiding AKR and AIN, two R genes that confer resistance to the bluegreen aphid (BGA) in Medicago. Pyramiding these genes in a single genotype results in enhanced resistance to BGA, while at the same time the pyramid eliminates deleterious effects associated with AIN. The analysis identified an important crosstalk among JA and SA signaling and regulation of reactive oxygen species (ROS) and cell death associated with R responses. In response to BGA attack, plants carrying only AIN deploy a strong hypersensitive response with cell death mediated by SA and ROS. On the other hand, the response in plants carrying AKR shows a unique JA-mediated response (not observed in basal defense or in the AIN genotypes). The pyramid genotype responded to aphids with a significant increase in JA and SA signaling that seem to act cooperatively to induce a stronger defense response against BGA. However, it is possible that JA also acts to reduce the hypersensitive response driven by AIN, thus eliminating negative effects associated with widespread cell death [16]. In a different study, Gao et al. [62] analyzed the effect of pyramiding Xa21 and xa5, two R genes that provide resistance to bacterial blight in rice. This study used a genome-wide approach to identify transcriptional changes associated with the presence of individual R genes or a pyramid in a single genotype; however, the analysis was conducted in the absence of the pathogen. Phenotypic analysis indicated that the pyramid had higher level of resistance to bacterial blight than either of the single R gene lines. Transcriptome data identified 2367, 2412, and 3596 DE genes in the NILs carrying Xa21, xa5, and Xa21 + xa5 respectively, when compared against the susceptible control line. Since no pathogen was present in these experiments, it is clear that in the case of the rice R genes there is a constitutive defense component associated with resistance. In addition, although there was a large number of genes (1633) DE only in the pyramid line, the focus of the study centered on the DE genes common to all R lines [62]. Thus, a deeper mechanistic understanding of the molecular bases for enhanced resistance in the pyramid was missed.

Pyramiding the Rag1 and Rag2 genes confers enhanced resistance to soybean aphids. Previous work indicated that, in plants carrying either of the two genes individually, a differential transcriptional response occurs early in the interaction with the insect [22,23,24]. We observed that the transcriptome response of the monogenic Rag1 and Rag2 genotypes were dissimilar to each other at 6 h after aphid infestation. Moreover, the Rag1/2 pyramid had its own strong, distinct response to aphids at this early time point. Interestingly, at 12 h, the R gene-mediated response of the monogenic genotypes and the pyramid were more similar. This suggests that events occurring during the first few hours after the initiation of aphid infestation are responsible for the enhanced resistance seen in the pyramid line. While some additive effects were encountered, the differential response observed in the pyramid is not purely the sum of the transcriptional responses observed in each monogenic genotype, indicating a synergistic interaction between Rag1 and Rag2 signaling. The cluster dendrogram showed that the Rag1/2 transcriptional response to aphids was more similar to that of the Rag2 genotype at 6 h but closer to the Rag1 response at 12 h, suggesting that modification of transcripts associated with the Rag2 resistance likely occurs earlier compared to transcripts involved in the Rag1 response in the pyramid genotype. This result suggests some hierarchical order in the Rag1/2 response, with Rag2 having more influence on the early response, although the majority of the DE genes observed at 6 h in the pyramid genotype are exclusively regulated in this genotype.

The Rag1/2 pyramid does not seem to follow the patterns observed in the R pyramid studies described above. We did not observe unique phytohormone signatures in the pyramid, and transcriptome differences between the pyramid and monogenic genotypes in the absence of aphids were minor. However, we observed a large transcriptional reprograming unique to the pyramid genotype. This synergistic response, defined as a set of genes DE only in the pyramid and not in either the Rag1 or Rag2 genotypes, is rapidly mounted, mostly evident at 6 h after the beginning of aphid infestation, and it affects a series of unique biological processes that may increase the defensive output in this genotype resulting in the enhanced resistance phenotype observed. As an example, our TF analyses suggest that, upon activation, Rag1 and Rag2 induce expression of different TF that cooperate on specific promoters to trigger this unique response. Many of the targets are also TF, creating an amplification of this synergistic effect.

Additionally, many TF act as homo or heterodimer complexes. Aphid-dependent induction of several highly similar soybean TF, close homologs to the same Arabidopsis TF, would provide greater complex diversity, creating novel target specificities. A similar scenario has been described for Arabidopsis allotetraploids created by crossing A. thaliana and A. arenosa autotetraploids [63]. In the allotetraploid progeny, differential accumulation of specific homoeologous TF (AtWRKY18, AaWRKY40, and AtWRKY60) was observed in response to bacteria and SA treatments. Changes in protein-protein interactions resulted in the formation of different heterodimers and these different complexes showed altered affinity for target promoters. These changes cause altered regulatory networks that could explain the enhanced resistance observed in the allotetraploids [63]. In the TF regulatory network identified in the Rag1/2 pyramid, two hubs with concerted action of “parental” TF, that is, TF DE in one of the monogenic genotypes and the pyramid, were predicted. One sector is centered on the interactions among homologs of MYB15, WRKY33, WRKY40, and WRKY46 (Sector 1, Fig. 6); the other is centered on the AFO-PHB hub (Sector 2, Fig. 6). In both cases, each node represents two or more soybean homologs of a unique Arabidopsis gene, suggesting the potential for a diversity of TF complexes. These hubs may in turn control the phenotypic outputs that increase the pyramid resistance against aphids, such as the biological processes uniquely regulated in the early (6 h) pyramid response that include defense responses related to increase in secondary cell wall deposition and suppression of photosynthesis.

While the TF network presented here is based on homology to Arabidopsis genes, and thus has limited power due to the uncertainty of assigning orthologs and potential functional divergency, it allows as to make predictions that can direct future functional studies.

The MYB-WRKY hub is likely to participate in the establishment of the early defense responses of the pyramid genotype. Cooperation of MYB and WRKY TF on individual promoters has been shown for the activation of defense against herbivores in soybean [64], and the roles of Arabidopsis MYB15, WRKY40, WRKY46, and WRKY33 in defense responses, including defense against aphids, are well-documented [65,66,67,68,69,70]. In fact, these four TF were among a set of TF specifically induced by cabbage aphid feeding in an experiment comparing Arabidopsis responses to insect and bacterial attacks [71], AtMYB15 is necessary for full harpin-induced resistance to green peach aphid in this plant [72], and silencing a wheat homolog of AtWRKY33 results in reduced resistance to the Russian wheat aphid [73]. Five soybean homologs of AtWRKY40 are induced in the early pyramid response to aphids in our dataset. Two of those (Glyma.17G222300 and Glyma.14G103100) are also induced in the Rag2 response, while the other three (Glyma.06G061900, Glyma.15G003300, and Glyma.07G023300) are only expressed in the synergistic set. We had previously observed increase of Glyma.17G222300 expression as part of the soybean PTI response to long-term aphid feeding [41]. In addition, Glyma.17G222300 and Glyma.14G103100 are also among the most highly induced genes in the resistance response to soybean cyst nematode (SCN) in wild soybean [74], and overexpression of Glyma.14G103100 resulted in decreased susceptibility to soybean cyst nematode in the susceptible soybean cultivar Williams 82 [75]. Two homologs of AtMYB15 are DE in the pyramid 6 h set, Glyma.20G209700 (also induced in Rag2) and Glyma.02G005600 (also induced in Rag1). Glyma.20G209700 is also upregulated in the PTI response to long term aphid feeding [41]. Glyma.02G005600 is also induced in the resistant response to SCN [74], suggesting that the MYB-WRKY hub also participates in the resistance response to other herbivores.

Cell wall modifications are a significant component of plant-insect interactions [76,77,78], and our analysis indicated that secondary cell wall biosynthesis is differentially regulated in the pyramid 6 h synergistic response. MYB TF are important regulators of the phenylpropanoid pathway and lignin biosynthesis [69]. AtMYB15 acts as a positive regulator of defense-induced lignification in Arabidopsis by binding to promoters of genes involved in G-lignin biosynthesis, and plants overexpressing AtMYB15 had increased lignin content [79]. Similarly, overexpression of CmMYB15 resulted in an increase in lignin accumulation in chrysanthemum and increased resistance to the aphid Macrosiphoniella sanborni [80]. Another important MYB TF in the regulation of secondary cell wall biosynthesis is AtMYB83 [81, 82], and two homologs of AtMYB83 (Glyma.11G133700 and Glyma.12G057900) are induced in the 6 h synergistic response and are part of the MYB-WRKY hub. Interestingly, homologs of several direct targets of AtMYB15 and AtMYB83 are differentially regulated in the pyramid at 6 or 12 h after aphid infestation, including homologs of PAL1, C4H, and COMT (MYB15 targets [79]), and MYB43, PAL1, MYB32, Family 8 and Family 43 glycosyltransferases, and EXPL1 (MYB83 targets [82]). Heterologous expression of Medicago truncatula WRKY108715 (AtWRKY33 homolog, and closest homolog to the synergistic gene Glyma.09G280200 in the MYB-WRKY hub) caused a large increase in lignin content in transgenic tobacco plants [83]. Thus, regulatory interactions between MYB and WRKY TF could be a significant contribution to cell wall strengthening which could limit the ability of aphids to penetrate leaf tissues to reach the phloem and in this way increase resistance in the pyramid genotype.

As mentioned, MYB TF are important regulators of the phenylpropanoid pathway which is also responsible for the biosynthesis of isoflavonoids in soybean, and quantitative trait locus analyses showed a correlation between resistance to soybean aphid and isoflavone content [84]. In a previous experiment, we showed that isoflavones can act as deterrents for soybean aphids [41]. Thus, another interesting outcome of MYB regulation is the production of isoflavones. One of the MYB15 homologs in MYB-WRKY hub, Glyma.20 g209700 (induced in Rag2 and Rag1/2) was identified in a genome-wide association study as one of the loci significantly associated with isoflavone concentration (called GmMYB29 in that study) and overexpression and RNAi of these gene resulted in increased and decreased isoflavone content respectively when compared to controls, in experiments using transgenic soybean hairy roots [85]. The other MYB15 homolog, Glyma.02G005600 (induced in Rag1 and Rag1/2) has also been implicated in the regulation of isoflavone biosynthesis. This TF, called GmMYB29A2 in the study, specifically controls glyceollin I biosynthesis in response to the oomycete Phytophthora sojae [86]. Moreover, several transcripts corresponding to the isoflavonoid branch of the phenylpropanoid pathway are also induced in the pyramid at 6 or 12 h, including isoflavone reductase (Glyma.11G070600), isoflavone 2′-hydroxylase (Glyma.11G051800), 2-hydroxyisoflavanone dehydratase (Glyma.10G250300), and isoflavone 4′-O-methyltransferase (Glyma.13G173300 and Glyma.13G173600). Thus, it appears that both the Rag1 and Rag2 pathways contribute to isoflavone production in the pyramid genotype and these phytoalexins could also be part of the defense output contributing to the increased-resistance phenotype observed for this genotype.

Suppression of photosynthesis and chlorophyll biosynthesis are also differentially regulated by aphid infestation in the synergistic response at 6 h. AtWRKY40 acts as a negative regulator of expression of nuclear genes that encode chloroplast proteins [87]. Specific AtWRKY40 targets include genes corresponding to light-harvesting chlorophyll a/b-binding (LHCB) proteins and a glutamyl-tRNA reductase involved in chlorophyll synthesis (HEMA1) [87, 88]. In addition to the WRKY40 homologs identified in the pyramid 6 h response, several photosynthesis-related WRKY40 targets are significantly repressed in the synergistic response, including homologs of LHCB1.4 (Glyma.16G165200, Glyma.16G165500, Glyma.16G165800, Glyma.16G162600), LHCB3 (Glyma.02G080800), LHCB6 (Glyma.05G119000), and HEMA1 (Glyma.06G091600, Glyma.14G185700, Glyma.02G218300). Thus, the MYB-WRKY hub identified in the 6 h pyramid response could also regulate the repression of photosynthesis observed in the early pyramid response. However, hubs in other sectors of the TF network identified in the 6 h pyramid response may also participate in the negative regulation of photosynthesis. A homolog of GNC (GATA, nitrate-inducible, carbon metabolism-involved) is repressed in the pyramid [Glyma.13G103900 (repressed in Rag2 and Rag1/2), Sector 2 in Fig. 6]. GNC controls expression of chloroplast genes and chlorophyll production in Arabidopsis, and gnc mutants have reduced levels of chlorophyll [89, 90]. Specific targets of GNC [91], such as LHCB1.4, are among the genes repressed in the 6 h synergistic set. Phytochrome A signal transduction (PAT1), a positive regulator of LHCB and HEMA1 expression [92] is also part of the sector 2 network (Glyma.12G216100, synergistic gene) and is repressed at 6 h. In addition, TF in sector 3 that include light-regulated and circadian clock-responsive factors could also participate in the photosynthesis regulation observed in the pyramid response [93].

Repression of photosynthesis is commonly associated with ETI responses [94,95,96,97]. Current models indicate that these changes in photosynthesis result in increases in reactive oxygen species (ROS) that are produced due to uncoupling of the electron transfer chain. The oxidative burst associated with these ROS is a key component of the ETI response, contributing to the hypersensitive response and other signaling functions, cell wall strengthening, and direct effects on pathogens [98]. Activation of the MPK3/MPK6 signaling cascade is critical for photosynthesis inhibition, ROS production, and ETI [97]. While the role of chloroplasts has not been well-characterized, ROS participation in defense responses in plant-aphid interactions have been previously reported [16, 22, 99,100,101,102]. In our analysis, a strong suppression of photosynthesis is observed in the early pyramid synergistic response. In addition, “regulation of plant-type hypersensitive response”, “regulation of hydrogen peroxide metabolism”, and “MAPK cascade” are among the biological processes overrepresented only in the Rag1/2 genotype dataset at 6 h (Fig. 4), and MPK3 (Glyma.U021800) is induced only in the Rag1/2 line at this time point. Thus, the suppression of photosynthesis observed in the Rag1/2 genotype may be associated with the production of ROS that could contribute to a hypersensitive-type response and also to cell wall-strengthening.

The differential Rag1/2 response appears transient, that is, a large number of genes DE at 6 h in this genotype return to basal levels by 12 h. A similar rapid and transient response to aphids was observed for Rag5-containing soybean plants when compared to the response of a susceptible NIL [43]. In that case, a large number of genes is DE in the infested Rag5 line at 6 or 12 h after aphid infestation when compared with the susceptible line, but only a few genes are DE between lines at 48 h post-infestation, indicating that the resistance response is rapid and may be established between 6 and 12 h after aphid infestation, as previously suggested [22]. Similarly, the susceptible and resistant response of peach plants to green peach aphid infestation show low level of similarity (23%) at 3 h post-infestation but a much larger overlap (59%) at 12 h [103]. It is also important to note that while gene expression may have subsided by 12 h, the proteins encoded by genes DE in the early pyramid response could have long half-lives, and many of the changes triggered by the pyramid response, such as strengthened cell walls or accumulation of phytoalexins, could persist for a longer period of time.


Pyramiding resistance genes is an effective method for crop improvement, but the molecular mechanisms underlying the enhanced resistance phenotype are not well-understood. We found that individual aphid resistance genes regulate different subsets of the defense repertoire in soybean, and that combining two R genes that confer resistance to soybean aphids in the same genotype results in an unexpected synergistic effect. The R gene-mediated response of monogenic genotypes and the pyramid are mostly different early (6 h after initial aphid infestation) but become more similar to each other and the susceptible response after 12 h, suggesting that the mechanisms that contribute to resistance are established very rapidly in this system. The unique, large synergistic response observed in the pyramid is likely responsible for the improved resistance phenotype expressed in this genotype, and it seems to be the result of a transcriptional network that combines Rag1- and Rag2-dependent TF in central interaction hubs that control target genes, including other TF, differentially regulated only in Rag1/2 plants. Our results suggest that the pyramid genotype elicits a defense response that is the combination of additive and synergistic effects. These results could explain some of the genetic variation observed in different soybean aphid biotypes. In order to overcome the unique response deployed by Rag1/2, aphids would need to evolve virulence factors that can overcome not only the individual Rag1 and Rag2 mechanisms, but also the Rag1/2 synergistic response. This is in part corroborated by the recent finding that the genome structure (SNPs) of soybean aphid biotype 4 (virulent on Rag1, Rag2, and Rag1/2) is not simply the additive effect of biotype 2 and biotype 3 (virulent on Rag1 and Rag2, respectively), but rather the majority of the sequence variation in its genome was unique to biotype 4 [104]. This may suggest that alternative mechanisms are needed for overcoming pyramided resistance. The observation of additive and synergistic responses seen in Rag1/2 soybeans is highly beneficial from a resistance management perspective, as soybean aphids need to overcome three defense mechanisms (Rag1, Rag2, and the synergistic response), which could greatly extend the durability of these varieties. This synergistic interaction of R gene signals causes photosynthesis perturbances that produce ROS and contribute to a rapid hypersensitive-type response, production of feeding-deterrent compounds, and strengthening of secondary cell walls that build up rapid barriers to aphid feeding and establishment of successful colonization. Understanding the specific subsets of defense responses activated by individual R genes and the way they interact when combined in the same plant genotype can lead to more effective approaches to biotechnological crop improvement and plant breeding as well as preservation of resistance traits.


Plant and insect material

Four soybean genotypes developed at Iowa State University were used for this study. The soybean genotypes were: IA3027RA12, the Rag1/2 soybean genotype that has both the Rag1 and Rag2 aphid resistance genes; IA3027RA1 contains Rag1 alone; IA3027RA2 has Rag2 alone; and IA3027, the susceptible control that has no aphid resistance gene. These soybean genotypes were developed using a backcrossing scheme in which the Rag1 and Rag2 donors were A08–123074 and LD08-89051a respectively, with the aphid-susceptible genotype, IA3027 used as the recurrent parent. Soybean genotypes IA3027, IA3027RA1, and IA3027RA12 are approximately 93.75% genetically identical [105, 106]. The genotype with Rag2 is an experimental line sharing 75% of it genes with the recurrent parent IA3027 [17]. For all experiments conducted in this study, soybean aphid biotype 1 was used, which is unable to colonize soybean genotypes with any known Rag gene. Aphids were obtained from a laboratory-maintained colony. Prior to use, aphids were raised on IA3027, the susceptible control. The growth chamber for the aphid colony was maintained at a constant temperature of 25 °C and photoperiod of 16 h light, 8 h darkness.

Experimental design, treatments and sample collection

Two main sets of experiments were conducted for this study: a) Phenotyping experiments to confirm expected aphid populations for each soybean genotype, and b) Transcriptome response to aphid feeding. All experiments were conducted in growth chambers maintained at a constant temperature of 25 °C and photoperiod of 16 h light, 8 h darkness. To eliminate seed or soil-borne pathogens, seed were surface sterilized with chlorine gas for 16 h as described by Paz et al. [107], and soil used for all experiments was steam-sterilized Metro Mix 900 growing mix (Sun Gro Horticulture, Vancouver, BC, Canada).

Phenotyping experiments

To confirm the aphid resistance/susceptibility of each soybean genotype, two no-choice growth chamber experiments arranged in a completely randomized design, with three replicates (pots) for each genotype were conducted. Seeds were individually grown in 5.5″ by 6″ plastic pots. Plants were fertilized once a week with a 1:1 mixture of all-purpose Scott’s Miracle-Gro Excel (21–5-20, The Scott’s Company LLC, Marysville, Ohio, USA) and Cal-Mag Miracle-Gro Professional (15–5-15, The Scott’s Co.), applied at a rate of 12.5 mL/L water. When plants reached the V3 growth stage [108], 30 mixed-age apterous aphids were placed on the abaxial side of the V3 middle leaflet of each plant using a fine paint brush and clip cages were placed on the infested leaflet to restrict aphid movement. One week later, aphid populations on each plant were determined. The phenotypic response of the different soybean genotypes was compared with Student’s t-test.

RNA-seq experiment and tissue collection

After confirmation of the expected aphid phenotypes, soybean plants of all four genotypes were grown in five growth chambers using a completely randomized design. All chambers were set to similar light, temperature and relative humidity conditions (Light: ~ 350 umol/m2/s, Temperature: 25 °C, Relative humidity: ~ 50%), and planting was done as described previously. Mock plants were grown in separate growth chambers from aphid-treated plants to avoid the potential effects of priming [23]. All plants were staked using bamboo stakes to keep them upright. When plants reached the V3 growth stage, clip cages were placed on the V3 middle leaflet of mock plants in each growth chamber, but without aphids. Thirty mixed-age apterous soybean aphids were placed on the abaxial side of the V3 middle leaflet of aphid-treated plants in each growth chamber and clip cages placed on the leaflet to restrict aphid movement. V3 middle leaflet samples (mock or aphid-treated) were individually collected from each plant after 6 h, gently removing aphids from each leaflet, and flash-frozen in liquid nitrogen. The collected leaflets of mock plants were gently brushed as well to mimic removal of aphids. Twelve hours after introduction of aphids to plants, leaf samples were collected from a separate set of plants. Three plants were used for each sample (genotype x treatment x time), and three replicates were used for the experiment. In total, 144 samples were collected for this study. Prior to processing, all samples were stored at − 80 °C.

RNA sample preparation and RNA-seq

All flash frozen leaf samples were ground using a mortar and pestle in liquid nitrogen. RNA was extracted from individual samples using the Qiagen® RNeasy® Plant Mini Kit (Qiagen®, Germantown, MD) with some modifications to the manufacturer’s protocol. Samples were incubated at 56 °C for two minutes and vortexed occasionally to disrupt tissue. Additionally, three rounds of RPE buffer washes were done instead of two. RNA samples were DNase treated using the Ambion® TURBO DNA-free kitTM (Ambion®, Austin, TX) to remove genomic DNA contamination. DNase-treated RNA samples were cleaned with Qiagen® RNeasy® MiniElute Cleanup Kit (Qiagen®, Germantown, MD). To ensure the integrity of RNA samples, all samples were checked using an Agilent® 2100 BioanalyzerTM (Agilent®, Santa Clara, CA) and samples with RNA Integrity Number (RIN) of ≥6.9 considered to be of good quality for RNA-seq. RNA from three leaflets (obtained from three different plants) was pooled in equal amounts for each RNA pool to be analyzed. A total of 48 RNA pools (3 replicates × 2 time points × 2 treatments × 4 soybean genotypes) were submitted to the DNA facility at Iowa State University for multiplex library preparation and single end sequencing using the Illumina HiSeq 2500 instrument. Each of the 48 cDNA libraries was sequenced at a read length of 100 base pairs.

Bioinformatics and statistical analysis

Each of the 48 sequenced libraries were processed using bioinformatics tools. First, the 100 base pair reads were trimmed to remove adapter sequences using scythe ( Next, 15 bases (sequencing artifacts) were removed from each read using FASTX trimmer ( Sickle ( was then used to remove bases with low quality scores and short reads. Using default settings, TopHat version 2.0.3 [109] was used to align the processed single end reads to the Williams 82 reference genome obtained from [42]. Samtools view [110] was used to remove unreliably mapped reads (mapping score < 1). The resulting mapping files (bamfiles) were then imported into the R statistical package (R Development Core Team, 2014; using the RSamtools statistical package in R ( The 48 bamfiles are available from the National Center for Biotechnology Sequence Read Archive under the BioProject accession number PRJNA478017.

The gene feature file for version 2 of the soybean genome was imported to R using rtracklayer [111] and the number of reads aligning to each gene for each sample was determined using GenomicAlignments [112]. Genes with counts per million < 1 were eliminated from further analysis. Data normalization across all treatments and genotypes was done for each time point using the Trimmed Mean of M (TMM) values [113] in the Bioconductor package edgeR [114]. Specifically, edgeR was used for single factor, pairwise comparisons to calculate normalization factors, estimate tagwise dispersion and determine differential gene expression. Initially, our data analysis pipeline used all samples in the experiment. Within the pipeline, we used the graphics package, ggplot2 [115] to visualize normalized gene expression of sample replicates for each treatment to ensure replicability. Due to considerable variation from other replicates, two samples from the Rag1/2 soybean genotype for the 6 h time point (one mock-treated sample and one aphid-treated sample) were removed from the analysis, leaving two replicates for each of these treatments for this genotype and time point. The design statement in the analysis pipeline was modified to reflect removal of these samples and the process from normalization to identification of DEGs was repeated. To identify genes differentially expressed in response to treatment, genotype or treatment x genotype interaction, our model took into account genotype and treatment (model.matrix(~genotype + trt, Design). To identify genes responding to treatment in a given genotype, our model grouped samples by type (model.matrix(~ 0 + Groups) and then we used contrast statements to facilitate comparisons. In all comparisons, differential expression for each gene was considered significant if false discovery rate (FDR) was < 0.05 (q-values < 0.05). Both treatment and genotype effects were examined in our statistical analysis, identifying genes that responded to aphid treatments in all four soybean genotypes and genotype-specific responses to soybean aphids. Genes with positive fold changes were induced by aphids while genes with negative fold changes were repressed.

Gene list comparisons and construction of Venn diagrams was done using Venny2.1 ( Cluster analysis was conducted to identify gene sets with similar expression patterns in response to genotype and treatment. Log counts per million for each sample were used for cluster analysis. DEGs identified at 6 and 12 H were included to examine expression across time, even though a given gene might only be significant at a single timepoint. The hclust package in R (, with the Pearson correlation and complete linkage methods, was used to construct the heatmap and cluster dendrogram, both of which showed similarities in expression patterns of DE genes across all comparisons. Z-scores for each DE gene were used during clustering.

To compare observed versus expected additive values for 6 h Rag1/2 synergistic genes, expected values were calculated as the sum of the logFC observed in the Rag1 and the Rag2 genotypes. Expected and observed values for the whole synergistic dataset or for induced and repressed subsets were compared using a χ2 test.

For comparison of the Rag1 locus between versions 1 and 2 of the Williams 82 genome assemblies, WebACT [116] was used to run BLASTN [117] alignments of the two regions, requiring 99% sequence identify for homology calls.

Annotation and GO enrichment of differentially expressed genes

Gene annotation was performed using the gene annotation lookup tool in SoyBase [118] ( The annotation result contained the best UniRef100 hit [119], the best homolog hit in Arabidopsis thaliana, and gene ontology (GO) information of the best Arabidopsis homolog (The Arabidopsis information resource, TAIR; The GO term enrichment tool in SoyBase ( was used to identify significantly overrepresented GO terms (biological process) in different gene sets. GO term enrichment was conducted using Fisher’s Exact Test and Bonferroni correction for multiple comparisons, identifying significantly overrepresented biological processes relative to the Williams 82 reference genome in each comparison (adjusted P ≤ 0.05).

Transcription factors (TF), TF expression network, and hormone signal analyses

To identify TF DE in the response to soybean aphids, we used the SoyDB transcription factor database [49]. TFF enrichment was conducted using Fisher’s Exact Test and Bonferroni correction for multiple comparisons, identifying significantly overrepresented TFFs relative to the Williams 82 reference genome in each comparison (Corrected P ≤ 0.05). Transcription factors that were involved in regulation of aphid-responsive biological processes were identified for each comparison for all soybean genotypes tested. Network analysis was performed using STRING with default parameters [120], but using the best Arabidopsis homolog for each soybean gene, identified in the annotation process. Manual curation of the output was carried out to annotate each gene with information of the datasets in which it was found to be differentially expressed, as well as the presence of multiple DE soybean genes corresponding to the same Arabidopsis gene. Changes in hormone signaling were evaluated using Hormonometer [45], also using the best Arabidopsis homolog for each soybean gene, identified in the annotation process, as previously described [46, 47].

Availability of data and materials

All data (raw reads and bamfiles) for this study are available in the National Center for Biotechnology Sequence Read Archive (SRA) under the BioProject accession number PRJNA478017.



bluegreen aphid


differentially expressed


effector-triggered immunity


fold change


false discovery rate


gene ontology




nucleotide-binding leucine-rich-repeat protein


PAMP-triggered immunity




reactive oxygen species




soybean cyst nematode


transcription factor


  1. Kourelis J, van der Hoorn RAL. Defended to the nines: 25 years of resistance gene cloning identifies nine mechanisms for R protein function. Plant Cell. 2018;30(2):285–99.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  2. Upson JL, Zess EK, Białas A, C-h W, Kamoun S. The coming of age of EvoMPMI: evolutionary molecular plant–microbe interactions across multiple timescales. Curr Opin Plant Biol. 2018;44:108–16.

    Article  PubMed  Google Scholar 

  3. Cui H, Tsuda K, Parker JE. Effector-triggered immunity: from pathogen perception to robust defense. Annu Rev Plant Biol. 2015;66(1):487–511.

    Article  PubMed  CAS  Google Scholar 

  4. Lachowiec J, Queitsch C, Kliebenstein DJ. Molecular mechanisms governing differential robustness of development and environmental responses in plants. Ann Bot. 2015;117(5):795–809.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  5. Wu C-H, Abd-El-Haliem A, Bozkurt TO, Belhaj K, Terauchi R, Vossen JH, et al. NLR network mediates immunity to diverse plant pathogens. Proc Natl Acad Sci. 2017;114(30):8113–8.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  6. Bentham AR, De la Concepcion JC, Mukhi N, Zdrzałek R, Draeger M, Gorenkin D, et al. A molecular roadmap to the plant immune system. J Biol Chem. 2020;295(44):14916–35.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. Bonardi V, Tang S, Stallmann A, Roberts M, Cherkis K, Dangl JL. Expanded functions for a family of plant intracellular immune receptors beyond specific recognition of pathogen effectors. Proc Natl Acad Sci. 2011;108(39):16463–8.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Windram O, Denby KJ. Modelling signaling networks underlying plant defence. Curr Opin Plant Biol. 2015;27:165–71.

    Article  PubMed  CAS  Google Scholar 

  9. Zhang W, Corwin JA, Copeland D, Feusier J, Eshbaugh R, Chen F, et al. Plastic transcriptomes stabilize immunity to pathogen diversity: the Jasmonic acid and salicylic acid networks within the Arabidopsis/Botrytis Pathosystem. Plant Cell. 2017;29(11):2727–52.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  10. Tsuda K, Somssich IE. Transcriptional networks in plant immunity. New Phytol. 2015;206(3):932–47.

    Article  PubMed  CAS  Google Scholar 

  11. Seo E, Choi D. Functional studies of transcription factors involved in plant defenses in the genomics era. Brief Funct Gen. 2015;14(4):260–7.

    Article  CAS  Google Scholar 

  12. Campbell MA, Fitzgerald HA, Ronald PC. Engineering pathogen resistance in crop plants. Transgenic Res. 2002;11(6):599–613.

    Article  PubMed  CAS  Google Scholar 

  13. Mundt CC. Durable resistance: a key to sustainable management of pathogens and pests. Infect Genet Evol. 2014;27:446–55.

    Article  PubMed  Google Scholar 

  14. Mundt CC. Pyramiding for resistance durability: theory and practice. Phytopathology. 2018;108(7):792–802.

    Article  PubMed  CAS  Google Scholar 

  15. MacIntosh GC. Gene pyramids and the balancing act of keeping pests at bay. J Exp Bot. 2019;70(18):4591–3.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  16. Kamphuis LG, Klingler JP, Jacques S, Gao L-L, Edwards OR, Singh KB. Additive and epistatic interactions between AKR and AIN loci conferring bluegreen aphid resistance and hypersensitivity in Medicago truncatula. J Exp Bot. 2019;70(18):4887–902.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Wiarda SL, Fehr WR, O'Neal ME. Soybean aphid (Hemiptera: Aphididae) development on soybean with Rag1 alone, Rag2 alone, and both genes combined. J Econ Entomol. 2012;105(1):252–8.

    Article  PubMed  CAS  Google Scholar 

  18. Gu K, Tian D, Qiu C, Yin Z. Transcription activator-like type III effector AvrXa27 depends on OsTFIIAγ5 for the activation of Xa27 transcription in rice that triggers disease resistance to Xanthomonas oryzae pv. Oryzae. Mol Plant Pathol. 2009;10(6):829–35.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  19. Ragsdale DW, Voegtlin DJ, O'Neil RJ. Soybean aphid biology in North America. Ann Entomol Soc Am. 2004;97(2):204–8.

    Article  Google Scholar 

  20. Tilmon KJ, Hodgson EW, O'Neal ME, Ragsdale DW. Biology of the soybean aphid, Aphis glycines (Hemiptera: Aphididae) in the United States. J Integr Pest Manag. 2011;2(2):A1–7.

    Article  Google Scholar 

  21. Neupane S, Purintun JM, Mathew FM, Varenhorst AJ, Nepal MP. Molecular basis of soybean resistance to soybean aphids and soybean cyst nematodes. Plants. 2019;8(10):374.

    Article  PubMed Central  CAS  Google Scholar 

  22. Li Y, Zou J, Li M, Bilgin DD, Vodkin LO, Hartman GL, et al. Soybean defense responses to the soybean aphid. New Phytol. 2008;179(1):185–95.

    Article  PubMed  CAS  Google Scholar 

  23. Studham ME, MacIntosh GC. Multiple phytohormone signals control the transcriptional response to soybean aphid infestation in susceptible and resistant soybean plants. Mol Plant-Microbe Interact. 2013;26(1):116–29.

    Article  PubMed  CAS  Google Scholar 

  24. Brechenmacher L, Nguyen THN, Zhang N, Jun T-H, Xu D, Mian MAR, et al. Identification of soybean proteins and genes differentially regulated in near isogenic lines differing in resistance to aphid infestation. J Proteome Res. 2015;14(10):4137–46.

    Article  PubMed  CAS  Google Scholar 

  25. Hill CB, Li Y, Hartman GL. Resistance to the soybean aphid in soybean germplasm. Crop Sci. 2004;44(1):98–106.

    Article  Google Scholar 

  26. Hill CB, Li Y, Hartman GL. A single dominant gene for resistance to the soybean aphid in the soybean cultivar Dowling. Crop Sci. 2006;46(4):1601–5.

    Article  Google Scholar 

  27. Hill CB, Li Y, Hartman GL. Soybean aphid resistance in soybean Jackson is controlled by a single dominant gene. Crop Sci. 2006;46(4):1606–8.

    Article  CAS  Google Scholar 

  28. Kim K-S, Bellendir S, Hudson K, Hill C, Hartman G, Hyten D, et al. Fine mapping the soybean aphid resistance gene Rag1 in soybean. TAG Theor Appl Genet. 2010;120(5):1063–71.

    Article  PubMed  CAS  Google Scholar 

  29. Mian MAR, Kang S-T, Beil S, Hammond R. Genetic linkage mapping of the soybean aphid resistance gene in PI 243540. TAG Theor Appl Genet. 2008;117(6):955–62.

    Article  CAS  Google Scholar 

  30. Hill CB, Kim K-S, Crull L, Diers BW, Hartman GL. Inheritance of resistance to the soybean aphid in soybean PI 200538. Crop Sci. 2009;49(4):1193–200.

    Article  CAS  Google Scholar 

  31. Kim K-S, Hill CB, Hartman GL, Hyten DL, Hudson ME, Diers BW. Fine mapping of the soybean aphid-resistance gene Rag2 in soybean PI 200538. Theor Appl Genet. 2010;121(3):599–610.

    Article  PubMed  CAS  Google Scholar 

  32. Hesler LS, Chiozza MV, O'Neal ME, MacIntosh GC, Tilmon KJ, Chandrasena DI, et al. Performance and prospects of rag genes for management of soybean aphid. Entomologia Exp et Applicata. 2013;147(3):201–16.

    Article  CAS  Google Scholar 

  33. Kim K-S, Hill CB, Hartman GL, Mian MAR, Diers BW. Discovery of soybean aphid biotypes. Crop Sci. 2008;48(3):923–8.

    Article  Google Scholar 

  34. Hill CB, Crull L, Herman TK, Voegtlin DJ, Hartman GL. A new soybean aphid (Hemiptera: Aphididae) biotype identified. J Econ Entomol. 2010;103(2):509–15.

    Article  PubMed  CAS  Google Scholar 

  35. Alt J, Ryan-Mahmutagic M. Soybean aphid biotype 4 identified. Crop Sci. 2013;53(4):1491–5.

    Article  Google Scholar 

  36. Ajayi-Oyetunde OO, Diers BW, Lagos-Kutz D, Hill CB, Hartman GL, Reuter-Carlson U, et al. Differential Reactions of Soybean Isolines With Combinations of Aphid Resistance Genes Rag1 , Rag2 , and Rag3 to Four Soybean Aphid Biotypes. J Econ Entomol. 2016;109(3):1431–7.

    Article  PubMed  CAS  Google Scholar 

  37. Zhang S, Wen Z, DiFonzo C, Song Q, Wang D. Pyramiding different aphid-resistance genes in elite soybean germplasm to combat dynamic aphid populations. Mol Breed. 2018;38(3):29.

    Article  CAS  Google Scholar 

  38. Varenhorst AJ, Pritchard SR, O’Neal ME, Hodgson EW, Singh AK. Determining the effectiveness of three-gene pyramids against Aphis glycines (Hemiptera: Aphididae) biotypes. J Econ Entomol. 2017;110(6):2428–35.

    Article  PubMed  CAS  Google Scholar 

  39. McCarville MT, O'Neal ME. Measuring the benefit of biological control for single gene and Pyramided host plant resistance for Aphis glycines (Hemiptera: Aphididae) management. J Econ Entomol. 2012;105(5):1835–43.

    Article  PubMed  CAS  Google Scholar 

  40. Diaz JEL, Ahsen ME, Schaffter T, Chen X, Realubit RB, Karan C, et al. The transcriptomic response of cells to a drug combination is more than the sum of the responses to the monotherapies. eLife. 2020;9:e52707.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Hohenstein JD, Studham ME, Klein A, Kovinich N, Barry K, Lee Y-J, et al. Transcriptional and chemical changes in soybean leaves in response to Long-term aphid colonization. Front Plant Sci. 2019;10.

  42. Schmutz J, Cannon SB, Schlueter J, Ma J, Mitros T, Nelson W, et al. Genome sequence of the palaeopolyploid soybean. Nature. 2010;463(7278):178–83.

    Article  PubMed  CAS  Google Scholar 

  43. Lee S, Cassone BJ, Wijeratne A, Jun T-H, Michel AP, Mian MAR. Transcriptomic dynamics in soybean near-isogenic lines differing in alleles for an aphid resistance gene, following infestation by soybean aphid biotype 2. BMC Genomics. 2017;18(1):472.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. O’Rourke JA, Graham MA. Gene expression responses to sequential nutrient deficiency stresses in soybean. Int J Mol Sci. 2021;22(3):1252.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  45. Volodarsky D, Leviatan N, Otcheretianski A, Fluhr R. HORMONOMETER: a tool for discerning transcript signatures of hormone action in the Arabidopsis transcriptome. Plant Physiol. 2009;150(4):1796–805.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  46. Tzin V, Fernandez-Pozo N, Richter A, Schmelz EA, Schoettner M, Schäfer M, et al. Dynamic maize responses to aphid feeding are revealed by a time series of transcriptomic and Metabolomic assays. Plant Physiol. 2015;169(3):1727–43.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  47. Bui H, Greenhalgh R, Ruckert A, Gill GS, Lee S, Ramirez RA, et al. Generalist and specialist mite herbivores induce similar defense responses in maize and barley but differ in susceptibility to Benzoxazinoids. Front Plant Sci. 2018;9.

  48. Yates-Stewart AD, Pekarcik A, Michel A, Blakeslee JJ. Jasmonic acid-isoleucine (JA-Ile) is involved in the host-plant resistance mechanism against the soybean aphid (Hemiptera: Aphididae). J Econ Entomol. 2020;113(6):2972–8.

    Article  PubMed  CAS  Google Scholar 

  49. Wang Z, Libault M, Joshi T, Valliyodan B, Nguyen HT, Xu D, et al. SoyDB: a knowledge database of soybean transcription factors. BMC Plant Biol. 2010;10(1):14.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  50. Thireault C, Shyu C, Yoshida Y, St. Aubin B, Campos ML, Howe GA: Repression of jasmonate signaling by a non-TIFY JAZ protein in Arabidopsis. Plant J 2015, 82:669–679, 4,

  51. Li Z, Luo X, Ou Y, Jiao H, Peng L, Fu X, et al. JASMONATE-ZIM DOMAIN proteins engage Polycomb chromatin modifiers to modulate Jasmonate signaling in Arabidopsis. Mol Plant. 2021;14(5):732–47.

    Article  PubMed  CAS  Google Scholar 

  52. Pauwels L, Goossens A. The JAZ proteins: a crucial Interface in the Jasmonate signaling Cascade. Plant Cell. 2011;23(9):3089–100.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  53. Kazan K, Manners JM. JAZ repressors and the orchestration of phytohormone crosstalk. Trends Plant Sci. 2012;17(1):22–31.

    Article  PubMed  CAS  Google Scholar 

  54. Shi A, Chen P, Li D, Zheng C, Zhang B, Hou A. Pyramiding multiple genes for resistance to soybean mosaic virus in soybean using molecular markers. Mol Breed. 2008;23(1):113–24.

    Article  CAS  Google Scholar 

  55. Huang N, Angeles ER, Domingo J, Magpantay G, Singh S, Zhang G, et al. Pyramiding of bacterial blight resistance genes in rice: marker-assisted selection using RFLP and PCR. Theor Appl Genet. 1997;95(3):313–20.

    Article  CAS  Google Scholar 

  56. Suh J-P, Jeung J-U, Noh T-H, Cho Y-C, Park S-H, Park H-S, et al. Development of breeding lines with three pyramided resistance genes that confer broad-spectrum bacterial blight resistance and their molecular analysis in rice. Rice. 2013;6(1):5.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Hittalmani S, Parco A, Mew TV, Zeigler RS, Huang N. Fine mapping and DNA marker-assisted pyramiding of the three major genes for blast resistance in rice. Theor Appl Genet. 2000;100(7):1121–8.

    Article  CAS  Google Scholar 

  58. Kim H-J, Lee H-R, Jo K-R, Mortazavian SMM, Huigen DJ, Evenhuis B, et al. Broad spectrum late blight resistance in potato differential set plants MaR8 and MaR9 is conferred by multiple stacked R genes. Theor Appl Genet. 2012;124(5):923–35.

    Article  PubMed  CAS  Google Scholar 

  59. Jo K-R, Kim C-J, Kim S-J, Kim T-Y, Bergervoet M, Jongsma MA, et al. Development of late blight resistant potatoes by cisgene stacking. BMC Biotechnol. 2014;14(1):50.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Kloppers FJ, Pretorius ZA. Effects of combinations amongst genes Lr13, Lr34 and Lr37 on components of resistance in wheat to leaf rust. Plant Pathol. 1997;46(5):737–50.

    Article  Google Scholar 

  61. Koller T, Brunner S, Herren G, Hurni S, Keller B. Pyramiding of transgenic Pm3 alleles in wheat results in improved powdery mildew resistance in the field. Theor Appl Genet. 2018;131(4):861–71.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  62. Gao L, Fang Z, Zhou J, Li L, Lu L, Li L, et al. Transcriptional insights into the pyramided resistance to rice bacterial blight. Sci Rep. 2018;8(1):12358.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  63. Abeysinghe JK, Lam K-M, Ng DW-K. Differential regulation and interaction of homoeologous WRKY18 and WRKY40 in Arabidopsis allotetraploids and biotic stress responses. Plant J. 2019;97(2):352–67.

    Article  PubMed  CAS  Google Scholar 

  64. Wang Y, Wang H, Ma Y, Du H, Yang Q, Yu D. Identification of transcriptional regulatory nodes in soybean defense networks using transient co-transactivation assays. Front Plant Sci. 2015;6.

  65. Moreau M, Degrave A, Vedel R, Bitton F, Patrit O, Renou J-P, et al. EDS1 Contributes to Nonhost Resistance of Arabidopsis thaliana Against Erwinia amylovora. Mol Plant-Microbe Int. 2012;25:421–30.

    Article  CAS  Google Scholar 

  66. Amrine KCH, Blanco-Ulate B, Cantu D. Discovery of Core biotic stress responsive genes in Arabidopsis by weighted gene co-expression network analysis. PLoS ONE. 2015;10(3):e0118731.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  67. Chen X, Li C, Wang H, Guo Z. WRKY transcription factors: evolution, binding, and action. Phytopatho Res. 2019;1(1):13.

    Article  Google Scholar 

  68. Bakshi M, Oelmüller R. WRKY transcription factors. Plant Signal Behav. 2014;9(2):e27700.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  69. Liu J, Osbourn A, Ma P. MYB transcription factors as regulators of Phenylpropanoid metabolism in plants. Mol Plant. 2015;8(5):689–708.

    Article  PubMed  CAS  Google Scholar 

  70. Foyer CH, Verrall SR, Hancock RD. Systematic analysis of phloem-feeding insect-induced transcriptional reprogramming in Arabidopsis highlights common features and reveals distinct responses to specialist and generalist insects. J Exp Bot. 2015;66(2):495–512.

    Article  PubMed  CAS  Google Scholar 

  71. Barah P, Winge P, Kusnierczyk A, Tran DH, Bones AM. Molecular signatures in Arabidopsis thaliana in response to insect attack and bacterial infection. PLoS ONE. 2013;8(3):e58987.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  72. Liu R, Lü B, Wang X, Zhang C, Zhang S, Qian J, et al. Thirty-seven transcription factor genes differentially respond to a harpin protein and affect resistance to the green peach aphid in Arabidopsis. J Biosci. 2010;35(3):435–50.

    Article  PubMed  CAS  Google Scholar 

  73. Van Eck L, Schultz T, Leach JE, Scofield SR, Peairs FB, Botha A-M, et al. Virus-induced gene silencing of WRKY53 and an inducible phenylalanine ammonia-lyase in wheat reduces aphid resistance. Plant Biotechnol J. 2010;8(9):1023–32.

    Article  PubMed  CAS  Google Scholar 

  74. Zhang H, Kjemtrup-Lovelace S, Li C, Luo Y, Chen LP, Song B-H. Comparative RNA-Seq analysis uncovers a complex regulatory network for soybean cyst nematode resistance in wild soybean (Glycine soja). Sci Rep. 2017;7(1):9699.

    Article  PubMed  PubMed Central  Google Scholar 

  75. Yang Y, Zhou Y, Chi Y, Fan B, Chen Z. Characterization of soybean WRKY gene family and identification of soybean WRKY genes that promote resistance to soybean cyst nematode. Sci Rep. 2017;7(1):17804.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  76. Goggin FL. Plant-aphid interactions: molecular and ecological perspectives. Curr Opin Plant Biol. 2007;10(4):399–408.

    Article  PubMed  CAS  Google Scholar 

  77. Silva-Sanzana C, Estevez JM, Blanco-Herrera F. Influence of cell wall polymers and their modifying enzymes during plant–aphid interactions. J Exp Bot. 2020;71(13):3854–64.

    Article  PubMed  CAS  Google Scholar 

  78. Wang Y, Sheng L, Zhang H, Du X, An C, Xia X, et al. CmMYB19 over-expression improves aphid tolerance in Chrysanthemum by promoting lignin synthesis. Int J Mol Sci. 2017;18(3):619.

    Article  PubMed Central  CAS  Google Scholar 

  79. Chezem WR, Memon A, Li F-S, Weng J-K, Clay NK. SG2-type R2R3-MYB transcription factor MYB15 controls defense-induced lignification and basal immunity in Arabidopsis. Plant Cell. 2017;29(8):1907–26.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  80. An C, Sheng L, Du X, Wang Y, Zhang Y, Song A, et al. Overexpression of CmMYB15 provides chrysanthemum resistance to aphids by regulating the biosynthesis of lignin. Hortic Res. 2019;6(1):84.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  81. McCarthy RL, Zhong R, Ye Z-H. MYB83 is a direct target of SND1 and acts redundantly with MYB46 in the regulation of secondary Cell Wall biosynthesis in Arabidopsis. Plant Cell Physiol. 2009;50(11):1950–64.

    Article  PubMed  CAS  Google Scholar 

  82. Zhong R, Ye Z-H. MYB46 and MYB83 bind to the SMRE sites and directly activate a suite of transcription factors and secondary wall biosynthetic genes. Plant Cell Physiol. 2011;53(2):368–80.

    Article  PubMed  CAS  Google Scholar 

  83. Naoumkina MA, He X, Dixon RA. Elicitor-induced transcription factors for metabolic reprogramming of secondary metabolism in Medicago truncatula. BMC Plant Biol. 2008;8(1):132.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  84. Meng F, Han Y, Teng W, Li Y, Li W. QTL underlying the resistance to soybean aphid (Aphis glycines Matsumura) through isoflavone-mediated antibiosis in soybean cultivar ‘Zhongdou 27’. Theor Appl Genet. 2011;123(8):1459–65.

    Article  PubMed  CAS  Google Scholar 

  85. Chu S, Wang J, Zhu Y, Liu S, Zhou X, Zhang H, et al. An R2R3-type MYB transcription factor, GmMYB29, regulates isoflavone biosynthesis in soybean. PLoS Genet. 2017;13(5):e1006770.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  86. Jahan MA, Harris B, Lowery M, Infante AM, Percifield RJ, Kovinich N. Glyceollin transcription factor GmMYB29A2 regulates soybean resistance to Phytophthora sojae. Plant Physiol. 2020;183(2):530–46.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  87. Van Aken O, Zhang B, Law S, Narsai R, Whelan J. AtWRKY40 and AtWRKY63 modulate the expression of stress-responsive nuclear genes encoding mitochondrial and chloroplast proteins. Plant Physiol. 2013;162(1):254–71.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  88. Liu R, Xu Y-H, Jiang S-C, Lu K, Lu Y-F, Feng X-J, et al. Light-harvesting chlorophyll a/b-binding proteins, positively involved in abscisic acid signalling, require a transcription repressor, WRKY40, to balance their function. J Exp Bot. 2013;64(18):5443–56.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  89. Bi Y-M, Zhang Y, Signorelli T, Zhao R, Zhu T, Rothstein S. Genetic analysis of Arabidopsis GATA transcription factor gene family reveals a nitrate-inducible member important for chlorophyll synthesis and glucose sensitivity. Plant J. 2005;44(4):680–92.

    Article  PubMed  CAS  Google Scholar 

  90. Hudson D, Guevara D, Yaish MW, Hannam C, Long N, Clarke JD, et al. GNC and CGA1 modulate chlorophyll biosynthesis and glutamate synthase (GLU1/Fd-GOGAT) expression in Arabidopsis. PLoS ONE. 2011;6(11):e26765.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  91. Xu Z, Casaretto JA, Bi Y-M, Rothstein SJ. Genome-wide binding analysis of AtGNC and AtCGA1 demonstrates their cross-regulation and common and specific functions. Plant Direct. 2017;1(4):e00016.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  92. McCormac AC, Terry MJ. Light-signalling pathways leading to the co-ordinated expression of HEMA1 and Lhcb during chloroplast development in Arabidopsis thaliana. Plant J. 2002;32(4):549–59.

    Article  PubMed  CAS  Google Scholar 

  93. Wang P, Hendron R-W, Kelly S. Transcriptional control of photosynthetic capacity: conservation and divergence from Arabidopsis to rice. New Phytol. 2017;216(1):32–45.

    Article  PubMed  CAS  Google Scholar 

  94. Berger S, Benediktyová Z, Matouš K, Bonfig K, Mueller MJ, Nedbal L, et al. Visualization of dynamics of plant–pathogen interaction by novel combination of chlorophyll fluorescence imaging and statistical analysis: differential effects of virulent and avirulent strains of P. syringae and of oxylipins on a. thaliana. J Exp Bot. 2006;58(4):797–806.

    Article  PubMed  CAS  Google Scholar 

  95. Bonfig KB, Schreiber U, Gabler A, Roitsch T, Berger S. Infection with virulent and avirulent P. syringae strains differentially affects photosynthesis and sink metabolism in Arabidopsis leaves. Planta. 2006;225(1):1–12.

    Article  PubMed  CAS  Google Scholar 

  96. Zou J, Rodriguez-Zas S, Aldea M, Li M, Zhu J, Gonzalez DO, et al. Expression profiling soybean response to pseudomonas syringae reveals new defense-related genes and rapid HR-specific downregulation of photosynthesis. Mol Plant-Microbe Interact. 2005;18(11):1161–74.

    Article  PubMed  CAS  Google Scholar 

  97. Su J, Yang L, Zhu Q, Wu H, He Y, Liu Y, et al. Active photosynthetic inhibition mediated by MPK3/MPK6 is critical to effector-triggered immunity. PLoS Biol. 2018;16(5):e2004122.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  98. Kretschmer M, Damoo D, Djamei A, Kronstad J. Chloroplasts and plant immunity: where are the fungal effectors. Pathogens. 2020;9(1):19.

    Article  CAS  Google Scholar 

  99. Miller G, Schlauch K, Tam R, Cortes D, Torres MA, Shulaev V, et al. The Plant NADPH Oxidase RBOHD Mediates Rapid Systemic Signaling in Response to Diverse Stimuli. Sci Signal. 2009;2:ra45.

    Article  PubMed  Google Scholar 

  100. Moloi MJ, van der Westhuizen AJ. The reactive oxygen species are involved in resistance responses of wheat to the Russian wheat aphid. J Plant Physiol. 2006;163(11):1118–25.

    Article  PubMed  CAS  Google Scholar 

  101. Kerchev PI, Fenton B, Foyer CH, Hancock RD. Infestation of potato (Solanum tuberosum L.) by the peach-potato aphid (Myzus persicae Sulzer) alters cellular redox status and is influenced by ascorbate. Plant, Cell Environ. 2012;35:430–40.

    Article  CAS  Google Scholar 

  102. Rasool B, Karpinska B, Konert G, Durian G, Denessiouk K, Kangasjärvi S, et al. Effects of light and the regulatory B-subunit composition of protein phosphatase 2A on the susceptibility of Arabidopsis thaliana to aphid (Myzus persicae) infestation. Front Plant Sci. 2014;5.

  103. Niu L, Pan L, Zeng W, Lu Z, Cui G, Fan M, et al. Dynamic transcriptomes of resistant and susceptible peach lines after infestation by green peach aphids (Myzus persicae Sülzer) reveal defence responses controlled by the Rm3 locus. BMC Genomics. 2018;19(1):846.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  104. Coates BS, Hohenstein JD, Giordano R, Donthu RK, Michel AP, Hodgson EW, et al. Genome scan detection of selective sweeps among biotypes of the soybean aphid, Aphis glycines, with differing virulence to resistance to A. glycines (Rag) traits in soybean, Glycine max. Insect Biochem Mol Biol. 2020;124:103364.

    Article  PubMed  CAS  Google Scholar 

  105. Brace RC, Fehr WR. Impact of combining the Rag1 and Rag2 alleles for aphid resistance on agronomic and seed traits of soybean. Crop Sci. 2012;52(5):2070–4.

    Article  Google Scholar 

  106. Varenhorst AJ, McCarville MT, O’Neal ME. Reduced fitness of virulent Aphis glycines (Hemiptera: Aphididae) biotypes may influence the longevity of resistance genes in soybean. PLoS ONE. 2015;10(9):e0138252.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  107. Paz MM, Martinez JC, Kalvig AB, Fonger TM, Wang K. Improved cotyledonary node method using an alternative explant derived from mature seed for efficient Agrobacterium-mediated soybean transformation. Plant Cell Rep. 2006;25(3):248.

    Article  CAS  Google Scholar 

  108. Fehr WR, Caviness CE. Stages of soybean development. Cooperative Extension Service, Agriculture and Home Economic Experiment Station, Iowa State University. Special Rep. 1977:80.

  109. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  110. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  111. Lawrence M, Gentleman R, Carey V. rtracklayer: an R package for interfacing with genome browsers. Bioinformatics. 2009;25(14):1841–2.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  112. Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, Gentleman R, et al. Software for computing and annotating genomic ranges. PLoS Comput Biol. 2013;9(8):e1003118.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  113. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3):R25.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  114. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  PubMed  CAS  Google Scholar 

  115. Wickham H: ggplot2: elegant graphics for data analysis. Verlag New York: Springer; 2009.

  116. Abbott JC, Aanensen DM, Rutherford K, Butcher S, Spratt BG. WebACT—an online companion for the Artemis comparison tool. Bioinformatics. 2005;21(18):3665–6.

    Article  PubMed  CAS  Google Scholar 

  117. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.

    Article  CAS  PubMed  Google Scholar 

  118. Grant D, Nelson RT, Cannon SB, Shoemaker RC. SoyBase, the USDA-ARS soybean genetics and genomics database. Nucleic Acids Res. 2010;38(suppl_1):D843–6.

    Article  PubMed  CAS  Google Scholar 

  119. The UniProt Consortium. UniProt: the universal protein knowledgebase. Nucleic Acids Res. 2017;45(D1):D158–69.

    Article  CAS  Google Scholar 

  120. von Mering C, Jensen LJ, Snel B, Hooper SD, Krupp M, Foglierini M, et al. STRING: known and predicted protein–protein associations, integrated and transferred across organisms. Nucleic Acids Res. 2005;33(Database issue):D433–7.

    Article  CAS  Google Scholar 

Download references


The authors appreciate the undergraduate research assistants Savanah McDermott, Victoria Ridout, and Siau Ting Lau for helping with sample collection and tissue processing. We also thank Dr. Walter Fehr for providing seed for the four soybean genotypes used in this study, Dr. Jamie O’Rouke for providing laboratory space for RNA preparation and storage, and Dr. Adrienne Moran Lauter for her advice on RNA preparation protocols.


This study was funded by the Soybean Checkoff through a grant from the Iowa Soybean Association, a U.S. Department of Agriculture (USDA) AFRI/NIFA grant (# 2019–67013-29352), a grant from the Roy J. Carver Charitable Trust (Muscatine, Iowa), and by the USDA Agricultural Research Service Project 5030–21220-006-00D. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the USDA. The findings and conclusions in this are those of the authors and should not be construed to represent any official USDA or U.S. Government determination or policy. USDA is an equal opportunity provider and employer. Funding bodies had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Authors and Affiliations



MIN, JDH, MAG, AKS, and GCM conceived and designed the experiments. MIN, JDH, CEM, MAG, and GCM carried out the data generation. MIN, JDH, CEM, MAG, YQ and GCM performed the data analysis. MIN, CEM, MAG and GCM prepared the figures. MIN and GCM wrote the manuscript. All authors revised the manuscript draft and read and approved the final manuscript.

Corresponding author

Correspondence to Gustavo C. MacIntosh.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional File 1.

Datasets corresponding to all differentially expressed genes in response to aphids for the two time points analyzed (Tab 1: full dataset; tab2: 6 h synergistic response)

Additional File 2.

List of differentially expressed genes in comparisons among mock-treatments for all genotypes

Additional File 3.

Supplemental Figs. S1, S2 and S3. Fig. S1. Comparison of expected versus observed log2FC for the Rag1/2 6 h synergistic genes. Fig. S2. Comparison of the transcriptional response of each genotype at different time points. Fig. S3. Plant hormone signatures inferred by Hormonometer from the transcriptome responses of each genotype

Additional File 4.

List of overrepresented GO-terms (biological processes) for each genotype response to aphids for the two time points analyzed

Additional File 5.

List of WRKY transcription factors differentially expressed in response to aphids in all genotypes

Additional File 6.

Analysis of the enrichment of TF families in the 6 h synergistic response

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Natukunda, M.I., Hohenstein, J.D., McCabe, C.E. et al. Interaction between Rag genes results in a unique synergistic transcriptional response that enhances soybean resistance to soybean aphids. BMC Genomics 22, 887 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Gene pyramiding
  • Soybean
  • Soybean aphids
  • Aphid resistance
  • Rag genes
  • RNA sequencing
  • Synergistic effect