Ten soybean genotypes were selected as a representation of the range of seed protein content observed in Canada. Soybean lines are listed from the lowest to highest seed protein content, with line 1 having the lowest average protein content and line 10 having the highest average protein content. The soybean lines used in this study were all developed at the Ottawa Research and Development Centre by Agriculture and Agri-Food Canada: X5583-1-041-5-5 (line 1), AC Harmony (Voldeng et al., 1996a) (line 2), AAC Halli (line 3), 90A01 (Cober et al., 2006) (line 4), Maple Amber (line 5), OT13-08 (line 6), OT14-03 (line 7), AAC Springfield (line 8), Jari (line 9), and AC Proteus (Voldeng et al., 1996b) (line 10).

Planting was performed in 2019 in replicated trials across four locations: Ottawa Ontario (latitude 45.39°, longitude −75.72°), Morden Manitoba (49.18°, −98.08°), Brandon Manitoba (49.86°, −99.98°), and Saskatoon Saskatchewan (52.15°, −106.57°). Seeds were planted in quadruplicate at the mid-end of May in a 4 × 5 rectangular lattice arrangement at a density of 50 seeds per m², and appropriate crop management practices were taken at each site. For additional information on planting, see Cober et al. (2023).

Young trifoliate leaf tissue was collected in triplicate from soybeans at the R5 (Pedersen and Licht, 2014) stage of maturity from otherwise healthy-looking plants. Tissue was flash-frozen in liquid nitrogen in the field immediately upon harvest, and samples were stored at −80°C. Western samples were shipped overnight on dry ice and immediately stored at −80°C until RNA extraction. From each plot, measurements for total seed protein and oil contents were performed using an Infratec 1241 Grain Analyzer (FOSS North America, Eden Prairie, MN, USA) at the Agriculture and Agri-Food Canada Ottawa Research and Development Centre. For additional phenotypic information, see [25].

RNA extractions using SPLIT Total mRNA Extraction Kit (Lexogen, Vienna, Austria) were performed on approximately 200 mg of crushed leaf tissue from each sample according to the manufacturer’s instructions. RNA quality was tested using a NanoDrop™ 2000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), agarose gel electrophoresis (1%), TapeStation 4200 RNA ScreenTape (Agilent, Santa Clara, CA, USA), and 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) at Génome Québec (Montréal, QC, Canada) and the Ottawa Research and Development Centre (Ottawa, ON, Canada). RNA integrity number (RIN) values of at least 6.5 and a Q30 score of at least 36 were selected for library preparation. Spike-in RNA variants (SIRVs) (Lexogen, Vienna, Austria) were integrated into the RNA samples as controls to monitor and compare key parameters (such as sensitivity and quantification); the E0 SIRV mix was used, containing 69 different isoform variants with known sequences at equal molar concentrations.

Paired-end sequencing was carried out using the Illumina HiSeq 4000 platform (Illumina, San Diego, CA, USA) at Génome Québec to create cDNA libraries for each sample. RNA-seq data were assessed using dupRadar (Sayols et al., 2016) (v3.16, Biberach an der Riß, Germany; Bioconductor, R) for duplication rate quality control. Normalization of reads was carried out at the individual sample level using edgeR (Robinson et al., 2010) (v3.16, Parkville, Victoria, Australia). Exploratory data analysis of normalized reads was performed using R.

QualiMap (García-Alcalde et al., 2012) (v2.2.1, Berlin, Germany) was used as a quality control step for sequence data feature alignment (genes and transcripts). Preseq (Daley et al., 2020) (v3.1.1, Los Angeles, CA, USA; Bioconductor, R) was used to estimate the number of distinct reads from each RNA-seq library. RSeQC (Wang et al., 2012) (v4.0.0, Nanjing, China; Bioconductor, R) was the primary tool used for comprehensive evaluation of the RNA-seq read data through the calculation of semantic read distribution of a sample, distance between reads, duplication presence, and junction saturation.

The Canadian Centre for Computational Genomics uses an in-house framework program, GenPipes (Bourgey et al., 2019), to perform the following major processing steps. Trimmomatic (Bolger et al., 2014) (v0.36, Jülich, Germany) was used to remove adaptor sequences and low quality score bases (phred score <30). Trimmed reads were aligned to the soybean genome (Glycine_max_v2.1, INSDC Assembly GCA_000004515.4, Jul 2018), using STAR (Dobin et al., 2013) (v2.7.7a, Menlo Park, CA, USA) under the command –runMode alignReads after generating index files from the aforementioned genome. HTSeq (Anders et al., 2015) (v0.12.3, Heidelberg, Germany) was used to obtain read counts using the following options: “-m intersection-nonempty”.

