Molecular mechanisms underpinning quantitative resistance to Phytophthora sojae in Glycine max using a systems genomics approach

Expression of quantitative disease resistance in many host–pathogen systems is controlled by genes at multiple loci, each contributing a small effect to the overall response. We used a systems genomics approach to study the molecular underpinnings of quantitative disease resistance in the soybean-Phytophthora sojae pathosystem, incorporating expression quantitative trait loci (eQTL) mapping and gene co-expression network analysis to identify the genes putatively regulating transcriptional changes in response to inoculation. These findings were compared to previously mapped phenotypic (phQTL) to identify the molecular mechanisms contributing to the expression of this resistance. A subset of 93 recombinant inbred lines (RILs) from a Conrad × Sloan population were inoculated with P. sojae isolate 1.S.1.1 using the tray-test method; RNA was extracted, sequenced, and the normalized read counts were genetically mapped from tissue collected at the inoculation site 24 h after inoculation from both mock and inoculated samples. In total, more than 100,000 eQTLs were mapped. There was a switch from predominantly cis-eQTLs in the mock treatment to an almost entirely nonoverlapping set of predominantly trans-eQTLs in the inoculated treatment, where greater than 100-fold more eQTLs were mapped relative to mock, indicating vast transcriptional reprogramming due to P. sojae infection occurred. The eQTLs were organized into 36 hotspots, with the four largest hotspots from the inoculated treatment corresponding to more than 70% of the eQTLs, each enriched for genes within plant–pathogen interaction pathways. Genetic regulation of trans-eQTLs in response to the pathogen was predicted to occur through transcription factors and signaling molecules involved in plant–pathogen interactions, plant hormone signal transduction, and MAPK pathways. Network analysis identified three co-expression modules that were correlated with susceptibility to P. sojae and associated with three eQTL hotspots. Among the eQTLs co-localized with phQTLs, two cis-eQTLs with putative functions in the regulation of root architecture or jasmonic acid, as well as the putative master regulators of an eQTL hotspot nearby a phQTL, represent candidates potentially underpinning the molecular control of these phQTLs for resistance.

Few studies have examined the physiological and molecular mechanisms of quantitative resistance in soybean to P. sojae; however, all have concluded that this is a complex trait.Although only a few mechanisms for QDR have been substantiated in plant systems (Nelson et al., 2018), numerous hypotheses have been developed for resistance to P. sojae in soybeans from several previous studies.These hypotheses include plant hormone signal transduction, including auxin acting as a susceptibility factor (Wang et al., 2010;Wang et al., 2012b;Stasko et al., 2020); suberin playing a role in slowing P. sojae hyphal infection in epidermal walls and middle lamellae (Thomas et al., 2007;Ranathunge et al., 2008); and the phenylpropanoid pathway acting as a positive regulator of P. sojae infection by increased content of glyceollin, daidzein, genistein, and salicylic acid (SA) (Abbasi et al., 2001;Mohr and Cahill, 2001;Graham et al., 2003;Lygin et al., 2013;Zhang et al., 2017).Components of the isoflavonoid pathway have been implicated in acting as antioxidants to reduce reactive oxygen species (ROS) and enhance QDR to P. sojae (Xu et al., 2012;Wong et al., 2014;Cheng et al., 2015;Dastmalchi et al., 2017).QDR in P. sojae has also been associated with increased expression of genes coding for pathogenesis-related 1a protein (PR1a), matrix metalloproteinases, basic peroxidases, and b-1,3-endoglucanases (Vega-Sańchez et al., 2005), as well as ubiquitination, plant cell structural modifications, serine-threonine kinase, and basal resistance (Wang et al., 2012b;Karhoff et al., 2022).Recently, a gene annotated as a major latex protein expressed in the roots and associated with biotic stress was implicated in QDR (de Ronne et al., 2022).In addition to these pathways, there are other well-documented pathways involved in plant defense against pathogens, including mitogen-activated protein kinase (MAPK) cascades and plant-pathogen interaction pathways (e.g., pathogen-associated molecular patterns (PAMP) and effectortriggered immunity (ETI) pathways), which may also play a role in QDR (Ausubel, 2005;Glazebrook, 2005;Chisholm et al., 2006;Jones and Dangl, 2006;Boller and Felix, 2009;Dodds and Rathjen, 2010;Gassmann and Bhattacharjee, 2012;Spoel and Dong, 2012;Zipfel, 2014;Cui et al., 2015;Nelson et al., 2018).Taken together, these reports emphasize the complexity of QDR mechanisms and responses to infection.Zhou et al. (2009) showed that by 5 days postinoculation with P. sojae, from tissue collected in front of the advancing lesion, 97% of the genes in the soybean genome responded to infection or genetic variation based on microarray data.Amid this genome-wide transcriptional reprogramming, it is difficult to identify the specific causal mechanisms, pathways, and putative candidate genes.Thus, a more robust approach is needed to elucidate the molecular mechanisms underlying QDR.
Numerous studies in plant-pathogen interactions have utilized a systems genomics approach to identify the molecular mechanisms that are controlled by a complex of genes, pathways, and networks contributing to the overall expression of a phenotype (Jansen and Nap, 2001;Kliebenstein, 2009;Druka et al., 2010;Feltus, 2014).This approach maps both phQTLs and expression (eQTLs) and combines these data with a gene network analysis to identify the specific alleles that control or contribute to the overall expression of resistance during a plant-pathogen interaction (Kliebenstein, 2009;Druka et al., 2010;Feltus, 2014).Using this approach, advancements were made towards understanding the mechanisms of QDR resistance in barley (Hordeum vulgare L.) to Puccinia hordei (Druka et al., 2008;Chen et al., 2010) and in maize (Zea mays L.) to Cercospora zeina (Christie et al., 2017).In the barley-Puccinia hordei system, the total number of candidate genes underlying phQTLs was reduced at four different loci, with the identification of a histidine kinase as a novel resistance gene at one locus (Druka et al., 2008).In a later study, the total number of candidate genes for QDR to Puccinia hordei was reduced to six candidates for barley Rphq11 (Chen et al., 2010).Co-expression of coronatine-insensitive 1 (COI1) and jasmonate responses in maize were correlated with resistance to C. zeina, while pathogen manipulation of the host plant through the diterpenoid biosynthesis pathway was associated with susceptibility (Christie et al., 2017).
Due to the nature of the P. sojae-soybean pathosystem response, in which a small proportion of the total phenotype is contributed by each locus and a large number of loci encompassing numerous genes respond to P. sojae infection, we have taken a systems genomics approach to elucidate the molecular mechanisms of QDR in soybean toward P. sojae.In this study, we aimed to (1) understand the transcriptional reprogramming that occurs earlier in the infection process during the transition from biotrophic to necrotrophic between the well-studied soybean cultivars, Conrad, which has high levels of QDR, and Sloan, which is moderately susceptible; (2) map the genetic control of the differential transcriptional response to inoculation with P. sojae; (3) identify functional enrichment of genes within eQTL hotspots and coexpression modules associated with disease; and (4) identify candidate genes.The emphasis in this study is placed on identifying factors that elicit expression of quantitative resistance and not those expressed during the R-gene (Rps gene) response, which have also been recently studied using a transcriptomic approach (Lin et al., 2014;Hale et al., 2023a).
Three types of analyses were used, including phQTL mapping, eQTL mapping, and weighted gene co-expression analysis (WGCNA).The integration of the three analyses allowed for a greater understanding of the relationships between gene expression and the resulting disease phenotypes.Ultimately, using the expression data of 93 RILs during the early stages of the infection process [24 hours after inoculation (hai)] as the switch from biotrophy to necrotrophy is occurring (Moy et al., 2004), we were able to map more than 100,000 eQTLs, identifying co-regulated gene modules associated with disease and five putative candidate genes for three phQTLs.Putative master regulators were identified for 16 key eQTL hotspots.This study is the first to our knowledge to elucidate the specific candidate genes that may regulate the extensive changes that occur during transcriptional reprogramming due to pathogen infection by utilizing eQTL mapping in inoculated and non-inoculated tissues within the same population.

Phenotyping and RNA extraction
A subset of 93 F 9:11 RILs from the full population of 316 RILs derived from a cross between the cultivars Conrad (resistant) and Sloan (susceptible) were selected randomly for the eQTL study.The parents and RILs were inoculated with P. sojae 1.S.1.1 as previously described (Wang et al., 2012b).Briefly, the roots of 7-day-old plants that were grown for one week in 29.5-ml Styrofoam cups containing vermiculite (Perlite Vermiculite Packing Industries, Inc., North Bloomfield, OH, USA) were washed in tap water.Seedlings were placed on a plastic tray on top of a cotton wicking pad and a polyester cloth.A wound on the main tap root of each plant was made with a scalpel approximately 2.5 cm below the crown and covered with a mycelial slurry of seven-day-old 1.S.1.1 of P. sojae grown on lima bean (Phaseolus lunatus L.) agar.These trays were placed in buckets and kept in a growth chamber at 25°C, 20% relative humidity (RH), on a 14-h light:10-h dark cycle.The experimental design was a randomized incomplete block design with three replications with each block containing 50-100 RILs.There were three trays for each RIL in each block.The first two trays were for RNA isolation where each tray was either inoculated with 1.S.1.1 or mock inoculated with lima bean agar (without P. sojae).
No RNA was isolated from the third tray which was inoculated with 1.S.1.1,and the lesion length was measured 7 days after inoculation (dai) from the top of the inoculation site to the leading edge of the lesion margin to ensure that inoculations had been successful.A schematic for the experimental design is shown in Supplementary Figure S1.
All agar was gently removed using a Kimwipe, and a 1-cm section of tissue from the inoculation site was collected 24 hai from the mock and inoculated trays.Tissue from eight to 10 plants from a single RIL was collected, and tissue from all three replicates of each RIL was pooled for RNA extraction for a total of ~24-30 plants per RIL.Tissue was ground in liquid nitrogen, and RNA was extracted using a Qiagen RNeasy Plant Mini Kit (Hilden, Germany) following the manufacturer's protocols.
RNA was quality and quantity checked using the TapeStation (Agilent Technologies, Santa Clara, CA, USA) and Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA), respectively.RNA quality and quantity standards were met for both the mockinoculated and inoculated treatment for each of the 93 RILs.Sequence data from this study can be found in the NCBI BioProject: PRJNA478334.

Genotypic and phenotypic data
Phenotypic QTLs have been previously reported in the Conrad × Sloan recombinant inbred population of 316 individuals (Stasko et al., 2016), a subset of which was used for this study.The methods for inoculation, experimental design, and statistical analysis for mapping were described in greater detail (Wang et al., 2010;Wang et al., 2012b;Stasko et al., 2016).To construct a genetic map to carry out both eQTL and phQTL mapping, genotypic data was obtained from Stasko et al. (2016).While included in WGCNA, RIL 12280 was removed from the phQTL and eQTL analyses due to missing genotypic data.The normality of the data was evaluated using the Shapiro-Wilk test and visually assessed by histograms and QQ plots.
To account for environment variation, best linear unbiased predictor (BLUP) values for RILs were extracted from the model by adjusting mean lesion lengths of 92 RILs to the checks (cultivars OX20-8 (Rps1a, no partial resistance), Williams 82 (Rps1k, moderate partial resistance), Conrad (parent, high partial resistance), Sloan (parent, moderately susceptible)) and removing environmental variation in R version 3.5.1 (R Core Team, 2018) with the package -'lme4' using the lmer function in version 1.1-19 (Bates et al., 2014).The model applied was Y ijk = µ + R i + T j + G(C) jk + ϵ ijk where µ = overall mean, R i = effect of the ith replication, T j = effect of the jth type of entry, G(C) jk = effect of the kth genotype within the type for RIL only, and ϵ ijk = experimental error.The RILs in this study represent a subset of all possible genetic combinations and were treated as random (Green and Tukey, 1960).The same principle was applied to the random effect of replication, while the assigned category of type of entry (a category for RILs as well as each of the checks) was treated as fixed.Variance components were estimated using the maximum likelihood method.Composite interval mapping (R/qtl; Broman et al., 2003) was utilized to confirm phQTL peaks that were associated with lesion length BLUP values in the 316 RIL population from Stasko et al. (2016).