DE analysis was performed using DESeq2 (Love et al., 2014) (v3.16, Heidelberg, Germany) with negative Binomial GLM fitting and Wald statistics: nbinomWaldTest. To transform expression data to be expressed as log₂FC, “ashr” (Stephens, 2017) was used. All datasets were trimmed to an adjusted p-value <0.01. For DE analysis, Ottawa samples were used as the control, and the three western locations were each used as the experimental data; the log₂ fold change (log₂FC) difference in expression data reflects a change occurring in the western-grown relative to eastern-grown cultivars. Identical genotypes were compared for each DE analysis, and comparisons were not made across different genotypes.

All DE datasets were amended with the corresponding information from the SoyBase Genome Annotation Source v2.0 (https://soybase.org/genomeannotation/), which includes annotation data from BLASTP, TAIR10, GO, Panther, PFAM, and KOG for all genes. In this study, both top-down and bottom-up analyses were used to assess DE genes between East and West. “Top-down” and “bottom-up” are descriptive terms for the direction of data analysis. To clarify, the top-down analysis uses the holistic dataset and investigates the DE genes from a bird’s-eye perspective without any specific functional selection—i.e., which genes (and their ontologies) are DE between East and West at the given cut-off criteria (p-value <0.01, |log2FC| ≥ 1.5). The top-down approach was used to holistically search the DE data for the most consistent DE genes between East and West; DE was cross-compared across all 30 datasets for most of the consistently (30/30 datasets) DE genes. The bottom-up analysis describes a different approach to the data: an investigation in which we search only for genes involved in one specific pathway of interest (the Asp-Ala-Glu pathway). The bottom-up approach was used to search through the DE data to identify genes with “asparagine”, “aspartate”, “alanine”, “glutamate”, and “oxaloacetate” as a component of their annotation (BLASTP, TAIR10, GO, Panther, PFAM, or KOG). With the use of a short bash script, all 30 DE datasets were searched for any gene with these keyword identifiers and were short-listed for pathway-specific analysis. The purpose of this bottom-up analysis was to provide a comprehensive investigation of DE genes within this pathway and provide insight into the underlying molecular mechanisms influencing low seed protein phenotypes in western Canada. A p-value <0.01 and a log₂FC change in expression of at least 1.5 (|log2FC| ≥ 1.5) were considered significantly DE for both the top-down and bottom-up approaches. Genes that are DE in 15 (50%) or more datasets were selected for downstream analysis. Because the SoyBase annotation database uses BLAST descriptions from 2014, an up-to-date BLAST search was carried out on all resultant top-down and bottom-up genes and is included alongside the SoyBase annotation output in Supplementary Table 1 and Supplementary Table 2 .

For the top-down analyses, GO term enrichment was assessed using the SoyBase GO Term Enrichment Tool (https://soybase.org/goslimgraphic_v2/dashboard.php) for the genes commonly upregulated and commonly downregulated across all 30 DE datasets. Enrichment was calculated from the ratio of expressed DE genes for a particular GO term to the expected number of DE genes for said term based on the total number of known associated genes in the full GO database.

Heatmapping was performed using Heatmapper (Babicki et al., 2016) using log-normalized read counts across all samples in this study calculated using R. Clustering was calculated using average linkage, and distance matrices were calculated using Pearson’s coefficient. The row Z-score normalizes expression data to improve visualization of heatmap data trends; this score is calculated by (gene expression value in sample of interest) − (mean expression across all samples)/(standard deviation) (Anders and Huber, 2010).

Gene IDs were converted to their corresponding National Center for Biotechnology Information (NCBI) ID numbers, which were then mapped using the Kyoto Encyclopedia of Genes and Genomes (KEGG) release v106.0 (https://www.kegg.jp/) for pathway enrichment using the Mapper Search functions, with Glycine max (gmx) as the organism identifier.

Chromosome positioning data for all 20 G. max chromosomes were extracted from SoyBase GWAS-based QTL database for all seed protein and oil QTLs (https://soybase.org/GWAS/list.php). The positional information for each of the 59 bottom-up genes of interest was also extracted from SoyBase and used to determine which genes fall within large QTL regions or in very close proximity to single-point QTLs. With the use of MapChart v2.32 (Voorrips, 2002), QTL and gene positional information were mapped to each chromosome and color-coded to be most informative.

There was a total of 4,047,045,039 reads over all 87 RNA-seq datasets (10 lines, four locations, and three replicates per line). Missing replicates are Saskatoon line 1 replicate 1, Saskatoon line 4 replicate 2, and Saskatoon line 10 replicate 3, which did not pass quality control (QC) upon repeated attempts. Table 1 summarizes the average seed protein and oil contents as a percentage of the total seed content. Included in Table 1 is the cumulative read depth across the three replicates per sample. Figure 1 shows the principal component analysis (PCA) of the transcriptome data for each replicate in East and West locations, organized by line (1–10). Across all samples, PC1 described a 66% variance, and PC2 described a 7% variance ( Figure 1 ). From Figure 1 , there is a clear distinction between East (cyan) and West (green, purple, and red) RNA-seq variability, indicating the RNA-seq data from the East are highly different than the data from the West, and the data from the three West locations cluster closely together, suggesting similar variability.

The top-down analysis showed significant enrichment of biosynthesis of secondary metabolites (gmx01110), which led to the downstream investigation of select sub-pathways. Among these sub-pathways, one pathway of interest was the alanine, aspartate, and glutamate (Ala-Asp-Glu) metabolism pathway (gmx00250), which was of particular interest because of the known relationship between nitrate assimilation during development and seed protein content at maturity (Hernández-Sebastià et al., 2005; Pandurangan et al., 2012). Asparagine synthesis and hydrolysis are components of the alanine, aspartate, and glutamate (Ala-Asp-Glu) metabolism pathway.

For the bottom-up analyses, significance criteria were loosened to include genes DE in a minimum of 15 of the 30 DE datasets, rather than all 30 datasets as used in the top-down analysis. This was to expand the DE data to include genes outside of the top-down lists. Significance criteria were maintained at a log₂FC of at least 1.5 and p-value <0.01. DE data were searched for any DE genes with annotations including the terms “alanine”, “aspartate”, “glutamate”, “asparagine”, and “oxaloacetate”, which fit these significance criteria. Supplementary Table 2 summarizes the log₂FC in expression for all significantly DE genes with these annotations. Figure 5 shows the relative expression data across all samples in this study as a heatmap; Pearson’s coefficient relationship between genes in each list was used to organize the heatmap. These heatmaps provide a visual summary of the expression data and the relationships between the genes, while the information in Supplementary Table 2 provides the specific log₂FC DE data (|log₂FC| > 1.5).

The PCA plots in Figure 1 show a clear separation between variability of expression in the East from variability of expression across all three western locations, implicating that all lines across the western-grown soybeans are behaving similarly to the others from the West and all eastern-grown soybeans are behaving similarly to others in the East. These plots were constructed using RNA-seq variability as the principal components, and with these results, we observe a clear difference in expression variability between the two geographic areas (East and West). Considering the fact that RNA-seq variability was used to assess identical genotypes grown in four different environments, it can be confidently concluded that across all 10 lines, soybeans in the West show differential transcriptomics than eastern-grown counterparts. Each plot in Figure 1 represents a different genotype; colored data points represent replicates from each location, and gray data points correspond with the data on all other plots (all other genotypes). It is evident that all samples from the West are behaving more or less the same across all 10 genotypes, and the same can be said for the East. Indeed, many factors, both genetic and environmental, cumulatively influence resultant seed protein content (Wang et al., 2019); thus, we used big-data ontology, pathway, and QTL analyses to refine genes that were consistently DE between East and West across 30 individual datasets. The clustering of the transcriptome data from the three West locations is clearly separate from the data from the East across all 10 genotypes, indicating that location biases are minimized.

Among enriched pathways, biosynthesis of secondary metabolites (gmx01110) was prominently DE between samples grown in the two regions. GO and KEGG analyses identified different amino acid biosynthetic pathways were differently regulated between East and West, including aromatic amino acids (GO:0009072; gmx00400), sulfur-containing amino acids (GO:0000096; gmx00270), glycine (GO:0006546), and glutamate aspartate and asparagine (GO:0006537; GO:0033345; GO:0008734; gmx00250). A focus on the Ala-Asp-Glu biosynthesis pathway (gmx00250) was chosen because of the relationship between nitrogen assimilates, asparagine, and seed protein at maturity. Previous investigations into the DE between vegetable soybean and grain soybean found that the Ala-Asp-Glu metabolic pathway was highly enriched, as well as fatty acid biosynthesis and metabolism, carbon (starch and sucrose) metabolism and transport, arginine and proline metabolism, and glycolysis/gluconeogenesis, all of which influence the attributes (including protein) of the resulting seed (Chen et al., 2022).

During embryo development in plants, sucrose provides a source of carbon, and glutamine and asparagine are the main nitrogenous assimilate sources (Rainbird et al., 1984). Asparagine plays a key role in the source (root nodules)–sink (seeds, mainly) translocation relationship (Lam et al., 1996). Asparagine has a relatively high nitrogen:carbon ratio and is biochemically stable, which make it ideal for nitrogen transport and storage. A careful balance exists between asparagine biosynthesis and degradation to maintain asparagine concentration. Asparagine represents up to 50% of the total free amino acids in the developing cotyledon (Hernández-Sebastià et al., 2005). AS is the major enzyme responsible for synthesizing asparagine. Typically, two or more AS genes are found in plants. High AS activity in the cotyledons of the germinating seed, as well as in mature root nodules, supports the idea that asparagine acts as a nitrogen transport system in legume plants (Lam et al., 1996). AS generates asparagine from aspartate by using glutamine or ammonia as a substrate for the transfer of the amide group to aspartic acid in an ATP-dependent reaction catalyzed by magnesium (Lea et al., 2007). AS proteins are categorized as either AS-A (a.k.a. AsnA) or AS-B (a.k.a. AsnB); AS-B family proteins, found in both prokaryotes and eukaryotes, can use both ammonia and glutamine as a nitrogen donor but prefer glutamine (Manhas et al., 2014). Glutamine-dependent AS is the main asparagine biosynthesis pathway in plants (Lam et al., 1996). ASPGs are ubiquitous across all domains of life. ASPG breaks down the isoaspartyl peptide bond in asparagine to aspartate and ammonia, which are then reassimilated through the glutamine synthase/glutamate synthase cycle (Gomes and Sodek, 1984; Haga and Sodek, 1987). In soybeans, ASPG activity is directly associated with a reduction in free asparagine by up to 18% while simultaneously increasing the amount of aspartate by up to 60% (Pandurangan et al., 2012). ASPG activity was also associated with a reduction in total nitrogen by 9%–13% and an increased concentration of seed oil by 5%–8% in soybeans (Pandurangan et al., 2012).

From this study, we observe differences in the expression of genes related to asparagine metabolism between soybeans grown in eastern and western Canada. The western-grown soybeans show downregulation of AS compared to eastern-grown soybeans ( Figure 5A ; Supplementary Table 2 ). This indicates that soybeans in the West are not synthesizing asparagine to the same degree as soybeans in the East, which may be directly attributed to the seed protein content at maturity. The western soybeans showed consistent upregulated expression of ASPG compared to the eastern counterparts. One of the most persistently upregulated ASPG genes, Glyma.17G238200, falls within a large protein QTL, seed protein 9-g5 ( Figure 7 ; Supplementary Table 4 ). Further, the most persistently downregulated asparagine-related gene, Glyma.15G073100 (stem-specific protein TSJT1), is within the major oil QTL, seed oil 11-g5 on chromosome 15 ( Figure 7 ; Supplementary Table 4 ), another chromosome highly enriched for seed protein and oil QTLs. This might suggest that these genes are linked to seed protein and oil contents. AS and TSJT1 both have annotations that include “asparagine synthetase”; however, the sequence-based analysis uncovered a major difference in their protein sequences: TSJT1 protein sequences do not contain the GATase activity position 2 Cys, and the Wali7 domain-containing protein sequence has a Glu at amino acid position 2 ( Supplementary Figure 2 ). AS genes Glyma.11G171400, Glyma.11G170300, and Glyma.18G061100 are extremely downregulated in Brandon and Morden, but not found to be downregulated in Saskatoon ( Supplementary Table 2 ). One TSJT1-encoding gene, Glyma.15G073100, was downregulated across all 30 datasets, while the other TSJT1-encoding genes were found to be downregulated mostly in Brandon and Morden (Glyma.13G279200) or Saskatoon (Glyma.07G257700).

The concentration of protein in mature soybeans is strongly associated with free asparagine in the plant during development, making it an ideal pathway for further investigation (Pandurangan et al., 2012). With the observations made in this study that ASPG is highly upregulated in the West and AS is highly downregulated in the West, it is entirely plausible that differences in asparagine metabolism are influencing the seed protein accumulation difference between eastern- and western-grown soybeans. It would be of interest to soybean breeding programs to consider increasing AS expression as an engineering target when designing high-protein soybean lines. AS1 overexpression in A. thaliana resulted in increased free asparagine levels and increased seed protein concentration (Lam et al., 2003). An increase in AS1 expression in soybean leaves had a positive correlation to seed protein concentration (Wan et al., 2006). In soybean roots, increased AS1 expression was correlated with an increased ratio of asparagine:aspartate in xylem sap headed to shoots, implicating more asparagine being transported to aerial tissues (Antunes et al., 2008).

Increased expression of asparagine aminohydrolases in western soybeans indicates that these plants are breaking down available asparagine to recycle the components, most specifically the nitrogen. When nitrogen is limited, hydrolyzing asparagine provides a source of nitrogen to be redirected into other processes. Because of the central intermediary relationship between alanine and/or serine and asparagine transamination (both amino acids can act as a substrate), increased ASPG expression, which leads to a reduction in freely available asparagine (Pandurangan et al., 2012), could logically be associated with equally proportionate increases in serine and alanine. Interestingly, in the data presented in this study, an increase in the gene encoding a serine–glyoxylate aminotransferase 2 (Glyma.18G116900) was observed, while multiple genes encoding alanine–glyoxylate aminotransferases (Glyma.01G129400, Glyma.03G040600, and Glyma.18G021000) were highly downregulated ( Figure 5B ; Supplementary Table 2 ). This means that different activities between two different enzymes of the same enzyme component (E.C.2.6.1.44) are simultaneously ongoing in western-grown soybeans, as presented in Figure 6 by the half-blue-half-red box.

The asparagine synthesis pathway within the Ala-Asp-Glu metabolism pathway shows two other enzyme components that directly influence asparagine biosynthesis/metabolism: E.C.6.3.1.1 and E.C.3.5.1.38 ( Figure 6 ). E.C.6.3.1.1 is an AS-A family AS that is an aspartate–ammonia ligase and was not found to be DE within our data. This is expected because plant AS is of the AS-B family of AS proteins and primarily uses glutamate as a N donor source rather than aspartate (Manhas et al., 2014). The KEGG pathway depicted in Figure 6 KEGG shows enzyme components that are known to be in soybeans by highlighting respective boxes in green. E.C.3.5.1.38 appears to be prokaryotic in nature as indicated by available information on KEGG and BRENDA enzyme databases. It should be noted that as a result of the decrease in asparagine, western-grown soybeans could be compensating by increasing the expression of other nitrogen-rich amino acid (arginine and lysine) metabolizing enzymes (Pandurangan et al., 2012), which were not covered in our bottom-up analysis and might serve as an interesting area for further research.

Alanine is one of the central intermediates in amino acid metabolism and a substrate of asparagine transaminase, the enzyme responsible for transferring the α-amino group between asparagine and glycine, alanine, serine, and homoserine (Pandurangan et al., 2012; Zhang et al., 2013; Gaufichon et al., 2015). There is a central intermediary relationship between asparagine transamination and alanine and/or serine in that both amino acids can act as a substrate (Pandurangan et al., 2012). Overall, the results from the alanine investigation indicate the downregulation of alanine-related genes, with six downregulated genes and one upregulated gene common across at least 50% of the DE datasets ( Figure 5B ; Supplementary Table 2 ). Significant downregulation of Glyma.01G129400, Glyma.03G040600, and Glyma.18G021000, three alanine–glyoxylate aminotransferases, was observed in the West ( Figure 5B ; Supplementary Table 2 ). An increase in ASPG leads to a decrease in asparagine and potentially an increase in alanine, which may in part explain the downregulation of alanine-related genes as a whole ( Figure 5B ). In Figure 6 , E.C.2.6.1.44 is both up- and downregulated in the West (half-blue-half-red box). Western-grown soybeans appear to be upregulating the alanine-related gene, Glyma.18G116900, encoding serine–glyoxylate aminotransferase 2 (E.C.2.6.1.44), which is also involved in serine-pyruvate transaminase activity (GO:0004760). Simultaneously, these soybeans are downregulating the expression of two alanine–glyoxylate aminotransferases (Glyma.01G129400 and Glyma.03G040600) ( Figure 5B ; Supplementary Table 2 ). In addition to the Ala-Asp-Glu pathway, enrichment for the cyanoamino acid metabolism pathway (gmx00460) is a result of the two ASPG genes (Glyma.17G238200 and Glyma.04G042100) ( Supplementary Table 3 ). High-protein soybean genotypes were found to have higher amounts of free asparagine and alanine in developing embryos than in low-protein genotypes (Hernández-Sebastià et al., 2005). Further, freely available 3-cyanoalanine was found to be highly correlated with seed protein and/or oil in soybeans (Wang et al., 2019). In an investigation into the genetic shift in soybeans over 24 years, it was found that newer soybean cultivars had a decrease in seed protein, alanine, and serine (de Borja Reis et al., 2020). The relationship between alanine, serine, and asparagine in soybean seed protein accumulation remains elusive, and further investigations into the relationship between these amino acids and protein should be explored.

The aspartate-family amino acid sub-pathway functions as a regulatory metabolic link with the tricarboxylic acid (TCA) cycle, biologically significant under extreme stress conditions, which deplete cellular energy (Galili, 2011), making it an essential metabolite for plant growth and stress acclimation (Han et al., 2021). Aspartate-family amino acids (lysine, threonine, methionine, and isoleucine) are synthesized in plants using aspartate as a central amino acid (Galili, 2011). In this study, Glyma.05G007600 encoding l-aspartate oxidase (E.C.1.4.3.16) was significantly upregulated in the West. This gene is physically close on chromosome 5 to many known oil QTLs ( Figure 7 ; Supplementary Table 4 ). Chromosome 5 is highly enriched for QTL influencing protein and oil; however, the molecular mechanisms driven by these loci remain largely unknown (Wang et al., 2019). The closeness in proximity suggests that they are tightly linked, and following recombination, it would be advantageous for these genes to remain together in future progeny.

Glyma.19G159300 (lifeguard 4 protein in G. max) was upregulated in 25 of 30 datasets ( Figure 5D ); this gene is most closely related to the gene encoding the glutamate-binding protein in A. thaliana and is also closely related to the inhibitor of the apoptosis-promoting BAX1 protein in A. thaliana ( Supplementary Table 2 ). Further investigations into the glutamate binding potential of this gene and the putative role it plays in signaling would be of merit to understanding the reason(s) for significant upregulation of Glyma.19G159300. It is likely that this protein plays a role in signaling and proliferation, potentially inhibiting apoptosis. Increasing expression of genes related to glutamate signaling and apoptosis suggests that western-grown soybeans are differently controlling cell proliferation compared to those in the East.

Glutamate has a remarkably wide range of biological roles because of the central position it plays in metabolism. It is suggested that glutamate compensates for the reduction in freely available asparagine by serving as a metabolite (nitrogen) storage molecule and behaves as an organic nitrogen signal in seedlings (Gutiérrez et al., 2008; Pandurangan et al., 2012). This study uncovered 30 glutamate-related genes that are DE in at least 50% of the datasets. Within the glutamate data, genes for proline dehydrogenases (Glyma.13G049700, Glyma.19G042900, and Glyma.19G043000) were among the most upregulated in the West. Proline dehydrogenase catalyzes the oxidation of l-proline to δ1-pyrroline-5-carboxylate, which provides a source of free electrons for transport (Servet et al., 2012). The two proline dehydrogenase genes on chromosome 19 both fall within a major seed protein QTL, seed protein 9-g6 ( Figure 7 ). Proline dehydrogenase catabolizes proline while simultaneously playing roles in energy, shuttling redox potential, and production of reactive oxygen species (ROS) to reach cellular homeostasis, adapt to the environment, and carry out physiological and pathological processes (Servet et al., 2012). A methionine synthase (Glyma.16G038300) was also highly upregulated in the West, which ultimately influences the sulfur-containing amino acid content of the developing seed. Methionine synthase is responsible for catalysis of 5-methyltetrahydropteroyltri-l-glutamate and l-homocysteine into l-methionine + tetrahydropteroyltri-l-glutamate in the synthesis of methionine and glutamate (Whitfield et al., 1970), a unique feature of some organisms that are able to convert Cys to Met under specific circumstances (Brosnan and Brosnan, 2006). Upregulation of this methionine synthase gene may suggest that western-grown soybeans are increasing glutamate in an attempt to compensate for a lack of free asparagine. If, coincidently, this also results in more Met production from l-homocysteine, the overall protein quality (in terms of 11S and 7S globulins) could be improved. Improved seed protein quality is an important consideration for soybean agriculture, particularly in regions where environmentally influenced decreases in seed protein levels are prominent (i.e., western Canada), and western-grown soybeans were found to have higher 11S:7S values than eastern-grown soybeans (Cober et al., 2023).

Also within the glutamate data are a number of DE glutamate receptor genes both upregulated (Glyma.07G226400, Glyma.10G213200, Glyma.13G049700, Glyma.13G093300, Glyma.13G172100, and Glyma.13G233400) and downregulated (Glyma.13G233000, Glyma.13G233300, and Glyma.17G067200) in the West ( Figure 5D ; Supplementary Table 2 ). This points to specific genes involved in DE glutamate signaling between eastern- and western-grown soybeans. The specific ligands for these receptors would be an interesting area for further research on these genes.

As previously stated, the Ala-Asp-Glu metabolism pathway coordinates a metabolic link to the TCA cycle, directly influencing energy production or depletion (Galili, 2011). Oxaloacetate is both a product and a beginning component of the TCA cycle; levels of oxaloacetate give an indication of the ongoing level of energy metabolism. The data uncovered in this study indicate major upregulation of oxaloacetate-related genes in western-grown soybeans; of eight significantly DE genes, seven were upregulated, and one was downregulated ( Figure 5E ; Supplementary Table 2 ). All of the DE genes related to oxaloacetate are malic enzymes (E.C.1.1.1.40, pyruvate metabolism pathway gmx00620; Supplementary Table 3 ), with the exception of Glyma.15G055600. Glyma.15G055600, a 2-oxoglutarate/malate translocator-like protein, was found within the oxaloacetate-related gene list ( Supplementary Table 3 ) and the only oxaloacetate-related gene that does not map to any pathway using KEGG ( Supplementary Table 4 ). However, Glyma.15G055600 falls within the same major oil QTL as one of the most persistently downregulated asparagine-related genes (Glyma.15G073100), seed oil 11-g5 on chromosome 15 ( Figure 7 ; Supplementary Table 4 ).

Malic enzyme is one of the key enzymes linked to fatty acyl chain biosynthesis. l-Malic acid was found to be highly correlated with protein and oil contents in soybeans (Wang et al., 2019). A significant amount of pyruvate, the precursor of acetyl-CoA synthesis for lipid biosynthesis, is produced as a result of malic enzyme activity in soybeans (Allen and Young, 2013). In western-grown soybeans in this study, malic enzyme activity was highly upregulated, which suggests that these plants are likely producing higher levels of acetyl-CoA for fatty acid production. Because of this, malic enzyme activity is almost certainly one of the regulatory mechanisms underlying the inverse relationship between seed protein and oil contents in soybeans. This mechanism serves are an optimal target for genetic engineering/control of the decision between lipid and protein biosynthesis. Increasing expression of malic enzyme would make a molecular conduit to directing nitrogen and carbon toward lipid biosynthesis and away from protein biosynthesis, serving as a molecular switch for the accumulation of major seed storage biomolecules (Morley et al., 2023).

The ontologies sulfur amino acid metabolic process (GO:0000096), sulfur compound biosynthetic process (GO:0044272), sulfur compound metabolic process (GO:0006790), iron-sulfur cluster binding (GO:0051536), and more were all enriched and overrepresented within the top-down GO analysis ( Figure 3 ; Supplementary Table 1 ). Further, the sulfur-containing amino acid Cys and Met metabolism pathway (gmx00270) was enriched within the DE data. These observations that transcription of genes related to the biosynthesis of sulfur-containing amino acids (Cys and Met) are DE between East and West likely play a role in the differences in seed protein quality seen between western- and eastern-grown soybeans (Cober et al., 2023). Sulfur-containing amino acids are essential to the formation of 11S storage proteins (glycinins) in soybean. Sedimentation coefficients (0.5-M ionic strength) are used to categorize seed storage proteins into 2S, 7S, 11S, and 15S fractions, of which the 11S and 7S fractions account for the majority of seed storage proteins (40% and 30% of seed storage protein, respectively) (Peng et al., 1984; Tsukada et al., 1986). Cys and Met are limited resources in soybeans, and tight control of biosynthesis of these amino acids is advantageous to glycinin production in soybeans. Because Cys and Met are vital to glycinin biosynthesis, the genes influencing expression and accumulation of sulfur-containing amino acids very likely influence glycinin accumulation (nutrient reservoir activity) and therefore seed protein content.

Enriched among the data were many other pathways of interest That almost certainly influence seed protein content, including fatty acid metabolic processes (fatty acid biosynthetic process GO:0006633; lipid transport GO:0006869; lipid binding GO:0008289; gmx01212), circadian rhythm (circadian rhythm GO:0007623; regulation of circadian rhythm GO:0042752; gmx04712), nutrient storage (nutrient reservoir activity GO:0045735), and carbohydrate metabolism (carbohydrate metabolic process GO:0005975; carbohydrate-binding GO:0030246; regulation of carbohydrate metabolic process GO:0006109) ( Figures 3 , 4 ; Supplementary Table 1 ). Western locations in this study are also further North than the Eastern locations and as a result experience longer photoperiods; the circadian rhythm of soybeans in these locations is almost certainly going to exhibit differences. An insertion/deletion in Glyma.20G85100, a circadian clock gene, was found to nearly perfectly correspond with high/low protein alleles for a QTL on chromosome 20, cqSeed protein-003 (Fliege et al., 2022). While Glyma.20G85100 was not DE in our data, genes with circadian rhythm among their ontologies were overrepresented in both upregulated (20) and downregulated genes (23) ( Figure 3 ; Supplementary Table 1 ). Response to UV (GO:0009411) is upregulated in western-grown soybeans, likely in response to reduced cloud cover in the prairies. Ontologies for microtubule cytoskeleton organization (GO:0000226), microtubule motor activity (GO:0003777), cytokinesis by cell plate formation (GO:0000911), cell wall biogenesis (GO:0042546), cell cycle (GO:0007049), and spindle assembly (GO:0051225) are highly enriched in the upregulated genes in our study ( Figure 3 , Supplementary Table 1 ). Soybeans in western Canada were found to be significantly taller than eastern-grown soybeans (Cober et al., 2023), which is likely influenced by these genes, suggesting the plants are spending more energy increasing in height than producing/filling seeds. Genes related to water stress were also found to be DE between East and West ( Figure 3 ); response to desiccation (GO:0009269), water transport (GO:0006833), response to water depravation (GO:0009415), and water channel activity (GO:0015250) are all likely the result of the lower precipitation and lower relative humidity in the West.

Protein and oil compete for space in the seed, which results in a pushing/pulling relationship between these seed storage macronutrients (Breene et al., 1988; Clemente and Cahoon, 2009). Previous investigations into seed protein, oil, and yield identified an influential QTL on chromosome 20 between Satt496 and Satt239 (Chung et al., 2003); however, there were no DE genes found between these markers in our data. The relationship between seed storage and metabolism is entangled by the upstream carbohydrate metabolism decision-making steps that lead to protein and/or oil biosynthesis (40% and 20%, respectively) while also maintaining a proportion (~35%) of the seed space for stored carbohydrates (Liu, 1997). Seed storage proteins (glycinins, vicilins, and cupins) have the ontology nutrient reservoir activity (GO:0045735). Nutrient reservoir activity (GO:0045735) was overrepresented in the top-down downregulated genes but underrepresented (though present) in the upregulated GO enrichment ( Supplementary Table 1 ). Known genes for glycinins and β-conglycinins (Glyma.03G163500 GY1, Glyma.03G163500 GY2, Glyma.19G164900 GY3, Glyma.10G037100 GY4, Glyma.13G123500 GY5, Glyma.10G246300 CG-1 α′1, Glyma.20G148400 CG-2 α2, Glyma.20G148300 CG-3 α1, Glyma.20G146200 CG-4 β1, and Glyma.20G148200 CG-4 β2) were not found to be significantly DE across a majority of lines but were found to be DE in some lines ( Supplementary Table 1 ). β-Conglycinin genes are reported to only be expressed in seeds in early embryogenesis; transcription spikes at mid-maturation and decreases before dormancy (Harada et al., 1989). These genes are not expressed in cotyledons or at maturity; thus, it is reasonable that we do not see notable DE between East and West in the present data; leaf tissue at the R5 (seed filling) stage was used for RNA-seq. A similar investigation into RNA-seq analysis of soybean pod data would provide further information into the DE of glycinins and β-conglycinins. The findings in this study share similarities with a similar study conducted in China on high and low seed protein varieties, which found DE of genes involved in the biosynthesis of amino acids and secondary metabolites, carbon metabolism, lipid metabolism, phenylpropanoid biosynthesis, and plant hormone signal transduction (Xu et al., 2022), although we examined DE of individual varieties across geographies and not DE between varieties.