RNA sequencing
For each of the 184 samples consisting of ~24-30 plants per RIL per treatment, 1 µg of total RNA was converted to complementary DNA libraries (Molecular and Cellular Imaging Center, Ohio Agricultural Research and Development Center) and sequenced using Illumina HiSeq paired-end (PE) 100 base-pair reads.The average reads per sample was 3,126,863, ranging from 2,484,193 to 5,087,094 reads per sample.While only a single sample was sequenced for each RIL:treatment combination, the use of a highly inbred RIL population makes it such that each allele is represented approximately 46 times per treatment (50% of 92) and therefore provides replication within the genome-wide analyses.In brief, all reads were quality checked, adapter removed, and quality trimmed using bioinformatics tools FastQC (Andrews, 2010) and BBTools (BBMap;Bushnell, 2014).All reads were mapped to the reference genome (Glycine max Wm82.a2.v1) retrieved from Phytozome 1 , using CLCBio (CLC Genomics Workbench 9.5.3 2 .Reads mapped to the reference genome were counted using FeatureCounts (Liao et al., 2014).Data were normalized using edgeR (Robinson et al., 2010) using the trimmed mean of the Mvalues (TMM) normalization.

iBMQ eQTL mapping
The integrated hierarchical Bayesian model for multivariate eQTL mapping (iBMQ; Imholte et al., 2013) is a multivariate mapping method that implements a multi-loci approach to allow for the identification of complex traits that are being controlled by multiple genes.iBMQ parameters are estimated using a Markov Chain Monte Carlo (MCMC) algorithm that models concurrently all genes and SNPs.The current version of iBMQ executed is available online 3 .Normalized unfiltered expression data (RNAseq data) and genotyping data, including SNPs that were identified in the newly conducted linkage map for all RILs, were used as inputs into iBMQ.These data sets were formatted into an ExpressionSet and Snpset using Biobase in Bioconductor 3.7 (Huber et al., 2015).The eQTLs were mapped by using gene expression (etrait) as a phenotypic trait identifying the SNP association for mock and inoculated treatments and mapped separately with 100,000 iterations with a burn-in of 50,000 to produce a posterior probability of association (PPA).A false discovery rate (FDR) of 10% was used to determine the PPA significance cutoff (Imholte et al., 2013).
Significantly mapped eQTLs were classified as either cis-or trans-using the eqtlClassifier in the iBMQ package.Classification of cis-or trans-eQTLs was determined by physical position; if a SNP was within 5 Mb of the gene it controlled, it was classified as cis-, and if a SNP was farther than 5Mb from the gene controlled, it was classified as trans-.The physical positions of genes and SNPs were determined using the reference genome annotation, Williams 82 v.a2.v1 4 .
The function hotspot finder inside the iBMQ package was used to identify individual markers (SNPs) that are associated with the expression of several genes and indicate eQTL, hotspots as described by Imholte et al. (2013).A SNP was identified as a hotspot if it was associated with >20 genes.

Gene enrichment
Significant cis-and trans-eQTLs, gene suites at eQTL hotspots, and co-expression modules were subjected to GO enrichment.AgriGO (v2.0;Tian et al., 2017) was used for GO enrichment analysis.Targeted QDR pathway mechanisms were also explored by extracting genes via KEGG (Kanehisa and Goto, 2000 5 ) from NCBI BioSystems 6 .Gene models from the reference soybean genome that met read-mapping thresholds for inclusion in eQTL mapping were used as the reference set for enrichment analyses (Schmutz et al., 2010).Fisher's exact test (Fisher, 1935) and multitest Hochberg FDR adjustment (Benjamini and Hochberg, 1995) were used to determine significance.

Weighted gene co-expression analysis
WGCNA was carried out according to Langfelder and Horvath (2008) using modified R-scripts (WGCNA Methods M1).Prior to WGCNA analysis, RNA-sequencing counts from 92 RILs plus the parents of the population were normalized using R Bioconductor package edgeR (Robinson et al., 2010).First trimmed means of Mvalues (TMM) were used to calculate normalization factors for all sample counts and were subjected to a log2 transformation.Genes were filtered by the threshold of greater than 2 counts per million (cpm) across 50 samples.Ultimately, 27,666 genes were used in the input matrix.Due to sample count bias, a consensus approach was taken.Based on the smallest standard deviation representing susceptible and resistant individuals separately, 30 resistant and 30 susceptible individuals were selected, and the minimum standard deviation was determined (10,000 iterations) and used as a cutoff to select the subsample of 30 resistant and 30 susceptible individuals.
The process of subsampling individuals who met the standard deviation cutoff through generating modules was done in 20,000 iterations and generated for 20 network analyses.
To examine the physiological relevance of each module within the network analyses, phenotypic traits were correlated to the module genes' expression using log expression eigengene values for each module regressed against the lesion length values of the RILs.Consensus networks were constructed independently for the positively and negatively correlated modules and again related to phenotypic traits (WGCNA Methods M2).
To determine the driving factors of co-expression modules, modules that were significantly correlated to PRR disease (p-value ≤ 0.05) were also correlated to inoculated eQTL hotspots, adapting methods from Zhang and Horvath (2005) to integrate SNP-based significance (PPA of eQTL) with network properties (gene significance) of co-expression module (WGCNA Methods M3).

Co-localization of eQTLs with phQTLs
Co-localization of cis-eQTLs with phQTLs was based on the physical coordinates of the cis-eQTLs gene model being within the flanking marker positions for selected phQTLs reported from the present and previous studies (Table 1).All phQTLs that could not be positioned onto the genome based on marker positions were removed from the data set, for a total of six out of 28 phQTLs.Based on the physical coordinates of markers flanking phQTLs, overlapping phQTLs were consolidated.To determine if the number of co-localized cis-eQTLs and consolidated phQTLs was significantly greater than expected, cis-eQTLs were permuted across all gene models eligible for eQTL analysis.A threshold (a = 0.05) was set based on 1,000 permutations.
Methods to determine the co-localization of trans-eQTLs with phQTLs were adapted from Christie et al. (2017).Base pair positions of consolidated phQTLs and a window using the average linkage disequilibrium (LD) block size (242 kb) around the marker to which trans-eQTLs mapped were used to determine phQTL/trans-eQTL co-localization.Linkage disequilibrium blocks were determined using the Haploview four-gamete method (Barrett et al., 2005).To determine if the number of co-localized trans-eQTLs and phQTLs was significantly greater than expected, trans-eQTLs were permuted across available markers on the genetic map.A threshold (alpha = 0.05) was set based on 1,000 permutations.

Identification of putative master regulators
Methods for identifying master regulators for eQTL hotspots were adapted from Wang et al. (2017).Genes within an up-and downstream 242-kb window (average LD block size for this population) of the hotspot SNP that also had cis-eQTLs mapping were identified as initial candidate master regulators.If no genes within the window were associated with cis-eQTLs, genes with no mapping etraits were considered.From these cis-eQTLs or positional candidates, putative master regulators were selected according to predicted functions of transcription factor (TF) or signaling molecules (SM).

Genetic map reconstruction and quantitative disease resistance to Phytophthora sojae
While phQTLs for QDR have been previously mapped in the F 9:11 RIL population derived from a cross between Conrad (resistant) and Sloan (susceptible) (Stasko et al., 2016), in order to ensure that the phQTLs were directly relevant to the expression data collected in this study (Supplementary Figure S2), we mapped phQTLs using data from only the subset of 92 RILs for which RNA-seq, phenotypic data, and previous genotypic data were available (Stasko et al., 2016).A genetic map was constructed consisting of 1,122 markers assembled in 28 linkage groups, with most of the 20 soybean chromosomes represented as a single linkage group and chromosomes 5, 6, 7, 11, 12, 13, 17, and 19 each represented by two linkage groups (Supplementary Table S1).A suggestive phQTL, significant at the chromosome but not genome-wide level, on chromosome 18 and nonsignificant associations with regional peaks but nonsignificant logarithm of the odds (LOD) scores on chromosomes 1, 16, and 19a correspond to phQTLs previously identified in the full RIL population (Supplementary Figure S2; Stasko et al., 2016).In comparison to previous studies, reduced significance is expected due to the smaller subset of RILs.

Genetic architecture of gene expression in mock and inoculated treatments
To identify the loci contributing to variation in gene expression, eQTLs were mapped under both mock and inoculated conditions.Gene expression levels were interpreted as quantitative traits (e-traits), and the locus or loci for each e-trait were mapped separately for inoculated and mock treatments.A total of 114,197 eQTLs were mapped from e-traits for the inoculated treatment from root samples collected 24 hai from the site of inoculation, representing transcripts from 35,781 unique genes associated with 74 unique loci (Supplementary Table S2).In contrast, the number of eQTLs mapped in the mock treatment was far lower but distributed across more loci, with 794 eQTLs representing 788 unique genes across 234 unique loci.The eQTLs identified in the inoculated treatment were distributed across most chromosomes, the exceptions being chromosomes 9 and 10 (Figure 1A), while eQTLs from the mock treatments were relatively evenly distributed across all chromosomes in the genome (Figure 1B).Interestingly, only two eQTLs were in common between mock and inoculated treatments, including the SNP, ss715580997, associated with e-trait Glyma.02G107800(putatively coding an uncharacterized protein), and the SNP O S U _ S N P _ G l y m a1 9 g 4 12 10 , as s o c i a t e d w i t h e -t r a i t , Glyma.19G224300 (putatively coding a protein involved in regulation of root development and cell wall).
The eQTLs associated with e-traits of nearby genes are considered cis-eQTLs, while eQTLs altering the expression of physically distant genes are classified as trans-eQTLs (Hubner et al., 2005;Kliebenstein, 2009).In the inoculated treatment, the vast majority (111,568; ~98%) of the eQTLs were in the trans configuration (Figure 1A; Supplementary Table S2).In contrast, in the mock treatment, 83% (657) of the eQTLs were in the cis configuration (Figure 1B; Supplementary Table S2).
Key regulatory genes are expected to influence the expression of many genes.Therefore, eQTL hotspots, genomic regions enriched for eQTLs, are likely localized to these key regulatory genes (Kliebenstein, 2009;Tian et al., 2016).Thirty-six hotspots distributed across 12 chromosomes were identified for the inoculated treatment, with seven located on chromosome 17.Of the 36 hotspots in the inoculated treatment, there were two very large hotspots, GM_1 and GM_15, associated with 23,919 and 25,727 e-traits, respectively.For the mock treatment, there were three hotspots located on chromosomes 13, 16, and 18.The mock and inoculated eQTL hotspots did not overlap with each other (Table 2).

Term enrichment and functional annotations of genes associated with each hotspot
Genes that are co-regulated may be associated with specific coordinated functions or resistance mechanisms.To identify a commonality of functions among e-traits mapping to each eQTL hotspot, GO term enrichment analysis was done.The e-traits associated with 10 of the 36 hotspots for the inoculated treatment were significantly enriched (p-value ≤ 0.05) for biological and/or cellular processes and molecular function terms.We were particularly interested in the enrichment of biological process GO terms, as enrichment of these terms may provide clues to the mechanisms of QDR.The biological process GO terms were significant for the e-traits associated with six of the 36 hotspots (Supplementary Table S3).A total of 16 and 31 biological processes were enriched among e-traits mapping to the large hotspots GM_1 (23,919 e-traits) and GM_15 (25,727 e-traits), respectively.Some of the most significantly enriched biological process terms were "intracellular transport," "gene expression," and "cellular processes" (Supplementary Table S3).Similar terms were also found for GM_2_B, GM_11, and GM13 as well as "intracellular signal transduction".Four additional hotspots (GM_17_D, GM_17_F, GM_18, and GM_18_A) were not enriched for biological process GO terms but were enriched for cellular and molecular function GO terms.Overall, biological process GO term enrichment provided evidence of a functional relationship among etraits within 10 hotspots and hinted at cell-to-cell signaling and protein modification roles but did not elucidate their involvement in a specific mechanism or pathway for QDR.GO term enrichment for the e-traits associated with the three mock eQTL hotspots was also identified but revealed very general terms (Supplementary Table S4).While e-traits associated with each of the three mock eQTL hotspots were enriched for GO terms, only e-traits associated with GM_18_M were enriched for biological process terms.These included the terms "regulation of cellular process" and "biological regulation."Any significant enrichment of GO terms indicates functional relationships among the genes at a hotspot; however, these general terms do little to inform the specific role of the e-traits at these hotspots.
The e-traits associated with each eQTL hotspot were also examined for enrichment of genes predicted to be involved in specific mechanisms of QDR in soybeans by P. sojae.The e-traits for six of 36 hotspots from the inoculated treatment (GM_1, GM_1_B, GM_2_B, GM_11, GM_13, and GM_15) were significantly enriched for genes in the plant-pathogen interaction pathway, and one hotspot (GM_11) was additionally significantly enriched for genes in the isoflavonoid pathway (Table 3).These findings suggest that these six hotspots are most likely regulating genes involved in PAMP-triggered immunity, defense-related gene induction, and/or programmed cell death, implicating leucine-rich repeat (LRR) encoding and other PAMP-triggering genes as likely candidate genes.
We identified no enrichment for hypothesized resistance mechanisms for e-traits mapping to the three hotspots from the mock treatment nor for the remaining 30 of 36 hotspots from the inoculated treatment.While these findings may indicate a lack of concerted functional relationships among e-traits mapping to these hotspots, it may also be due to a lack of statistical power in the smaller hotspots or that the associated e-traits may be involved in unknown or untested mechanisms of QDR.

Weighted gene co-expression network analysis in the Conrad × Sloan RIL population
To confirm and further characterize the transcriptional reprogramming that occurs following infection with P. sojae of  Fisher, 1935;Benjamini and Hochberg, 1995).
this RIL population, we constructed expression networks, or weighted gene co-expression networks, determined by pairwise correlation of gene expression profiles (Langfelder and Horvath, 2007).This network analysis resulted in a total of six robustly defined modules (Supplementary Figures S4-S6).Module eigengene expression values from three modules, honeydew1, dark red, and grey, comprised of 20,976, 2,785, and 385 genes, respectively, had significant positive correlations with both BLUP values and lesion length for the PRR disease phenotype (Figure 2), indicating that as the expression of genes within these modules increased, susceptibility to PRR increased.

Term enrichment and functional annotations of co-expression modules
Following similar methods as to how we assessed the putative functions of eQTL hotspots, to ascertain the potential roles of the modules in susceptibility to PRR, GO term enrichment and functional annotations were assessed for genes within the honeydew1, dark red, and grey modules.While no significant GO enrichment was observed for the smaller dark red and grey modules, for the honeydew1 module, the most significantly enriched GO terms for biological processes included "translation," "protein localization," and "macromolecule localization" (Supplementary Table S5).
Genes within the co-expression modules were also evaluated for gene enrichment within the six known QDR pathways, for which we had previously evaluated e-traits in each hotspot for enrichment.The honeydew1 module was enriched for five of the six evaluated pathways: plant hormone signal transduction, phenylpropanoid biosynthesis, isoflavonoid biosynthesis, MAPK signaling, and plantpathogen interaction pathways (Table 4).The dark red module was also significantly enriched for plant hormone signal transduction, MAPK signaling, and plant-pathogen interaction pathways.None of the six pathways were significantly enriched in the grey module, possibly due to the relatively smaller size of this module or indicating that there may be other pathways associated with resistance.

Genetic architecture of co-expression modules
We integrated genetic markers associated with inoculated eQTL hotspots and co-expression modules (Zhang and Horvath, 2005) to further investigate the genetic loci influencing the three coexpression modules that were significantly correlated with PRR (honeydew1, dark red, and grey).The expression of each of the three module eigengenes was significantly correlated with the SNP (s) marking one or more of the eQTL hotspots.The correlation of an eQTL hotspot and a co-expression module to the same SNP suggests regulation of both by a common genetic mechanism (Table 2; Supplementary Figures S7A-C).The grey module was significantly correlated to the genetic marker for hotspot GM_5 (423 e-traits).The expression of both the honeydew1 and dark red module eigengenes was significantly correlated to SNPs from hotspots GM_13_A (20 e-traits) and GM_17_D (303 e-traits).

Co-localization of phQTLs and eQTLs
In order to understand the regulation of QDR by P. sojae, colocalizations of eQTLs or eQTL hotspots with phQTLs were identified.Among all of the 8 phQTLs mapped to unique locations in populations derived from Conrad × Sloan for P. sojae isolate 1.S.1.1 (Wang et al., 2010;Wang et al., 2012a;Stasko et al., 2016; Table 1), three phQTLs, which mapped to chromosomes 1, 18, and 19, co-localized with three hotspots and nine eQTLs (three inoculated trans-and six inoculated cis-eQTLs (Table 1).For both cis-and trans-eQTLs, this is significantly more co-localization than expected by random chance.
Co-localized with phQTL_1 were four cis-eQTLs for Glyma.01G160600,Glyma.01G162600,Glyma.01G170600, and Glyma.01G171300 and two trans-eQTLs, both representing Glyma.10G026500, which had e-traits mapping independently to two markers within this phQTL region (Table 1).The e-traits of genes mapping in cis are towards genes annotated as a vacuolar iron transporter homolog, two-component response regulator ARR2, and an uncharacterized expressed sequence.The predicted protein of Glyma.10G026500,mapping to the two trans-eQTLs, was also uncharacterized.Interestingly, just outside (138 kb) of the average LD block size window (242 kb), which we used to consider co-localization between mapped trans-eQTL and phQTL, is the hotspot GM_1_B.
The e-traits for the paralogs Glyma.18G270900 and Glyma.08g249200 are both co-localized with phQTL_18b in cis and trans, respectively (Table 1).The gene sequences are present in syntenic blocks from the recent duplication within the soybean genome (Schmutz et al., 2010).Both genes are predicted to encode a malectin/receptor-like protein kinase.
Covering 4.7 Mb, phQTL_19b represented the largest genomic size of the phQTL we considered.As such, it co-localized with all three hotspots on chromosome 19 (GM_19, GM_19_A, and GM_19_B).However, these hotspots, all relatively small and possessing between 30 and 67 e-traits, were not enriched for any function or GO annotation, providing little evidence of a coordinated function.However, in addition to co-localization with the three hotspots, phQTL_19b also co-localized with both the mock and inoculated cis-eQTLs for the Glyma.19G224300e-traits.

Identification of master regulators controlling the expression of downstream suites of genes
Master regulators control the expression of a suite of downstream genes.To identify master regulators in QDR for P. sojae, we examined the key eQTL hotspots, which we defined as those that were significantly enriched for e-traits, correlated with modules, or colocalized with phQTL.From the 36 total eQTL hotspots, our primary and meta-analyses identified 16 key eQTL hotspots through the significant pathway or GO term enrichment of e-traits within a hotspot (GM_1, GM_1_B, GM_2_B, GM_11, GM_13, GM_15, GM_17, GM_17_D, GM_17_F, GM_18, GM_18_A), by the significant correlation of co-expression modules (GM_5, GM_13_A, GM_17_D), and/or the co-localization with a phQTL for resistance to P. sojae (GM_19, GM_19_A, GM_19_B).To better understand the regulation of these key hotspots, we identified their putative master regulators by integrating genetic position and putative gene function (methods adapted from Wang et al., 2017).We identified genes within the hotspot regions predicted to encode transcription factors or signaling molecules as putative master regulators (Table 5).
For nine of the 16 hotspots evaluated for master regulators (GM_1, GM_1_B, GM_2_B, GM_5, GM_11, GM_13, GM_13_A, GM_15, GM_17), we identified co-localized cis-eQTLs putatively encoding signal molecules and/or transcription factors (Table 5).To note, numerous hotspots were associated with multiple putative master regulators with genes for cis-eQTLs encoding both signaling molecules and transcription factors.For the four remaining eQTL hotspots (GM_17_D, MD_17_F, GM_18, and GM_18_A), while no genes  within the LD window (242 kb) had cis-eQTLs, genes within the LD windows were putatively encoding transcription factors and/or signaling molecules and identified as candidate master regulators of these hotspots.In total, 15 genes putatively encoding transcription factors and 24 genes  The hotspot SNP is located within this putative master regulator.
putatively encoding signaling molecules were identified as candidate master regulators for the 16 key eQTL hotspots (Table 5).

Discussion
Many phQTLs have been mapped through several generations in this Conrad × Sloan population (Wang et al., 2010;Wang et al., 2012a;Wang et al., 2012b;Stasko et al., 2016); however, identifying the mechanisms underpinning these phQTLs has, thus far, been unsuccessful.Studies using the resistant parent 'Conrad' have identified many putative mechanisms of quantitative resistance, demonstrating the complex nature of QDR and the potential that these mechanisms could be interacting (Vega-Sańchez et al., 2005;Thomas et al., 2007;Ranathunge et al., 2008;Zhou et al., 2009;Wang et al., 2012b).The application of a systems genomics approach in this study has allowed us to disentangle the complex genetic architecture of gene expression related to QDR for P. sojae using multiple approaches, including eQTL mapping, co-expression network analysis, and the colocalizations of phQTLs, eQTLs, and co-expression modules.The approaches taken in this study confirmed hypothesized mechanisms as well as provided evidence to suggest potential novel mechanisms of QDR for this pathosystem.
4.1 Inoculation with P. sojae causes transcriptional reprogramming that occurs in a trans-regulatory manner Expression QTLs were successfully mapped in both the inoculated and mock treatments; however, there were 144-fold more eQTLs mapped for the inoculated treatment compared to the mock treatment.An average of two eQTLs per gene were mapped in the inoculated treatment, which is consistent with previous eQTL mapping studies (Schadt et al., 2003;Swanson-Wagner et al., 2009;Christie et al., 2017).Yet, in this study, only 794 eQTLs were identified for the mock, the majority of which were cis-eQTL.In stark contrast to the mock-inoculated treatment, nearly all eQTLs identified from the inoculated treatment in this study were trans-eQTLs (98%).
Gene expression has been shown in previous eQTL studies to be controlled by trans-or a combination of both trans-and ciselements (West et al., 2007;Christie et al., 2017;Sun et al., 2017;Li et al., 2018).However, the number of eQTLs mapped, as well as the proportion of trans-vs.cis-eQTLs, has not followed a specific trend between species and populations (Keurentjes et al., 2007;West et al., 2007;Potokina et al., 2008;Swanson-Wagner et al., 2009;Hammond et al., 2011;Christie et al., 2017).For example, the number of eQTLs varied approximately ninefold between RIL populations in Arabidopsis (Keurentjes et al., 2007;West et al., 2007).While some studies have reported nearly equal ratios of trans-and cis-eQTLs detected in both Arabidopsis and barley (Keurentjes et al., 2007;Potokina et al., 2008), other studies have revealed a predominance of trans-regulation of eQTLs in Arabidopsis (86% and greater trans-eQTLs) (West et al., 2007;Soltis et al., 2020), barley (70% trans-eQTLs) (Druka et al., 2008), Brassica rapa (77% trans-eQTLs) (Hammond et al., 2011), and maize (up to 80% trans-eQTLs) (Swanson-Wagner et al., 2009;Christie et al., 2017).These varying results have been attributed to statistical power to detect trans-eQTLs, the size of the mapping population, the high polymorphism rate among genotypes in the study, and true biological differences between systems and their overall genetic architecture (Kliebenstein, 2009;Soltis et al., 2020).To date, only a few studies have been completed to address these questions of differing detection of eQTL types across plant species and populations (Franceschini et al., 2012;Saha and Battle, 2018).
The majority of plant-based eQTL studies to date have focused on natural genetic variation within breeding populations, different stages of maturation, or specific production or accumulation of compounds, and few have focused on mechanisms of disease resistance.Specifically in soybean, previous eQTL analyses identified predominantly trans-eQTLs for the genetic architecture of immature soybean seed (86.6%, trans-eQTLs) (Bolon et al., 2014) and dissection of isoflavonoid accumulation in soybean seed (60.6%, trans-eQTLs) (Wang et al., 2014).
The large number of trans-eQTLs mapped in this study were primarily associated with only eight eQTL hotspots, indicating massive transcriptional reprogramming resulting from the inoculation of soybean with P. sojae.This confirms several previous studies for quantitative resistance (Zhou et al., 2009;Soltis et al., 2020) and Rpsgene-related responses (Lin et al., 2014;Hale et al., 2023b).The eQTL hotspots identified in this study at 24 hai may represent key regulatory hubs and the control of signaling networks specifically in response to infection by P. sojae.More importantly, none of the 36 hotspots identified from the inoculated treatment overlap with the three hotspots identified from the mock treatment.This suggests that these hotspots in the mock represent either constitutive differences in regulation between Conrad and Sloan or a response specific to the mock treatment at the 24-hai time point.Of the few eQTL studies that have mapped transcriptional responses to disease, only a fraction of these studies compared eQTLs mapped in disease versus non-disease conditions in plant systems.Moscou et al. (2011), using microarrays to assay transcripts in barley following both inoculations with Puccinia graminis and mock inoculation, had findings that differed from this study, with similar numbers of eQTLs mapped in both mock and inoculated samples and the majority classified as cis.Here, the differences in the number of eQTL mapped between treatments, the lack of concordance of the hotspots between the mock and inoculated treatments, the correspondence with phQTLs, and the functional enrichment of genes within hotspots together suggest that the changes in transcription are due to infection by P. sojae through a coordinated transcriptional response of multiple plant defense mechanisms.

Major eQTL hotspots and coexpression networks elucidated potential QDR mechanisms, including signal integration and defense action via cell wall strengthening
Expression QTL hotspots are a single polymorphism associated with the expression of numerous genes (Neto et al., 2012), and the genetic regions may harbor important regulatory genes.In this study, the four largest hotspots, mapping to chromosomes 1, 2, 11, and 15, accounting for more than 80,000 eQTLs (>70%) from the inoculated treatment, were each enriched for genes within the plant-pathogen interaction (PPI) pathways.Specifically, PPI pathway genes found within the hotspots were predicted to function throughout the pattern-triggered immunity (PTI) pathway.The separation of the PTI and ETI pathways within the context of plant resistance to oomycetes has recently come into question in favor of a three-layer plant immune system (consisting of the recognition, signal integration, and defense-action layers) describing both PTI and ETI for plant-pathogenic oomycete infection (Wang et al., 2019;Naveed et al., 2020).The signalintegration layer represents a complex network of pathway cascades including phosphorylation, ubiquitination, relocation, degradation, stabilization of proteins, transcriptional regulation, and chemical signaling (Wang et al., 2019).The QDR pathway enrichment within these major hotspots, along with the significantly enriched GO terms related to cell-to-cell signaling and protein modification, support the involvement of these hotspots in the signalintegration layer of defense.This aligns with previously implicated defense mechanisms in plant-oomycete interactions and specifically within the P. sojae-soybean pathosystem (Wang et al., 2012b;Wang et al., 2019).
These four eQTL hotspots represent genetic variation for transcriptional reprogramming resulting from inoculation with P. sojae, a phenomenon that has been previously reported (Zhou et al., 2009;Wang et al., 2010).Yet, these hotspots do not localize to the regions of any phQTLs identified in this study or previous studies.Samad-Zamini et al. (2017) had similar findings, where none of the hotspots co-localized with phQTLs during a time-course assay of Fusarium graminearum (FHB) infection of wheat (Triticum aestivum L.).This lack of co-localization may be due to residual genetic variation of e-traits not significantly attributed to phQTL; this variation may involve complex genetic interactions, including epistasis (Li, 2019).
In this study, we also identified a total of 24,146 genes within three co-expression modules that were significantly correlated to the PRR disease resistance response.The SNPs corresponding to two hotspots, GM_13_A and GM_17_D, were also both significantly correlated to the co-expression modules dark red and honeydew1.GM_17_D was enriched for GO terms including "hydrolase activity" and "cell wall structure," functions that align with the third layer of resistance, defense-action (Wang et al., 2019).Genes involved in the modification of the cell wall and hydrolase activity have been shown to be involved in plant defense responses (Smith et al., 1988;Smith et al., 1990;Walton, 1994;Minic, 2008).Hydrolase expression has been previously associated with quantitative resistance to P. sojae with specific hydrolases suppressed at 48 hai in the resistant parent Conrad and suppressed at 72 hai in the susceptible parent Sloan (Wang et al., 2012b).The GM_17_D hotspot may be involved in the coordinated regulation of these two modules, and importantly, these coexpressed genes may play a role in limiting pathogen penetration into the cell wall in response to pathogen presence.The dark red and honeydew1 co-expression modules, along with the third co-expression module, grey, were also each enriched for genes involved in the plant hormone signal transduction, MAPK signaling pathway, and PPI pathways, providing evidence for coordinated regulation of expression among each of these defense mechanisms.
The honeydew1 module was further described by its additional enrichment of genes from the phenylpropanoid biosynthesis and isoflavonoid pathways.The role of phenylpropanoid and isoflavonoid in the R-gene-mediated response has been well studied (Graham et al., 2007), including the recent identification of a transcription factor that modulates this response (Jahan et al., 2020).However, the specific role of phenylpropanoid and isoflavonoid pathways in QDR has been more elusive.Gene expression in the honeydew1 module is correlated with increased susceptibility, supporting recent evidence in the cross-talk that occurs between the pathways for R-gene-mediated and QDR.Previous studies showed both SA and JA increasing at inoculated roots, with JA further increased in later time points after inoculation (Stasko et al., 2020;Karhoff et al., 2022).This 24-hai time point could be a critical time as the pathogen switches from hemibiotrophy to the necrotrophic phase (Moy et al., 2004).Several genes in the phenylpropanoid pathway have been identified as playing a role in resistance to P. sojae in soybeans.For example, soybean cinnamate 4-hydrolase (GmC4H1; first hydroxylation step of the phenylpropanoid pathway) was induced at 24 hai in the resistant parent Conrad, and greater colonization of P. sojae was measured in GmC4H1-silenced plants (Yan et al., 2019).Recently, in the wheat-Fusarium graminearum system, it was reported that wheat genotypes with greater levels of resistance had a constitutive expression of genes for plant cell wall biogenesis and terpene biosynthesis (Buerstmayr et al., 2021).
4.3 Genetic regulation of trans-eQTLs in response to the pathogen is predicted to occur through TF and signaling molecules involved in PPI, plant hormone signal transduction, and novel mechanisms of resistance Numerous studies have proposed that cis-acting mechanisms (i.e., transcription factors) can affect the expression of e-traits at trans-eQTL hotspots (Albert and Kruglyak, 2015;Wang et al., 2017).Thus, cis-eQTLs located near regulatory genes have the potential to be master regulators for these e-traits in a given hotspot (Bryois et al., 2014;Albert and Kruglyak, 2015;Yao et al., 2015).We identified candidate master regulators for the key hotspots that had significant GO or pathway enrichment, were correlated to a co-expression module, and/or were co-localized with phQTL.Of these, four putative master regulators were predicted to function in the PPI pathway as LRR-RLKs, or EFhand motif proteins.
Two EF-hand motif proteins were each identified as candidate master regulators associated with the PPI pathway for hotspots GM_2_B and GM_15, respectively, with GM_2_B being a hotspot enriched for genes within the PPI pathway.Approximately 250 EFhand motifs have been identified in plants and are involved with Ca 2+ , which acts as a messenger that regulates responses to external stimuli, development, and hormones, including plant defense and stress response (Poovaiah and Reddy, 1993;Trewavas and Mahlo, 1998;Reddy and Reddy, 2001;Zielinski, 1998).The majority of Ca 2+ sensors in soybeans possess the EF-hand motif and have at least one or more hormone-or stress-response-related cis-elements in their promoter region (Zeng et al., 2017).These hormone-or stress-response-related elements have been characterized by functioning in the regulation of abscisic acid (ABA) signaling, auxin response, ethylene response, and phosphate starvation response.Of these signaling and responses potentially regulated by EF-hand motif encoding genes, ABA has been shown to be a negative regulator of R-gene-mediated resistance (Ward et al., 1989;MacDonald and Cahill, 1999), and auxin has been reported to enhance plant susceptibility to P. sojae in soybean (Stasko et al., 2020) and other pathogens (Wang et al., 2007;Domingo et al., 2009;Kidd et al., 2011).Auxin transporters and auxin-induced proteins have been upregulated in susceptible parents in the P. sojae-soybean pathosystem (Wang et al., 2012b).Auxin transport transcripts of GmPIN were higher in expression in the resistant Conrad following inoculation with P. sojae compared to mock, whereas in the susceptible, fewer GmPIN changed in expression levels (Stasko et al., 2020).Additionally, ethylene-responsive genes have also been known to induce resistance in the P. sojae-soybean pathosystem (Sugano et al., 2013;Zhao et al., 2017), as well as play a role in the regulation of pathogenesis-related gene expression (Lorenzo et al., 2003;Pieterse et al., 2009;Rehman and Mahmood, 2015).Taken together, the EF-hand motif-encoding genes are excellent candidate master regulators for GM_2_B and GM_15.
In addition to these candidate master regulators within the PPI pathway and those that overlap phQTLs, several MYB-TFs were identified as candidate master regulators for GM_11 and GM_19, GM_19_A, and GM_19_B.MYB transcription factors are one of the six major TF families functioning in plant defense (Ng et al., 2018), responding to both abiotic and biotic stresses, and functioning in primary and secondary metabolism (Stracke, 2001;Ambawat et al., 2013), including the regulation of the phenylpropanoid pathway (Liu et al., 2015).Additionally, GmMYB29A2 is essential for the Rgene response to P. sojae in soybeans, regulating the accumulation of glyceollin in Williams 82 (Jahan et al., 2020), and MYB transcripts were also detected by capture-seq from a transcriptome data set of the R-gene response in Williams 82 (Hale et al., 2023b).They are known to act as a positive regulator of hypersensitive response in PCD in response to fungal and bacterial pathogens (Vailleau et al., 2002).Thus, these may be putative master regulators involved with the positive regulation of PCD in the soybean-P.sojae pathosystem.
GM_17_D did not have any cis-eQTLs mapping to this hotspot.However, significant differential expression is not a requirement for a master regulator.For example, hunchback (hb), encoding a ZN-finger TF in Drosophila melanogaster, was identified as a candidate master regulator for mitigation of lead exposure, located near a trans-eQTL hotspot, yet the candidate itself had no e-traits mapped (Qu et al., 2018).Here, we note that the zinc finger TF (Glyma.17G200200) is a candidate master regulator of GM_17_D because it is not only physically located near the hotspot but is also within the honeydew1 co-expression network that is correlated with PRR disease.
In addition to the candidate master regulators functioning in known or hypothesized pathways for QDR, we also identified candidate master regulators that putatively influence novel pathways for QDR.These novel pathways for QDR in P. sojae included secondary metabolite biosynthesis, RNA transport, thioredoxin metabolism (GM_1_B), lysine degradation (Glyma_2_B), reticulon (Lee et al., 2011), starch and sucrose metabolism (GM_13), thymine degradation (GM_15), and a number of serine/threonine protein kinases and phosphatases that impact other metabolic pathways.These candidate master regulators of novel QDR pathways include TFs and signaling molecules that potentially regulate the expression of downstream genes related to hotspots.Further studies will be needed to determine if and how these pathways are playing a role in the P. sojaesoybean pathosystem.

Co-localization of phQTLs with eQTL points to causal candidate genes for QDR
To identify gene expression variation that may be causal to PRR disease resistance, we focused on those eQTLs that co-localized with phQTLs, indicating a strong link between transcriptional phenotype and the genes underpinning the disease resistance phenotype.Specifically, co-localized cis-eQTLs, the genes regulating colocalized trans-eQTLs, or trans-eQTL hotspots may be causal for resistance to PRR.
While none of the four cis-eQTLs co-localized with phQTL_1 have obvious functions in quantitative disease resistance, the hotspot GM_1_B, which neighbors phQTL_1, is enriched for genes functioning in PPI pathways, making regulators of this hotspot viable causal genes for phQTL_1.Three cis-eQTLs were identified as candidate master regulators for this hotspot: Glyma.01G156600,Glyma.01G156200, and Glyma.01G156700.These genes are predicted to encode a thioredoxin reductase, a membrane transport protein, and a hydroxymethylglutaryl-CoA reductase, respectively.The predicted thioredoxin reductase (Glyma.01G156600) is of interest given the role of thioredoxin in disease resistance and potentially QDR to P. sojae, with thioredoxin-encoding genes identified as candidate genes for several quantitative disease resistance loci (QDRL) towards P. sojae (Huang et al., 2016;Stasko et al., 2016).Additionally, a thioredoxin-encoding gene has been shown to be the causal gene for resistance at the Scmv1 phQTL for sugarcane mosaic virus in maize (Liu et al., 2017).The predicted membrane transport protein encoded by Glyma.01G156200 is an auxin efflux carrier.Auxin has been previously described in numerous studies as being involved in susceptibility to plant pathogens (Wang et al., 2007;Domingo et al., 2009;Kidd et al., 2011;Pieterse et al., 2012).Finally, hydroxymethylglutaryl-CoA reductase, predicted to be encoded by Glyma.01G156700, is involved in terpenoid and secondary metabolite biosynthesis (Antolin-Llovera et al., 2011).Thus, these potential master regulators for the Gm_1_B hotspot may represent the causal genetic variation for phQTL_1.
The phQTL_18b was co-localized with two eQTLs controlling the expression of Glyma.18G270900 and Glyma.08g249200 in cis and trans, respectively, each putatively encoding a malectin/ receptor-like protein kinase.In Arabidopsis, the homolog of these genes, FERONIA, has been experimentally shown to have multiple functions, including as a modulator of ethylene response (Deslauriers and Larsen, 2010) and in reactive oxygen species (ROS)-mediated root hair development (Duan et al., 2010).ROS is a well-known mediator of stress-induced responses and functions in growth and development (Werner, 2004;Swanson and Gilroy, 2010;Torres, 2010).FERONIA also functions to inhibit jasmonic acid (JA) signaling through phosphorylation of the transcription factor MYC2 in Arabidopsis (Guo et al., 2018).In soybean, a role for JA was proposed in the later stages of infection by P. sojae (Stasko et al., 2020).The JA pathway was suppressed in incompatible Rgene reactions to P. sojae (Lin et al., 2014), and JA accumulation significantly increased in P. sojae-inoculated susceptible lines in contrast to the mock-inoculated and to lines with quantitative resistance alleles (Karhoff et al., 2022).
In addition to co-localization with the three hotspots (GM_19, GM_19_A, and GM_19_B), phQTL_19b co-localized with both the mock and inoculated cis-eQTLs for the Glyma.19G224300e-traits.This inoculated cis-eQTL is also part of the GM_19_B eQTL hotspot.The etraits for Glyma.19G224300represented one of only two pairs of e-QTLs that were found under both mock and inoculated conditions, indicating possible constitutive control of both GM_19_B and of phQTL_19b.Glyma.19G224300 is predicted to encode a germin-like protein (GLP).Among their functions, GLPs can be involved in response to abiotic stress (Barman and Banerjee, 2015).In Arabidopsis, upregulation of the Glyma.19G224300homolog AT1G09560 results in reduced primary root and enhanced lateral root growth (Ham et al., 2012).Glyma.19G224300 may function in root architecture, providing constitutive quantitative resistance to P. sojae, with the differences in disease resulting from expression changes mapping to GM_19_B.

Concluding remarks
This vast transcriptional reprogramming due to pathogen infection compared to the nondisease state had not been previously explored through eQTL methodology using RNA-sequencing data in this hostpathogen system.Ultimately, this study identified gene co-expression modules associated with resistance and susceptibility to P. sojae in this RIL population.Clearly, the transcriptional response to this pathogen is complex, as there were more than 100-fold greater number of eQTLs in the inoculated compared to the mock treatment, as well as a predominance of trans-eQTLs in the inoculated over the mock treatment.Further evidence supporting cell wall structure, auxin response, jasmonic acid signaling, and PPI receptor and signaling genes as mechanisms of resistance are provided, as well as several new potential mechanisms for regulating resistance as well as potential susceptibility factors.Further confirmation of the candidate genes regulating trans-eQTLs and/or acting as the causal variation of phQTLs will need to be explored through functional studies.The development of this large dataset and analyses through co-expression networks, eQTLs, and phQTLs have the potential to be expanded to elucidate more biologically relevant information on P. sojae infection as well as constitutive differences between two cultivars.
FIGURE 1 Expression quantitative trait loci (eQTL) were mapped in 92 Conrad × Sloan F 9:11 RIL population at false discovery rate (FDR) of 10% with 100,000 iterations.(A) Significantly mapped inoculated treatment eQTL at a significance posterior probabilities of association (PPA) cutoff of 0.786 for the inoculated treatment.(B) Significantly mapped mock treatment eQTL at a significance PPA cutoff of 0.688.Base-pair positions are shown with zero at the bottom of the y-axis and right of the x-axis.Markers in a vertical formation indicate trans-eQTL, and markers in a diagonal formation indicate cis-eQTL.Chromosomes with no data contain nonclassified eQTL due to missing physical location data.

FIGURE 2
FIGURE 2Consensus-weighted gene co-expression network analysis (WGCNA) for positively correlated modules.The best linear unbiased predictors (BLUPs) and lesion length values positively correlated to co-expression modules are in red, and those negatively correlated are in green.The first number represents the correlation coefficient, the number in parenthesis indicates the p-value.N/A indicates no consensus co-expression module was preserved.

TABLE 1
Co-localization of cisand trans-eQTL as well as hotpots with consolidated phenotypic quantitative trait loci (phQTLs) mapped in the subpopulation and previous phQTL mapped in multiple generations of the Conrad × Sloan RIL population.

TABLE 2
Summary of inoculated and mock causal hotspots for number and regulation of expression quantitative trait loci (eQTLs) positively correlated to resistance towards Phytophthora sojae.

TABLE 3
Pathway enrichment of expression quantitative trait loci hotspots.

TABLE 4
Enrichment of genes within KEGG pathways for each co-expression module related to susceptibility.

TABLE 5
Candidate master regulators of key hotspots.
b Candidate master regulators were identified as genes within a 242-kb (average linkage disequilibrium (LD) block size) upstream or downstream of the hotspot SNP, and genes with cis-eQTL mapping were identified as initial candidate master regulators.Genes were then selected as candidate master regulators if their putative function included transcription factors, signaling molecules, or known associations involved in the enriched quantitative disease resistance (QDR) pathways.Gene ID based on the Wms82.a2.v1 sequence (soybase.org).cDescription, PFAM, Panther, and pathway were retrieved from https:// phytozome.jgi.doe.gov(accessedJuly 2023).df i Co-localizes with phenotypic quantitative trait loci (phQTL).jIntragenic hotspot SNP.k