
==== Front
Plant Cell
Plant Cell
plcell
The Plant Cell
1040-4651
1532-298X
Oxford University Press US

39056474
10.1093/plcell/koae202
koae202
Research Article
AcademicSubjects/SCI01270
AcademicSubjects/SCI01280
AcademicSubjects/SCI02286
AcademicSubjects/SCI02287
AcademicSubjects/SCI02288
Transcriptomic response to nitrogen availability reveals signatures of adaptive plasticity during tetraploid wheat domestication
https://orcid.org/0000-0002-1950-954X
Pieri Alice Department of Agricultural, Food and Environmental Sciences, Marche Polytechnic University, via Brecce Bianche, Ancona 60131, Italy

https://orcid.org/0000-0002-2880-5033
Beleggia Romina Council for Agricultural Research and Economics (CREA), Research Centre for Cereal and Industrial Crops (CREA-CI), Foggia 71122, Italy

https://orcid.org/0000-0001-8980-3034
Gioia Tania School of Agricultural, Forestry, Food and Environmental Sciences, University of Basilicata, Potenza 85100, Italy

https://orcid.org/0000-0003-4856-1336
Tong Hao Bioinformatics Department, Institute of Biochemistry and Biology, University of Potsdam, Potsdam 14476, Germany
Systems Biology and Mathematical Modeling Group, Max Planck Institute of Molecular Plant Physiology, Potsdam 14476, Germany

https://orcid.org/0000-0003-3796-8587
Di Vittori Valerio Department of Agricultural, Food and Environmental Sciences, Marche Polytechnic University, via Brecce Bianche, Ancona 60131, Italy

https://orcid.org/0000-0003-0114-7314
Frascarelli Giulia Department of Agricultural, Food and Environmental Sciences, Marche Polytechnic University, via Brecce Bianche, Ancona 60131, Italy

https://orcid.org/0000-0002-2152-4634
Bitocchi Elena Department of Agricultural, Food and Environmental Sciences, Marche Polytechnic University, via Brecce Bianche, Ancona 60131, Italy

https://orcid.org/0000-0002-1249-5165
Nanni Laura Department of Agricultural, Food and Environmental Sciences, Marche Polytechnic University, via Brecce Bianche, Ancona 60131, Italy

https://orcid.org/0000-0002-1636-0719
Bellucci Elisa Department of Agricultural, Food and Environmental Sciences, Marche Polytechnic University, via Brecce Bianche, Ancona 60131, Italy

https://orcid.org/0000-0001-8775-1541
Fiorani Fabio Institute of Biosciences and Geosciences (IBG-2): Plant Sciences, Forschungszentrum Julich GmbH, Julich 52428, Germany

https://orcid.org/0000-0003-1704-2541
Pecchioni Nicola Council for Agricultural Research and Economics (CREA), Research Centre for Cereal and Industrial Crops (CREA-CI), Foggia 71122, Italy

https://orcid.org/0009-0005-3014-0869
Marzario Stefania School of Agricultural, Forestry, Food and Environmental Sciences, University of Basilicata, Potenza 85100, Italy

https://orcid.org/0000-0003-2885-5468
De Quattro Concetta Department of Biotechnology, University of Verona, Strada Le Grazie 15, Verona 37134, Italy

https://orcid.org/0000-0003-1570-8370
Limongi Antonina Rita Department of Biotechnology, University of Verona, Strada Le Grazie 15, Verona 37134, Italy

https://orcid.org/0000-0002-9573-0510
De Vita Pasquale Council for Agricultural Research and Economics (CREA), Research Centre for Cereal and Industrial Crops (CREA-CI), Foggia 71122, Italy

https://orcid.org/0000-0002-6101-1550
Rossato Marzia Department of Biotechnology, University of Verona, Strada Le Grazie 15, Verona 37134, Italy

https://orcid.org/0000-0003-0369-8777
Schurr Ulrich Institute of Biosciences and Geosciences (IBG-2): Plant Sciences, Forschungszentrum Julich GmbH, Julich 52428, Germany

https://orcid.org/0000-0002-2925-2258
David Jacques L AGAP, Univ. Montpellier, CIRAD, INRAE, Institut Agro, Montpellier 34060, France

https://orcid.org/0000-0003-2671-6763
Nikoloski Zoran Bioinformatics Department, Institute of Biochemistry and Biology, University of Potsdam, Potsdam 14476, Germany
Systems Biology and Mathematical Modeling Group, Max Planck Institute of Molecular Plant Physiology, Potsdam 14476, Germany

https://orcid.org/0000-0001-9598-3131
Papa Roberto Department of Agricultural, Food and Environmental Sciences, Marche Polytechnic University, via Brecce Bianche, Ancona 60131, Italy

Author for correspondence: r.papa@univpm.it (R.P.), romina.beleggia@crea.gov.it (R.B.)
The author responsible for distribution of materials integral to the findings presented in this article in accordance with the policy described in the Instructions for Authors (https://academic.oup.com/plcell/pages/General-Instructions) is: Roberto Papa (r.papa@univpm.it).

Conflict of interest statement. None declared.

9 2024
26 7 2024
26 7 2024
36 9 38093823
10 10 2023
17 7 2024
31 7 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of American Society of Plant Biologists.
2024
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The domestication of crops, coupled with agroecosystem development, is associated with major environmental changes and provides an ideal model of phenotypic plasticity. Here, we examined 32 genotypes of three tetraploid wheat (Triticum turgidum L.) subspecies, wild emmer, emmer, and durum wheat, which are representative of the key stages in the domestication of tetraploid wheat. We developed a pipeline that integrates RNA-Seq data and population genomics to assess gene expression plasticity and identify selection signatures under diverse nitrogen availability conditions. Our analysis revealed differing gene expression responses to nitrogen availability across primary (wild emmer to emmer) and secondary (emmer to durum wheat) domestication. Notably, nitrogen triggered the expression of twice as many genes in durum wheat compared to that in emmer and wild emmer. Unique selection signatures were identified at each stage: primary domestication mainly influenced genes related to biotic interactions, whereas secondary domestication affected genes related to amino acid metabolism, in particular lysine. Selection signatures were found in differentially expressed genes (DEGs), notably those associated with nitrogen metabolism, such as the gene encoding glutamate dehydrogenase (GDH). Overall, our study highlights the pivotal role of nitrogen availability in the domestication and adaptive responses of a major food crop, with varying effects across different traits and growth conditions.

Nitrogen had a pivotal role in tetraploid wheat domestication and adaptive responses, leading to diverse impacts on gene expression and selection signatures during primary and secondary domestication.

Università Politecnica delle Marche 10.13039/501100005758 Agrofood Science and Safety European Plant Phenotyping Network 284443 FP7 Research Infrastructures 10.13039/100011260 Horizon 2020 10.13039/501100007601 862201
==== Body
pmcIntroduction

Domestication influences the genetic diversity of animals and plants as they adapt to agroecosystems and undergo selection to meet human preferences and needs. This process is typically associated with the genome-wide loss of nucleotide diversity due to the combined consequences of selection and genetic drift, which is known as the domestication bottleneck. The loss of genetic diversity has been documented in many domesticated species by comparing them with wild relatives (Bitocchi et al. 2017). A parallel effect is the reprogramming of gene expression and the loss of expression diversity, which was first reported in the common bean (Phaseolus vulgaris; Bellucci et al. 2014) and subsequently in other domesticated plants and animals (Sauvage et al. 2017; Liu et al. 2019; Burgarella et al. 2021). Similar observations have been reported at the level of metabolic diversity (Beleggia et al. 2016).

Changes in nucleotide and gene expression diversity during the domestication of tetraploid wheat (Triticum turgidum L., 2n = 4x = 28; AABB genome) are not fully understood. Tetraploid wheat was domesticated in two well-defined phases. Primary domestication from wild emmer (Triticum turgidum ssp. dicoccoides) to emmer (Triticum turgidum ssp. dicoccum) started ∼12,000 yrs ago in the Fertile Crescent. This was followed by secondary domestication from emmer to durum wheat (Triticum turgidum ssp. durum), which started 8,000–10,000 yrs ago in the Near East and gave rise to durum wheat, the most important form of tetraploid wheat and currently the most widespread Mediterranean crop (Gioia et al. 2015; Taranto et al. 2020; Levy and Feldman 2022).

The transition from wild environments to early farming and eventually to modern high-input agroecosystems had profound ecological consequences. Throughout history, humans have employed various methods to enhance soil fertility, such as soil preparation to facilitate organic matter mineralization and the use of livestock manure, with evidence dating back to Neolithic early farming sites ∼7,900 yrs ago (Bogaard et al. 2013). However, the scale and intensity of fertilizer use have escalated over time, especially with the advent of the Haber-Bosch industrial process, which heavily relies on nonrenewable fossil fuels. Notably, the widespread overreliance on nitrogen (N) fertilizers in modern industrial agriculture can be traced back to the Donald model (Donald 1968). This model aims to optimize crop yields by minimizing intraspecific competition and providing substantial agronomic inputs, including fertilizers (Fréville et al. 2022). Today, the global application of N fertilizers, particularly to cereal crops, exceeds 80 million tons annually (Ludemann et al. 2022). N is an essential macronutrient whose availability is directly linked to crop yield and grain quality (protein content) (Barneix 2007; Howarth et al. 2008; Laidò et al. 2013), but it is also directly harmful to humans and the environment. Indeed, excess N from agricultural sources is one of the major fresh water pollutants, causing the eutrophication of aquatic ecosystems (Rockström et al. 2009). The production of industrial fertilizers contributes ∼3% of global CO2 and is a primary source of N2O (Wood and Cowie 2004). Understanding genetic variations in N acquisition, assimilation, and metabolism can therefore provide strategies for crop improvement to meet the United Nations Sustainable Development Goals (SDGs) (Plett et al. 2018; Hawkesford and Griffiths 2019).

Environmental changes that accompanied the domestication of crops over thousands of years can be tolerated by organisms that exhibit phenotypic plasticity, defined as the ability of a genotype to exhibit changes in a specific trait across different environments, and through the modulation of gene expression (Bradshaw 1965; Laitinen and Nikoloski 2019). Understanding the molecular basis of phenotypic plasticity in crops and their wild relatives can help to address the challenges faced by modern agriculture. In tetraploid wheat, phenotypic differences in below-ground and above-ground growth traits related to N availability primarily arose during secondary domestication (Gioia et al. 2015), but the relationship between N metabolism and changes in gene expression plasticity during domestication is unclear.

Here, we analyzed 32 wild emmer, emmer, and durum wheat genotypes by RNA-Seq in contrasting N availability scenarios, to investigate the potential role of N during the domestication of tetraploid wheat. Our study elucidates the subspecies-specific responses to nitrogen, as well as nucleotide and gene expression diversity during both primary and secondary domestication phases. Our results provide insight into the pivotal role of N during the domestication and adaptive plasticity of one of our major food crops.

Results and discussion

Choice of the reference genome: mapping accuracy and comparative analysis across available tetraploid wheat references

The inclusion of genotypes from three distinct tetraploid wheat subspecies (i.e. wild emmer, emmer, and durum wheat) within our panel presented challenges when comparing results across groups. In this context, the choice of the reference genome was of crucial importance. Given that reference genomes were only available for two of the three subspecies examined in our study, namely wild emmer and durum wheat, we considered using the A and B subgenomes of bread wheat (Triticum aestivum) as an outgroup reference closely related to all subspecies in the panel. This option aimed to mitigate potential biases that might arise from favoring any subspecies, ensuring a balanced comparison across the genotypes. Indeed, the A and B subgenomes of bread wheat serve as a good intermediate between emmer (wild and domesticated) and durum wheat. Bread wheat shares ancestral A and B subgenomes with the wild and cultivated emmer, originating from the same founding population. On the other hand, the A and B subgenomes of durum wheat underwent differentiation coinciding with the origin of bread wheat (Haudry et al. 2007; Levy and Feldman 2022).

However, we rigorously tested these assumptions and validated our choice by also mapping reads to currently available tetraploid wheat reference genomes: wild emmer Zavitan (Zhu et al. 2019) and durum wheat Svevo (Maccaferri et al. 2019).

We prepared 128 RNA-Seq libraries from the 4-week-old leaves of 32 tetraploid wheat genotypes representing wild emmer, emmer, and durum wheat, grown in two contrasting N conditions (N starvation and optimal N availability) (Supplementary Data Set 1). The mapping frequency across the entire genome was consistent among the three references (86% to 87%) with an average of 6.8 million mapped reads per genotype (Supplementary Data Set 1A to C). However, while the proportion of reads mapping to genic regions was similar when comparing bread wheat and wild emmer wheat (average 73%), it was lower when using the durum wheat reference genome (average 52%) primarily due to the absence of untranslated regions (UTRs) in the Svevo reference annotation (Supplementary Data Set 1A to C).

Pairwise genetic distances between the three reference genomes were computed using an alignment-free method based on the MinHash technique (Ondov et al. 2016). This method compresses large genomic sequences (the 3 entire genomes, in our case) into sketch representations, allowing for rapid similarity estimations with bounded error. We found that the Mash distance (D), which is an approximation of the mutation rate (Ondov et al. 2016), between the bread wheat Chinese Spring A and B subgenomes was ∼0.014 (P-value < 10−10) when compared with wild emmer Zavitan, and ∼0.008 (P-value < 10−10) when compared with durum wheat Svevo. The distance between durum wheat Svevo and wild emmer Zavitan was the same as that between bread wheat and wild emmer Zavitan (D ∼0.014, P-value < 10−10). From the resulting distances, the average nucleotide identity (ANI) can be extracted as D ∼1—ANI (Ondov et al. 2016), obtaining a value for ANI of ∼99% between each genome, thus confirming their close relationship.

Computing raw read counts and filtering out genes with weak expression (see Materials and Methods) resulted in 32,358 genes from the bread wheat Chinese Spring reference genome. We performed the same process using the wild emmer Zavitan and durum wheat Svevo reference genomes, resulting in 33,586 and 29,784 genes, respectively. A comprehensive comparison of the sequences of the three gene sets also showed an ANI exceeding 98%.

Overall, despite minor differences in the outcomes across different references, the overall patterns remain consistent, and the use of an outgroup species for reference might help mitigate biases and ensure fair representation of all subspecies, maintaining mapping accuracy and coverage of gene regions, making it a suitable choice for our analysis.

A greater loss of nucleotide diversity occurred during the secondary domestication of tetraploid wheat

Variant calling on the RNA-Seq of the whole panel of 32 genotypes produced 800,996 high-quality single-nucleotide polymorphisms (SNPs) including “population-specific” SNPs found only in one of the three subspecies. The number of polymorphic sites was similar in wild emmer (617,128) and emmer (613,509), but was much lower in durum wheat (425,513). Site frequency spectra for each subpopulation are provided in Supplementary Fig. S1. SNPs were categorized as “private” if found exclusively in a single subspecies or “shared” if they were distributed across two or three subspecies. We identified 190,377 common SNPs shared by all three taxa. As expected, wild emmer and emmer shared the highest percentage of SNPs (33%, 206,578). In contrast, durum wheat shared only 11% (46,352) of its SNPs with wild emmer and 17% (71,147) with emmer (Supplementary Data Set 2).

SNP principal component analysis (PCA) revealed the broad genetic structure of the three wheat taxa (Fig. 1) and confirmed that secondary domestication had a greater impact than primary domestication in differentiating the durum wheat subspecies. The 12 durum wheat genotypes are genetically very similar, forming a dense cluster that is clearly distinguishable from the wild emmer and emmer genotypes. In contrast, the wild emmer and emmer genotypes are loosely clustered, indicating a greater genetic admixture. These results are consistent with previous genetic studies on the origins of domesticated tetraploid wheat and reflect the multiple stages of domestication (Haudry et al. 2007; Luo et al. 2007; Civáň et al. 2013; Oliveira et al. 2020).

Figure 1. Principal component analysis (PCA) of 32 wheat genotypes based on single-nucleotide polymorphisms (SNPs). The first two principal components (PC1 and PC2) are shown. The three colors represent different taxa. Labels show the accession name of each genotype.

Nucleotide diversity estimates (π and θ) revealed the expected substantial loss of diversity during domestication, highlighting the greater impact of secondary domestication (Table 1). Considering π, the average nucleotide diversity of durum wheat was ∼17% lower than that of domesticated emmer, which in turn was ∼11% lower than that of wild emmer. The cumulative effect of primary and secondary domestication was a ∼26% reduction in the nucleotide diversity of durum wheat compared to its wild ancestor (Table 1).

Table 1. Nucleotide diversity estimates and diversity loss for the three wheat taxa

	Wild emmer	Emmer	Durum wheat	Loss of nucleotide diversity (%)	
Lpd	Lsd	Both	
Π	0.0050	0.0045	0.0037	11.4	16.8	26.3	
θ	0.0047	0.0040	0.0029	15.3	27.2	38.3	
Diversity loss is shown during primary domestication (wild emmer to emmer, Lpd), secondary domestication (emmer to durum wheat, Lsd), and both processes (wild emmer to durum wheat), based on average π and θ estimates of nucleotide diversity.

To ensure that our results were not biased toward the chosen reference genome, we also repeated the variant calling, PCA based on SNPs, and nucleotide diversity estimates for each subspecies using the wild emmer and durum wheat reference genomes to allow comparison with the bread wheat reference. The numbers of SNPs, nucleotide diversity estimates, and diversity loss estimates are summarized in Supplementary Data Set 2. Both the wild emmer and durum wheat references yielded fewer SNPs (604,479 and 544,406, respectively) compared to using the bread wheat reference (800,996). However, the similarity in the number of polymorphic sites between wild emmer and emmer along with the lower number of durum wheat SNPs, as well as the ratio between private and shared SNPs, was reaffirmed with both of these alternative references (Supplementary Data Set 2). The equivalence in the utilization of the three references was reinforced by identical PCA results obtained from the two sets of SNPs derived from the calls using wild emmer and durum wheat references (Supplementary Fig. S2).

In addition, the estimates of π and θ for wild emmer and emmer exhibited lower subspecies-specific values compared to the bread wheat reference (Supplementary Data Set 2). Considering π, the values obtained using the wild emmer reference were 0.0027 for wild emmer, 0.0027 for emmer, and 0.0024 for durum, whereas similar but lower values were obtained when using the durum wheat reference, especially for durum wheat (π = 0.0011). Even so, the overall trend of diversity loss confirmed previous findings using the bread wheat reference, indicating a more pronounced impact of secondary domestication. The percentage losses of π nucleotide diversity observed using the bread wheat reference (11.4 for primary domestication and 16.8 for secondary domestication) were in the range obtained using wild emmer (3.2 for primary domestication and 10.1 for secondary domestication) and durum wheat (12.5 for primary domestication and 48.1 for secondary domestication).

In summary, although we observed slight variations in SNP detection and nucleotide diversity estimates across various reference genomes, the general trends of genetic variation and diversity loss remained consistent. By detailing our methodology for testing available reference genomes, we aim to provide guidance for other researchers encountering similar challenges. In the case of tetraploid wheat, we propose that the use of an outgroup species as a reference genome does not introduce bias in gene identification and subsequent analyses. On the contrary, it may mitigate biases and ensure fair representation for all subspecies in the study.

The variability of gene expression during domestication was influenced by N availability

To quantify the diversity of gene expression in each subspecies, we calculated evolvability scores under high- and low-N availability conditions. Evolvability was estimated using the additive coefficient of variation (CVA) in read counts (Supplementary Data Set 3). In contrast to heritability, CVA is a standardized measure of additive genetic variation that is not influenced by other sources of variance (Houle 1992; Hansen et al. 2011) and is therefore well suited for comparative analysis (Garcia-Gonzalez et al. 2012; Gioia et al. 2015). We found that the CVA was the highest in wild emmer, followed by emmer and then durum wheat, indicating a decline during domestication under both N availability conditions (Fig. 2A and B; Table 2). The loss of diversity in gene expression has been observed across the domestication process in other crops, such as common bean (Bellucci et al. 2014), tomato (Solanum lycopersicum) (Sauvage et al. 2017), and sorghum (Sorghum bicolor) (Burgarella et al. 2021) as well as domesticated animal species (Liu et al. 2019). However, we observed a higher mean CVA across all three subspecies under low-N compared to high-N conditions (Fig. 2A and B; Table 2). This suggests that higher N availability promotes a more uniform gene expression pattern, whereas higher variability (plasticity) is observed during N starvation.

Figure 2. Density plots of the additive coefficient of variation (CVA) in the three wheat taxa. Comparison of the estimated density functions of the CVA in gene expression, calculated using all 32,358 genes. A) Low-N conditions. B) High-N conditions. Dashed lines represent the averaged CVA value, colored according to the different taxa.

Table 2. Mean additive coefficient of variation (CVA) in gene expression and loss of expression diversity for the three wheat taxa

	Wild emmer	Emmer	Durum wheat	Loss of expression diversity (%)	
Lpd	Lsd	Both	
CVA high N	0.062	0.056	0.048	9.1	14.5*	22.3	
CVA low N	0.076	0.063	0.056	17.6*	11.1*	26.7	
Diversity loss is shown during primary domestication (wild emmer to emmer, Lpd), secondary domestication (emmer to durum wheat, Lsd), and both processes (wild emmer to durum wheat), based on averaged CVA values calculated for all 32,358 genes.

* P < 0.001, Mann–Whitney U-test for difference between Lpd and Lsd within each N condition and difference between high N and low N within Lpd and within Lsd.

We used the contrasting N conditions of our samples to examine whether the loss of expression diversity is associated with the specific aspects of the cultivation environment, causing primary and secondary domestication to have a substantially different impact. Under high-N conditions, we observed a ∼9% loss in expression diversity in emmer compared to wild emmer (effect of primary domestication) and a ∼15% loss in durum wheat compared to emmer (effect of secondary domestication). In contrast, these losses were ∼18% and 11% under N starvation conditions, revealing twice the loss of expression diversity during primary domestication, but a lower value during secondary domestication (Table 2). All four values differed significantly from each other (Mann–Whitney U-test, P < 0.001). The opposing expression diversity profiles during domestication under high-N and low-N conditions were observed not only for overall gene expression but also for the subgroup comprising all differentially expressed genes (DEGs) and the subgroup comprising all unmodulated genes (Supplementary Table S1). The loss of expression diversity among the DEGs due to primary domestication was ∼9% and ∼15% under high-N and low-N conditions, respectively, whereas the corresponding losses due to secondary domestication were ∼18% and ∼14% (Supplementary Table S1). The loss of expression diversity among the unmodulated genes was similar to the values for overall gene expression (Supplementary Table S1).

A phenotypic study of the same accessions used in the present work has already shown that secondary domestication reduced the phenotypic diversity under high-N conditions, but the reduction was smaller and not significant under N starvation (Gioia et al. 2015). In the case of durum wheat, selection has apparently enhanced the growth response to N availability, indicating a putative focus on improving N uptake and utilization efficiency. Our expression diversity results indicate that selection may have favored specific traits and thus led to a more uniform set of cultivars, as also suggested in an earlier study based on morphological traits (Gioia et al. 2015).

Domestication and nitrogen availability played a role in the divergence of tetraploid wheat

Genetic differentiation among the three subspecies was estimated by calculating the pairwise fixation index (FST) for every gene locus in our dataset. As shown in Fig. 3A, the lowest genetic differentiation was observed between wild emmer and emmer (median FST, ∼0.08), whereas much higher genetic differentiation was found between emmer and durum wheat (median FST, ∼0.24) and, similarly, between wild emmer and durum wheat (median FST, ∼0.26). These values align with earlier findings that examined broad collections of tetraploid wheat accessions (Luo et al. 2007), suggesting an indication of the representativeness of the utilized genotypes.

Figure 3. F ST and QST distributions. A) Boxplots showing the gene locus FST distribution for every subspecies pairwise comparison. B) Boxplots showing the gene expression QST distribution for every subspecies pairwise comparison under low-N and high-N conditions, represented by empty and hatched grayscale bars, respectively. The borders of the box represent the 25th and 75th percentiles. The horizontal line in the middle of the box represents the median. Whiskers extend to the minimum and maximum values, unless a point exceeds 1.5 times the interquartile range, in which case the whisker extends to this value and values beyond are plotted as individual points (outliers).

Divergence at the transcriptomic level was estimated by calculating QST, the quantitative analog of FST, taking N availability into account as an environmental variable. Under both N conditions, we observed the same trend shown for FST (Fig. 3B). Specifically, secondary domestication had a stronger impact on differentiation (emmer vs durum wheat; median QST LN, ∼0.04; median QST HN, ∼0.15) than primary domestication (wild emmer vs emmer, median QST LN, ∼4.7 × 10−9, median QST HN ∼0.018). Interestingly, the QST distributions of every pairwise comparison showed higher values under high-N conditions compared to N starvation (Fig. 3B), suggesting that the response to N availability has been under selection during domestication and breeding substantially contributed to the differentiation of gene expression in tetraploid wheat in response to different agroecosystems.

We therefore dissected the transcriptome as a multidimensional plastic phenotype, in which the abundance and expression patterns of thousands of genes in different environmental conditions can be processed as phenotypic traits, reflecting variability in levels of gene expression. Studying the transcriptome as a plastic phenotype is a powerful approach to understand the molecular basis of plasticity and its evolutionary potential (Leinonen et al. 2013; Oostra et al. 2018). By assessing how gene expression patterns over an environmental gradient vary within populations (reflecting phenotypic plasticity) and among populations (reflecting genetic differentiation), it is possible to identify candidate genes and regulatory pathways that may be under selection for adaptation (Oleksiak et al. 2002; Whitehead and Crawford 2006a, 2006b). In tetraploid wheat, the QST–FST comparison method has been used to detect selection signatures for metabolites, treated as molecular phenotypic traits (Beleggia et al. 2016). Until now, despite its suitability for the analysis of gene expression data, this method has been rarely adopted in transcriptomics studies (Roberge et al. 2007; Kohn et al. 2008; Aykanat et al. 2011).

We implemented a methodology based on QST distributions and QST–FST comparisons to perform a “selection scan,” seeking genes whose expression was potentially under selection. To ensure the reliability of gene selection and minimize the risk of false positives, we established two thresholds based on the distributions of gene expression heritability (Supplementary Fig. S3), while acknowledging that we developed this methodology as a proof of concept with a limited number of genotypes. Consequently, we aimed to establish a high confidence level to ensure the robustness of our approach. Genes with H2 < 0.7 were removed in order to retain only candidates in the top 15% of the distribution for which we can ascertain that their variation is predominantly genetic. However, to consider those genes whose expression was strongly influenced by N availability, we also evaluated the percentage of the species × environment (S × N) variance component (i.e. every species subgroup × N condition), retaining those genes meeting at least the 20% threshold (Supplementary Fig. S3). Altogether, 5,868 genes (∼18% of the total number) met these criteria and only those with QST values in the 5% right-hand tail of the distributions were considered candidates for selection. The QST and FST values of the filtered genes were then compared (Supplementary Fig. S4) to confirm that their divergent expression (high QST values) was caused by directional selection (QST > FST) and not by genetic drift (QST ≈ FST) or stabilizing selection (QST < FST) (Bonnin et al. 1996; Whitlock and Guillaume 2009). After removing FST values < 0.01, we observed that all the resulting 967 genes satisfied the criterion QST > FST, indicating that their expression was probably subjected to directional selection in at least one of the evolutionary contexts examined herein (i.e. primary and/or secondary domestication under high-N and/or low-N availability conditions). Notably, ∼280 genes were consistently detected in each of the six comparisons (Supplementary Data Set 4). These genes exhibited high QST values, indicating substantial divergence in expression levels, coupled with lower FST values, suggesting limited nucleotide-level divergence. This suggests that the observed signals likely arise from upstream regulatory mechanisms influencing gene expression, rather than mutations within the gene's coding regions, potentially leading to altered gene products. Consequently, we conducted SNP annotation on the 967 genes under selection, accounting for both synonymous and non-synonymous mutations. We then compared these results by annotating variants present in another set of 967 genes randomly selected from the total pool of 32,358 genes. As detailed in Supplementary Table S2, for each taxon, we detected an average lower number of both synonymous and non-synonymous mutations in the genes under selection compared to the randomly selected group. Specifically, the ratios of non-synonymous/synonymous mutations were ∼0.78 for selected genes and ∼0.81 for random ones in wild emmer and emmer wheat, showing comparable values. However, in durum wheat, this ratio was significantly lower: ∼0.69 for selected genes and ∼0.76 for random ones (Kolmogorov–Smirnov two-sided test P < 2.2e-16). These findings bolster our identified signals of selection, suggesting that our 967 candidate genes, precisely because they are putatively under directional selection, may have also undergone a process of purifying selection, thereby preventing the accumulation of deleterious mutations.

Gene Ontology (GO) enrichment analysis was performed on the six groups of genes. Not all six comparisons yielded significant results, but significant and distinct GO categories emerged specifically during primary and secondary domestication under high-N (Supplementary Fig. S5). However, examining the entire evolutionary process (from wild emmer to durum wheat), categories associated with the amino acid biosynthesis were enriched also in genes showing selection signatures when comparing wild emmer and durum wheat under low-N conditions (Supplementary Fig. S5). Our results suggest that primary and secondary domestication involved different natural and artificial selection pressures affecting distinct sets of genes and phenotypes influenced by N availability, indicating that certain genetic pathways may have been particularly important for the adaptation of wheat to different environmental conditions during its domestication history.

Among the genes under selection during primary domestication, categories linked to “defense-related programed cell death, modulated by biotic interactions” were enriched, suggesting an enhanced hypersensitive response to pathogens. As wild genotypes transition to agroecosystems characterized by dense crop monocultures, they encounter increased disease pressure from crop-specific pathogens (Savary et al. 2019). This prompts a hypersensitive response, potentially leading to programmed cell death and necrosis as a defense mechanism. Pathogen defense mechanisms in plants often intersect with the regulation of beneficial symbiotic interactions, suggesting a tradeoff between symbiosis-associated traits and innate immunity (Porter and Sachs 2020). Moreover, domesticated crops are less able to leverage microbial interactions compared to wild counterparts, as shown by comparative studies involving bread wheat landraces and old vs modern varieties (Valente et al. 2023). This reduced capacity may in part reflect the widespread adoption of high-input agricultural practices, wherein the availability of fertilizers diminishes the need for plants to invest in symbiotic relationships (Martín-Robles et al. 2018).

Among the genes under selection during secondary domestication, we observed the enrichment of categories associated with amino acid metabolism, particularly those related to the “lysine catabolic process” (Supplementary Fig. S5). This included genes encoding the bifunctional enzyme lysine ketoglutarate reductase/saccharopine dehydrogenase (LKR/SDH), which breaks down lysine via the saccharopine pathway (SACPATH). The structure and transcription of the LKR/SDH gene have been investigated in durum wheat, revealing species-dependent differences in expression and lineage-specific variations between monocots and dicots (Anderson et al. 2010). Lysine is a limiting essential amino acid in cereal grains, and efforts have been made to enhance its content in crops like maize (Zea mays) and rice (Oryza sativa) by targeting the catabolic pathway (Houmard et al. 2007; Frizzi et al. 2008; Long et al. 2013). However, lysine-rich proteins generally do not accumulate to high levels in cereal seeds, which instead stockpile prolamins (such as gliadin in wheat). The SACPATH appears to direct lysine toward the production of glutamic acid, a precursor of proline, which is abundant in gluten (Arruda and Barreto 2020).

Evolutionary metabolomics has revealed signatures of selection affecting amino acid metabolism during secondary domestication (Beleggia et al. 2016). Changes in amino acid metabolism during domestication have been observed in crops such as sunflower, maize, and common bean based on nucleotide data (Chapman et al. 2008; Swanson-Wagner et al. 2012; Bellucci et al. 2014). In durum wheat, domestication has been associated with the selection of specific protein compositions, reducing the diversity of gliadin and glutenin subunits, thus affecting grain yield and gluten properties (Laidò et al. 2013, 2014). Moreover, the SACPATH is upregulated in response to drought stress in spring wheat genotypes, particularly in drought-tolerant varieties, suggesting a role in stress adaptation (Michaletti et al. 2018). Proline, derived from this pathway, may serve as a major constituent of storage proteins and also a key osmoprotectant produced in response to stress (Kavi Kishor et al. 2022). Overall, selection during wheat domestication may have influenced the expression of SACPATH genes, favoring not only protein composition but also abiotic stress tolerance.

Changes in nitrogen availability trigger gene expression, resulting in a twofold increase in the number of differentially expressed genes in durum wheat compared to emmer and wild emmer wheat

We identified DEGs in each subspecies that discriminated between high-N conditions and N starvation using a stringent pipeline and strict thresholds (Padj < 0.001) to reduce the number of false positives. We found 3,326 DEGs in wild emmer, 3,305 in emmer and 5,901 in durum wheat, with more upregulated than downregulated genes in all three subspecies. Durum wheat had the highest percentage of private DEGs (∼42%, 2,479), whereas similar numbers were found in wild emmer (∼14%, 458) and emmer (∼15%, 486). Wild emmer and emmer shared ∼23% (749) and ∼21% (700), respectively, of their DEGs with durum wheat. The percentage of DEGs shared only between wild emmer and emmer was 4% (146), but almost 60% of wild emmer and emmer DEGs and ∼33% of durum wheat DEGs were shared by all three taxa (Fig. 4A). The proportions of private and shared DEGs were preserved when we separated them into upregulated and downregulated subsets (Fig. 4B and C). In all three taxa, most DEGs were located on chromosomes 2A, 2B, 3A, 3B, 5A, and 5B, each carrying >7.5% of the DEGs, whereas chromosomes 6A and 6B each contained only ∼5% of the DEGs (Supplementary Fig. S6).

Figure 4. Venn diagrams showing differentially expressed genes (DEGs) when comparing high-N and low-N conditions within each subspecies. A) Total set of DEGs, B) upregulated DEGs only, and C) downregulated DEGs only.

GO enrichment analysis of the DEGs meeting the threshold FDR < 0.05 revealed 23 macro-categories in wild emmer, 21 in emmer, and 25 in durum wheat (Supplementary Fig. S7). The main differences between the three subspecies were observed for categories related to “signaling,” “regulation of biological process,” “developmental process,” and “metabolic process” (Supplementary Fig. S7). We observed the uniform enrichment of GO categories associated with upregulated genes in all three subspecies, including terms linked to N and amino acid metabolism as well as carbon (C) metabolism and photosynthesis (Supplementary Fig. S7). In contrast, the enrichment of GO categories associated with downregulated genes was more selective, with some GO categories related to N metabolism enriched only in durum wheat, including GO:0006807 and GO:0034641 (N compound and cellular N compound metabolic process, respectively) and GO:0006536 “glutamate metabolic process” (Supplementary Data Set 5). Functional annotations of the most strongly modulated genes (top 5% |log2FC| values) are reported in Supplementary Data Set 6.

Our data confirm, on a larger set of samples, earlier observations on the response of wheat to N starvation based on transcriptomics and metabolomics data. These earlier studies included one emmer and one durum wheat genotype also present in our sample set (Beleggia et al. 2021), but also considered the durum wheat cultivar Svevo (Curci et al. 2017) and various bread wheat cultivars (Sultana et al. 2020). As expected, genes involved in N metabolism were modulated during N starvation. Among the key genes for N assimilation, those encoding asparagine synthetase and nitrite reductase were upregulated in every taxon, whereas those encoding glutamate carboxypeptidase and glutamate decarboxylase were downregulated. We observed contrasting profiles for genes encoding ureide permease (a ureide transporter), which were strongly upregulated in all three subspecies in response to N stress, whereas genes encoding nitrate transporters were strongly downregulated. The modulated genes also included transporters of amino acids and other nutrients.

N starvation also influenced other metabolic pathways, revealing many further DEGs involved in C metabolism, especially fatty acid metabolism, glycolysis, photosynthesis, and the tricarboxylic acid (TCA) cycle. About 10% of the highest-ranking DEGs represented transcription factors and protein kinases. The most common functional category (accounting for 17% of annotated DEGs) reflected the general stress response to N starvation, including the mitigation of oxidative stress and detoxification. Examples included genes encoding cytochrome P450s, glutaredoxin family proteins, glutathione S-transferases, and peroxidases (Supplementary Data Set 6).

To compare gene expression between the three taxa while taking the environmental effects into account, we also identified DEGs between each pair of subspecies under all N conditions. Accordingly, we compared emmer vs wild emmer (primary domestication, high- and low-N), durum wheat vs emmer (secondary domestication, high- and low-N), and durum wheat vs wild emmer (cumulative effect, high- and low-N) (Supplementary Fig. S8). The wild emmer vs emmer comparison revealed few DEGs regardless of N availability (12 and 11 DEGs under high-N and low-N conditions, respectively), whereas the emmer vs durum wheat comparison revealed 41 DEGs associated with high N and 29 associated with N starvation, and the wild emmer vs durum wheat comparison revealed 46 DEGs associated with high-N and only 10 associated with N starvation. These data indicate that the number of DEGs increases during domestication but only when there is a sufficient N supply (Supplementary Fig. S8). Interestingly, there were more upregulated than downregulated genes in all pairwise comparisons under high-N conditions (∼65%) but the proportion increased under N starvation, particularly for the comparison of wild emmer vs durum wheat (90%). The preponderance of upregulated genes during domestication has also been observed in maize (Lemmon et al. 2014), whereas domestication was shown to increase the proportion of downregulated genes in common bean (Bellucci et al. 2014), eggplant (Solanum melongena) (Page et al. 2019), and sorghum (Burgarella et al. 2021) landraces compared to wild relatives. The absence of consistent patterns suggests that the evolution of domesticated phenotypes is driven by specific processes that are unique to each crop.

Among the 102 DEGs found in at least one of the six pairwise comparisons between subspecies (Supplementary Data Set 7), 35 were also found among DEGs identified between contrasting N conditions and 24 of these were proposed to be under selection. Overall, six genes were identified in all three experiments (i.e. differentially expressed between subspecies and between contrasting N conditions, and showed evidence of selection).

Selection shaped the expression profiles of genes modulated by nitrogen availability

The 6,991 DEGs found in at least one species when comparing contrasting N conditions included 101 putatively under selection, which are candidates for the adaptive response to N availability. The expression profiles of these selected genes, in all three subspecies under both N conditions, are shown in Supplementary Fig. S9. We applied PCA to the normalized read counts in order to determine whether the different genotype groups can be separated based on their gene expression. Initially, we incorporated all 6,991 DEGs (Fig. 5A and B) before focusing on the subset of 101 DEGs that were also putatively under selection (Fig. 5C and D). When considering all DEGs, PC1 did not completely separate the durum wheat genotypes from the other taxa, in contrast to the clear separation observed for the SNP data (Fig. 1), and this was particularly evident during N starvation (Fig. 5A). There was also a moderate degree of overlap between the wild emmer and emmer genotypes along PC2. However, when we focused on the DEGs under selection, PC1 separated the durum wheat genotypes into a densely clustered group (as observed for the SNP data) under both N conditions, and PC2 separated the wild emmer and emmer genotypes more clearly, especially under high-N conditions (Fig. 5C and D). For an overview of overall gene expression patterns in the two N conditions, please refer to Supplementary Fig. S10, which illustrates the PCA conducted on the complete set of genes (32,358).

Figure 5. Principal component analysis (PCA) of 32 wheat genotypes based on expression data of differentially expressed genes (DEGs) when comparing high-N and low-N conditions within each subspecies. A, B) Plots based on all 6,991 DEGs (not filtered): A low-N conditions and B high-N conditions. C, D) Plots based on 101 DEGs that are also putatively under selection: C low-N conditions and D high-N conditions. Samples are represented by taxon-specific colored dots. Labels show the accession name of each genotype.

The selection signatures (based on QST–FST values) allowed us to identify genes that have putatively diverged between tetraploid wheat subspecies due to selection pressure. By integrating this information with classical differential expression analysis, which identified genes that are differentially expressed under contrasting N conditions, it was possible to pinpoint 101 candidate genes that are both genetically differentiated and functionally relevant to the environmental factor of interest: N availability during the domestication and diversification of cultivated wheat. Functional annotation (Supplementary Data Set 8) revealed upregulated genes associated with C metabolism as well as some encoding transcription factors and transporters, as well as both upregulated and downregulated genes associated with general stress responses and N metabolism, specifically those encoding enzymes involved in amino acid metabolism such as methionine aminopeptidase, aspartokinase and glutamate dehydrogenase (GDH). The latter is particularly noteworthy because, in addition to its modulation in response to different N conditions and the presence of selection signatures, the GDH gene was also upregulated in the comparison between wild emmer and durum wheat under high-N conditions (Supplementary Data Set 7). GDH is a key enzyme involved in N metabolism and N/C balance (Miflin and Habash 2002). This is supported by the co-localization of quantitative trait loci for GDH activity and physiological traits associated with the flag leaf lamina, such as soluble protein and amino acid content, as well as flag leaf area and dry weight (Fontaine et al. 2009). Selection signatures were also identified in the GDH gene when comparing landraces with old and modern durum wheat cultivars (Taranto et al. 2020). Our results confirm that N metabolism has been a key driver during the evolutionary history of wheat, particularly the central role of glutamate in the process of domestication. This was also suggested by a combined transcriptomics and metabolomics study showing that glutamate and γ‐aminobutyric acid (mainly synthesized from glutamate) are central to the genotype-specific response of emmer and durum wheat to N starvation (Beleggia et al. 2021).

We have shown that significant changes occurred at the nucleotide and gene expression levels during the domestication of tetraploid wheat, taking into account the environmental variable of N availability. We observed that more nucleotide diversity has been lost during secondary domestication compared to primary domestication. In addition, we found a parallel trend in the loss of gene expression diversity associated with the domestication process, with a stronger effect due to secondary domestication. The observed loss of expression diversity may be related to N availability in the durum wheat selection environment. Our findings suggest that selection may have operated in different directions during primary and secondary domestication, the former involving changes related to biotic interactions and the latter related to amino acid metabolism.

Despite the limited number of genotypes available for our study, the innovative combination of RNA-Seq analysis and the estimation of quantitative genetic parameters allowed us to develop a pipeline for the identification of selection signatures and phenotypic plasticity in gene expression data based on evolvability and QST–FST scores. Emphasizing the pioneering nature of our work, we opted to introduce stringent and high confidence thresholds for the considered parameters, aiming to present our methodology as a proof of concept. While presenting promising results, we acknowledge the potential for further refinement and adjustment of the methodology in experiments employing larger genotype datasets. The set of genes we identified with underlying selection signatures will facilitate the development of innovative strategies to improve resource use efficiency and environmental sustainability in crop management.

Materials and methods

Plant material and experimental design

The study included 32 tetraploid wheat genotypes, comprising 10 accessions of wild emmer (T. turgidum ssp. dicoccoides), 10 accessions of emmer (T. turgidum ssp. dicoccum), and 12 accessions of durum wheat (T. turgidum ssp. durum) (Supplementary Data Set 1). The samples analyzed in our study were selected from a larger experiment conducted in October 2012, as previously described (Gioia et al. 2015) and were chosen as representative of the majority of the diversity within the panel. Briefly, wheat genotypes were grown for 4 wk under high-N and nitrogen starvation (low-N) conditions in the Phytec Experimental Greenhouse at the Institute of Biosciences and Geosciences (IBG-2), Plant Sciences Institute, Forschungszentrum Jülich GmbH, Germany (50°54′36′′ N, 06°24′49′′ E). Seeds of uniform size and mass were visually selected, surface sterilized with 1% (w/v) NaClO for 15 min and pre-germinated. After germination, seedlings showing uniform growth (seminal root length, 1 to 2 cm) were transferred to soil-filled rhizoboxes, which were placed in the automated GROWSCREEN-Rhizo phenotyping system available at IBG-2. We used a Type 0 manually sieved peat soil (Nullerde Einheitserde; Balster Einheitserdewerk, Frondenberg, Germany), which provided low nutrient availability (ammonium N and nitrate N concentrations of <1.0 and <1.0 mg l−1, respectively). All plants were watered twice daily with 400 ml tap water and were supplied three times per week with 200 ml modified Hoagland solution (Hoagland and Arnon 1950), adapted for optimal N and N starvation conditions. Stock solution contained 5 mm KNO3, 5 mm Ca(NO3)2, 2 mm MgSO4, 1 mm KH2PO4, and trace elements. For N starvation conditions, KNO3 and Ca(NO3)2 were replaced with K2SO4 and CaCl2·6(H2O), respectively. The experiment was carried out under natural lighting in the greenhouse, with an air temperature of 18 to 24 °C and a relative humidity of 40% to 60%. For each N treatment, we used two replicates of each genotype with two plants per replicate (four plants per genotype in total). After 4 wk, leaves were pooled from two plants of the same genotype growing in the same rhizobox. Accordingly, four independent biological replicates (two replicates per N condition) were produced for each genotype, with the exception of wild emmer IG 46504, PI 233288, PI 466991, PI 538656, emmer MG 5293/1, and durum wheat Creso, Pedroso, and Trinakria, for which only three replicates were available, and emmer Molise Sel. Colli and durum wheat Simeto, for which eight replicates were available. The tissues were immediately frozen in liquid N2 and stored at −80 °C. Further details of the experiment and growth conditions are provided elsewhere (Gioia et al. 2015).

RNA extraction and sequencing

RNA was extracted from 100 mg of frozen ground leaves per replicate using the Spectrum Plant Total RNA kit (Sigma-Aldrich, St Louis, MO, USA) followed by treatment with RNase-free DNase using the On-Column DNase I Digestion Set (Sigma-Aldrich). RNA integrity and purity were assessed by agarose gel electrophoresis and a Bioanalyzer 2100, respectively (Agilent/Bonsai Technologies, Santa Clara, CA, USA). Only RNA samples with an RNA integrity number >8.0 were considered suitable for analysis.

Library construction and RNA sequencing were carried out using the Illumina mRNA-Seq platform at the Montpellier Genomix sequencing facility (http://www.mgx.cnrs.fr) as previously described (David et al. 2014). Briefly, RNA samples were processed using TruSeq RNA sample preparation kits v2 (Illumina, San Diego, CA, USA). Libraries were quantified by RT-qPCR using the KAPA Library Quantification Kit for Illumina Sequencing Platforms (Roche, Basel, Switzerland), followed by quality control using a DNA 100 Chip on a Bioanalyzer 2100. Cluster generation and sequencing were carried out using the Illumina HiSeq 2000 instrument and TruSeq PE Cluster Kit v3, following the Illumina PE_Amp_Lin_Block_V8.0 recipe, and Illumina TruSeq PE Cluster v3-cBot-HS kits with the 2 × 100 cycles, paired-end, indexed protocol, respectively (David et al. 2014).

RNA-Seq library processing and mapping

We pre-processed 128 raw paired-end RNA-Seq libraries (David et al. 2014). Cutadapt (Martin 2011) was then used to remove adaptor sequences and trim the end of reads with low-quality scores (parameter -q 20) while keeping reads with a minimum length of 35 bp. Reads with a mean quality score < 30 were discarded, and orphan reads (whose mates were discarded in the previous filtering steps) were also removed (David et al. 2014). The final quality of trimmed and filtered reads was assessed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).

The bread wheat (Triticum aestivum cv. Chinese Spring) genome assembly IWGSC RefSeq v2.1, along with the corresponding genome annotation, were downloaded from the IWGSC data repository hosted by URGI-INRAE (https://wheat-urgi.versailles.inra.fr/) and used as a reference to map each cleaned library to the A and B subgenomes. To validate our choice of reference genome, reads were also mapped to the available tetraploid wheat reference genomes: wild emmer accession Zavitan (https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_002162155.3/) and durum wheat cultivar Svevo (https://www.interomics.eu/durum-wheat-genome). Pairwise genetic distances between the three reference genomes were computed using Mash v2.3 (Ondov et al. 2016).

STAR v2.7.0e (Dobin et al. 2013) was used for read mapping with the –quantMode TranscriptomeSAM and –quantTranscriptomeBan Singleend options. The output alignments were translated into transcript coordinates (in addition to alignments in genomic coordinates), allowing insertions, deletions, and soft-clips in the transcriptomic alignments. The transcriptomic alignments were used as inputs for salmon v1.6.0 (Patro et al. 2017) to quantify gene expression. Raw read counts were computed for all genes in each sample and, to filter out weakly expressed transcripts, only genes with at least 1 count per million (CPM) in at least 10 samples (of the same subspecies) were retained. This was calculated separately in each of the three subspecies and the raw counts of the filtered genes in each subspecies were then combined for downstream analysis (Supplementary Data Set 3).

Variant identification

Variants were called by applying BCFtools v1.15 (previously SAMtools) (Danecek et al. 2021) to the alignment bam files. The “bcftools mpileup” command was used to determine the genotype likelihoods at each genomic position, with a minimum alignment quality of 20 and a minimum base quality of 30. The actual calls were obtained using the “bcftools call” command. The resulting VCF file was filtered using the “bcftools view” command, removing indels and keeping only sites covered by at least three reads in all genotypes. Subsequently, only biallelic SNPs with maximum values of 50% missingness and a 1% minor allele frequency were retained. To identify private and shared SNPs among the different subspecies, every possible comparison of the three subsampled VCF files (wild emmer, emmer, and durum wheat) was carried out using the “bcftools isec” command.

Population genetics analysis

Variants were filtered (one SNP per 500 kb) using the VCFtools v0.1.17 –thin 500000 option (Danecek et al. 2011) and then converted to ped format with PLINK v1.90p (Purcell et al. 2007). PLINK was also used to compute genetic distances between individuals with the –distance-matrix flag. The output matrix was used as input for PCA with the cmdscale function of R v4.2.1 (R Core Team 2022).

Genetic diversity statistics, including nucleotide diversity (π and θ) (Watterson 1975; Tajima 1983) were computed on the alignment bam files for each subspecies, from the folded site frequency spectra using ANGSD (Korneliussen et al. 2014). First, the doSaf function was used to estimate per-site allele frequencies (Saf); then, realSFS was used to get the site frequency spectra. The loss of diversity statistic (Vigouroux et al. 2002) was used to test the impact of primary and secondary domestication on the molecular diversity of the three subspecies. For primary domestication, the statistic was computed as [1−(xemmer/xwild)], where xemmer and xwild are the diversities in emmer and wild emmer, respectively, measured using π and θ. If xemmer was higher than xwild, then the parameter was calculated as [(xwild/xemmer)–1]. The loss of diversity due to secondary domestication in durum wheat versus emmer was calculated as [1−(xdurum/xemmer)], where xdurum and xemmer are the diversities in durum wheat and emmer, respectively. If xdurum was higher than xemmer, then the parameter was calculated as [(xemmer/xdurum)–1].

We calculated FST for each pair of populations using ANGSD (Korneliussen et al. 2014). Saf and 2D site frequency spectra were calculated for nucleotide diversity, and then, the fst index function was used to obtain the global estimate. To get an FST value for each gene in our dataset, we used the fst print function, which prints the posterior expectation of genetic variance between populations (called A), and total expected variance (called B) for every locus. We then computed the weighted FST as the ratio between the summed As and summed Bs for every gene region, using an ad hoc R script.

Expression profiles, heritability, and QST analysis

Raw read counts of the 32,358 genes were normalized using the vst method allowing the additive CVA (standard deviation/mean) to be calculated for the two N conditions in every subspecies, averaging the biological replicates of every genotype. The statistical loss approach (Vigouroux et al. 2002) was then applied to test the loss of expression diversity in the different groups, as previously reported (Bellucci et al. 2014). The statistical significance of the differences between each CVA value and the percentage loss of expression diversity was determined using the Mann–Whitney U-test in R v4.2.1 (R Core Team 2022) with the function wilcox.test.

To compute heritability, the raw counts of each subspecies under each condition were first normalized using the trimmed mean M-values normalization method in the R package edgeR (Robinson et al. 2010) and the voom normalization method in the R package limma (Smyth 2005). To determine the variance component of each factor and heritability, the following model was considered:

Yijkl=Si+Gj(i)+Nk+(S×N)ik+(G×N)jk(i)+εl(ijk)

where Yijkl is the normalized gene expression level, Si is the species factor, Gj(i) is the genotype factor nested in species, Nk is the N-level factor, (S×N)ik is the interaction between species and N levels, (G×N)jk(i) is the interaction between genotypes and N levels, and εl(ijk) is the residual error. All factors were treated as random effects in the model except the intercept, which was a fixed effect. The linear mixed models were fitted using the lmer function in R package lme4 based on the normalized data of each transcript (Bates et al. 2015). The heritability (H2) was calculated as H2=VS+VGVA, where VA=VS+VG+VN+VS×Nn+VG×Nn+Vεn, VS is the variance of species, VG is the variance of genotype, VN is the variance of N level, VS×N is the variance of species and N-level interaction, VG×N is the variance of genotype and N-level interaction, Vε is the residual variance, and n is the number of N levels. VS×N and VG×N represent the genotype × environment interaction variance components at the species and genotype (nested in species) levels, respectively.

Q ST was calculated between pairs of the three subspecies under low-N and high-N conditions separately. The wild emmer vs emmer comparison revealed the effects of primary domestication, the emmer vs durum wheat comparison revealed the effects of secondary domestication, and the wild emmer vs durum wheat comparison revealed the cumulative effect of domestication. To this end, the model can be reduced to Yijl=Si+Gj(i)+εl(ij) at each N level. The QST value was calculated as QST=VSVS+VG, that is, the ratio of between-species to within-species variances.

Q ST Distributions were used to perform a “selection scan” on a restricted number of genes. First, genes were filtered for H2 ≥ 0.7, and in order not to lose genes whose expression was strongly influenced by N availability, the species × environment (S × N) variance component was also evaluated (i.e. every species subgroup × N condition), retaining those genes meeting the threshold S × N ≥ 0.2 (Supplementary Fig. S3). Successively, we obtained six different QST value distributions (QST WILD EMMER VS EMMER, QST EMMER VS DURUM WHEAT, and QST WILD EMMER VS DURUM WHEAT, each for high-N and low-N conditions), and we retained the 5% right-hand tail of every distribution. Finally, we compared FST and QST values for every gene, discarding FST values < 0.01. We confirmed that every retained gene satisfied the condition QST > FST, allowing it to be classed as undergoing directional selection. For these selected genes, SNPs were annotated using SnpEff v.5.1d (Cingolani et al. 2012).

Differential expression analysis

Differential gene expression was assessed by analyzing the pre-processed raw count dataset (32,358 genes). We identified DEGs by comparing (i) two conditions (i.e. high-N and low-N levels) within each subspecies, and (ii) pairs of the three subspecies under the same N levels, which considered the genotypes nested in species. For the two scenarios, we used three different approaches to detect DEGs: one linear model-based approach implemented in the R package limma (Smyth 2005), and two Poisson model-based approaches implemented in the R packages edgeR (Robinson et al. 2010) and DESeq2 (Love et al. 2014). In all approaches, the normalization of raw counts was applied by default in the package before differential analysis. To reduce the number of false positives, the intersection of DEGs resulting from the three approaches was retained (Rapaport et al. 2013; Zhang et al. 2014), and the significance threshold was set to an adjusted P-value < 0.001. The DEGs between high- and low-N levels in at least one subspecies were used for PCA following the DESeq2 approach (Love et al. 2014), first using all the DEGs, then repeating the analysis on the DEGs considered to be under selection. At each step, counts were normalized using the vst method before the plotPCA function was applied to define principal components 1 and 2 for the two N levels separately. The expression patterns of the DEGs considered to be under selection were plotted from normalized read counts using the R package pheatmap (https://github.com/raivokolde/pheatmap).

GO enrichment analysis

Enriched terms in the DEGs and genes under selection were identified using agriGO v.2.0 (Tian et al. 2017). All annotated genes of bread wheat were used as background, and the following parameters were set: hypergeometric test, multiple hypothesis test adjustment according to the Hochberg FDR procedure at significance level < 0.05, and minimum number of mapping entries = 3.

Accession numbers

The RNA-Seq libraries generated and analyzed in this study have been deposited in the Sequence Read Archive (SRA) of the National Center of Biotechnology Information (NCBI) under BioProject number PRJNA1015013. The nucleotide sequences of all discussed genes are available on the GarinGenes Database (https://wheat.pw.usda.gov/cgi-bin/GG3/browse.cgi) and can be accessed using the Gene model IDs provided in the Supplementary Data Sets (e.g. TraesCS2A03G0941300).

Supplementary Material

koae202_Supplementary_Data

Acknowledgments

We acknowledge the CINECA award under the ISCRA initiative, for the availability of high-performance computing resources and support.

Author contributions

R.P. and R.B. conceived and designed the study. T.G. and F.F. carried out the experiments. J.L.D. carried out the library preparation and sequencing. A.P. performed the RNA-Seq analysis. A.P., H.T., and Z.N. performed the bioinformatics analysis and analyzed the data. C.D.Q., A.R.L., and M.R. provided technical support for RNA-Seq analysis. A. P., R.B., and R.P. wrote the paper. Z.N., J.L.D., U.S., V.D.V., G.F., E.Bi., L.N., E.Be., N.P., S.M., and P.D.V. reviewed and contributed to the editing of the manuscript. All authors have read and approved the manuscript.

Supplementary data

The following materials are available in the online version of this article:

Supplementary Figure S1. Folded site frequency spectra (SFS) of single-nucleotide polymorphisms (SNPs) in the three wheat taxa.

Supplementary Figure S2. Principal component analysis (PCA) of 32 wheat genotypes based on single-nucleotide polymorphisms (SNPs) using different reference genomes.

Supplementary Figure S3. Workflow of gene expression selection scanning.

Supplementary Figure S4. Comparison of QST and FST estimates of the 5,868 genes showing H2 ≥ 0.7 or S × N ≥ 0.2.

Supplementary Figure S5. Gene ontology (GO) categories of genes under selection.

Supplementary Figure S6. Genome-wide distribution of differentially expressed genes (DEGs) in the comparison between contrasting N conditions within each subspecies.

Supplementary Figure S7. GO Classification of DEGs in the comparison between contrasting N conditions within each subspecies.

Supplementary Figure S8. DEGs between subspecies.

Supplementary Figure S9. Expression profiles of the 101 DEGs putatively under selection in the three wheat taxa in high-N (N+) and low-N (N−) conditions.

Supplementary Figure S10. Principal component analysis (PCA) of 32 wheat genotypes based on expression data of all 32,358 genes in each subspecies.

Supplementary Table S1. Mean CVA in gene expression for the three wheat taxa and loss of expression diversity for two gene subgroups (6,991 differentially expressed and 25,367 non-differentially expressed genes).

Supplementary Table S2. Synonymous and non-synonymous mutations and non-synonymous/synonymous ratios of two gene groups: 967 genes under selection and 967 randomly selected genes.

Supplementary Data Set 1, A to C. List of the 128 samples and reads mapping results using three different reference genomes.

Supplementary Data Set 2. Results of single-nucleotide polymorphism (SNP) calling and nucleotide diversity estimates, using three different reference genomes.

Supplementary Data Set 3. Raw read counts of the 32,358 genes in the 128 RNA-Seq samples.

Supplementary Data Set 4. List of the 967 genes retained from the “selection scan”.

Supplementary Data Set 5, A to D. List of GO “Biological process” and “Molecular function” subcategories for upregulated and downregulated DEGs under different nitrogen N conditions for each subspecies.

Supplementary Data Set 6. Functional annotations of the DEGs between N conditions in each subspecies. Genes with the top 5% |log2FC| values are shown.

Supplementary Data Set 7. Functional annotations of the DEGs between subspecies under all N conditions.

Supplementary Data Set 8. Functional annotation of the 101 genes selected by the integration of selection signatures and differential expression analysis between N conditions.

Funding

This work was financially supported by Università Politecnica delle Marche, by the PON a3 PlASS Project (PlASS - Platform for Agrofood Science and Safety) and the Transnational Access capacities of the European Plant Phenotyping Network (EPPN, grant no. 284443) funded by the FP7 Research Infrastructures. H.T. and Z.N. were funded by the Horizon 2020 Research and Innovation Framework Programme through grant no. 862201.

Data availability

Scripts used in this study are available at Github https://github.com/PapaLab/transcriptomics_triticum.git.
==== Refs
References

Anderson OD , Coleman-DerrD, GuYQ, HeathS. Structural and transcriptional analysis of plant genes encoding the bifunctional lysine ketoglutarate reductase saccharopine dehydrogenase enzyme. BMC Plant Biol. 2010:10 (1 ):113. 10.1186/1471-2229-10-113 20565711
Arruda P , BarretoP. Lysine catabolism through the saccharopine pathway: enzymes and intermediates involved in plant responses to abiotic and biotic stress. Front Plant Sci. 2020:11 :587. 10.3389/fpls.2020.00587 32508857
Aykanat T , ThrowerFP, HeathDD. Rapid evolution of osmoregulatory function by modification of gene transcription in steelhead trout. Genetica. 2011:139 (2 ):233–242. 10.1007/s10709-010-9540-2 21190065
Barneix AJ . Physiology and biochemistry of source-regulated protein accumulation in the wheat grain. J Plant Physiol. 2007:164 (5 ):581–590. 10.1016/j.jplph.2006.03.009 16690166
Bates D , MächlerM, BolkerB, WalkerS. Fitting linear mixed-effects models using lme4. J Stat Softw. 2015:67 (1 ):1–48. 10.18637/jss.v067.i01
Beleggia R , OmranianN, HoltzY, GioiaT, FioraniF, NigroFM, PecchioniN, De VitaP, SchurrU, DavidJL, et al Comparative analysis based on transcriptomics and metabolomics data reveal differences between emmer and durum wheat in response to nitrogen starvation. Int J Mol Sci. 2021:22 (9 ):4790. 10.3390/ijms22094790 33946478
Beleggia R , RauD, LaidòG, PlataniC, NigroF, FragassoM, De VitaP, ScossaF, FernieAR, NikoloskiZ, et al Evolutionary metabolomics reveals domestication-associated changes in tetraploid wheat kernels. Mol Biol Evol. 2016:33 (7 ):1740–1753. 10.1093/molbev/msw050 27189559
Bellucci E , BitocchiE, FerrariniA, BenazzoA, BiagettiE, KlieS, MinioA, RauD, RodriguezM, PanzieraA. Decreased nucleotide and expression diversity and modified coexpression patterns characterize domestication in the common bean. Plant Cell. 2014:26 (5 ):1901–1912. 10.1105/tpc.114.124040 24850850
Bitocchi E , RauD, BellucciE, RodriguezM, MurgiaML, GioiaT, SantoD, NanniL, AtteneG, PapaR. Beans (Phaseolus ssp.) as a model for understanding crop evolution. Front Plant Sci. 2017:8 :722. 10.3389/fpls.2017.00722 28533789
Bogaard A , FraserR, HeatonTH, WallaceM, VaiglovaP, CharlesM, JonesG, EvershedRP, StyringAK, AndersenNH. Crop manuring and intensive land management by Europe's first farmers. Proc Natl Acad Sci U S A. 2013:110 (31 ):12589–12594. 10.1073/pnas.1305918110 23858458
Bonnin I , ProsperiJM, OlivieriI. Genetic markers and quantitative genetic variation in Medicago truncatula (Leguminosae): a comparative analysis of population structure. Genetics. 1996:143 (4 ):1795–1805. 10.1093/genetics/143.4.1795 8844165
Bradshaw AD . Evolutionary significance of phenotypic plasticity in plants. Adv Genet. 1965:13 :115–155. 10.1016/S0065-2660(08)60048-6
Burgarella C , BergerA, GléminS, DavidJ, TerrierN, DeuM, PotD. The road to Sorghum domestication: evidence from nucleotide diversity and gene expression patterns. Front Plant Sci. 2021:12 :666075. 10.3389/fpls.2021.666075 34527004
Chapman MA , PashleyCH, WenzlerJ, HvalaJ, TangS, KnappSJ, BurkeJM. A genomic scan for selection reveals candidates for genes involved in the evolution of cultivated sunflower (Helianthus annuus). Plant Cell. 2008:20 (11 ):2931–2945. 10.1105/tpc.108.059808 19017747
Cingolani P , PlattsA, WangLL, CoonM, NguyenT, WangL, LandSJ, LuX, RudenDM. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin). 2012:6 (2 ):80–92. 10.4161/fly.19695 22728672
Civáň P , IvaničováZ, BrownTA. Reticulated origin of domesticated emmer wheat supports a dynamic model for the emergence of agriculture in the fertile crescent. PLoS One. 2013:8 (11 ):e81955. 10.1371/journal.pone.0081955 24312385
Curci PL , Aiese CiglianoR, ZuluagaDL, JanniM, SanseverinoW, SonnanteG. Transcriptomic response of durum wheat to nitrogen starvation. Sci Rep. 2017:7 (1 ):1176. 10.1038/s41598-017-01377-0 28446759
Danecek P , AutonA, AbecasisG, AlbersCA, BanksE, DePristoMA, HandsakerRE, LunterG, MarthGT, SherryST, et al The variant call format and VCFtools. Bioinformatics. 2011:27 (15 ):2156–2158. 10.1093/bioinformatics/btr330 21653522
Danecek P , BonfieldJK, LiddleJ, MarshallJ, OhanV, PollardMO, WhitwhamA, KeaneT, McCarthySA, DaviesRM, et al Twelve years of SAMtools and BCFtools. Gigascience. 2021:10 (2 ):giab008. 10.1093/gigascience/giab008 33590861
David J , HoltzY, RanwezV, SantoniS, SarahG, ArdissonM, PouxG, ChouletF, GenthonC, RoumetP, et al Genotyping by sequencing transcriptomes in an evolutionary pre-breeding durum wheat population. Mol Breed. 2014:34 (4 ):1531–1548. 10.1007/s11032-014-0179-z
Dobin A , DavisCA, SchlesingerF, DrenkowJ, ZaleskiC, JhaS, BatutP, ChaissonM, GingerasTR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013:29 (1 ):15–21. 10.1093/bioinformatics/bts635 23104886
Donald CM . The breeding of crop ideotypes. Euphytica. 1968:17 (3 ):385–403. 10.1007/BF00056241
Fontaine JX , RavelC, PageauK, HeumezE, DuboisF, HirelB, Le GouisJ. A quantitative genetic study for elucidating the contribution of glutamine synthetase, glutamate dehydrogenase and other nitrogen-related physiological traits to the agronomic performance of common wheat. Theor Appl Genet. 2009:119 (4 ):645–662. 10.1007/s00122-009-1076-4 19513687
Fréville H , MontazeaudG, ForstE, DavidJ, PapaR, TenaillonMI. Shift in beneficial interactions during crop evolution. Evol Appl. 2022:15 (6 ):905–918. 10.1111/eva.13390 35782010
Frizzi A , HuangS, GilbertsonLA, ArmstrongTA, LuethyMH, MalvarTM. Modifying lysine biosynthesis and catabolism in corn with a single bifunctional expression/silencing transgene cassette. Plant Biotechnol J. 2008:6 (1 ):13–21. 10.1111/j.1467-7652.2007.00290.x 17725550
Garcia-Gonzalez F , SimmonsLW, TomkinsJL, KotiahoJS, EvansJP. Comparing evolvabilities: common errors surrounding the calculation and use of coefficients of additive genetic variation. Evolution (NY). 2012:66 :2341–2349. 10.1111/j.1558-5646.2011.01565.x
Gioia T , NagelKA, BeleggiaR, FragassoM, FiccoDBM, PieruschkaR, De VitaP, FioraniF, PapaR. Impact of domestication on the phenotypic architecture of durum wheat under contrasting nitrogen fertilization. J Exp Bot. 2015:66 (18 ):5519–5530. 10.1093/jxb/erv289 26071535
Hansen TF , PélabonC, HouleD. Heritability is not evolvability. Evol Biol. 2011:38 (3 ):258–277. 10.1007/s11692-011-9127-6
Haudry A , CenciA, RavelC, BataillonT, BrunelD, PoncetC, HochuI, PoirierS, SantoniS, GléminS, et al Grinding up wheat: a massive loss of nucleotide diversity since domestication. Mol Biol Evol. 2007:24 (7 ):1506–1517. 10.1093/molbev/msm077 17443011
Hawkesford MJ , GriffithsS. Exploiting genetic variation in nitrogen use efficiency for cereal crop improvement. Curr Opin Plant Biol. 2019:49 :35–42. 10.1016/j.pbi.2019.05.003 31176099
Hoagland DR , ArnonDI. The water-culture method for growing plants without soil. Circ Calif Agric Exp Stn. 1950:347 :32.
Houle D . Comparing evolvability and variability of quantitative traits. Genetics. 1992:130 (1 ):195–204. 10.1093/genetics/130.1.195 1732160
Houmard NM , MainvilleJL, BoninCP, HuangS, LuethyMH, MalvarTM. High-lysine corn generated by endosperm-specific suppression of lysine catabolism using RNAi. Plant Biotechnol J. 2007:5 (5 ):605–614. 10.1111/j.1467-7652.2007.00265.x 17553105
Howarth JR , ParmarS, JonesJ, ShepherdCE, CorolDI, GalsterAM, HawkinsND, MillerSJ, BakerJM, VerrierPJ. Co-ordinated expression of amino acid metabolism in response to N and S deficiency during wheat grain filling. J Exp Bot. 2008:59 (13 ):3675–3689. 10.1093/jxb/ern218 18791197
Kavi Kishor PB , SuravajhalaP, RathnagiriP, SreenivasuluN. Intriguing role of proline in redox potential conferring high temperature stress tolerance. Front Plant Sci. 2022:13 :867531. 10.3389/fpls.2022.867531 35795343
Kohn MH , ShapiroJ, WuCI. Decoupled differentiation of gene expression and coding sequence among Drosophila populations. Genes Genet Syst. 2008:83 (3 ):265–273. 10.1266/ggs.83.265 18670138
Korneliussen TS , AlbrechtsenA, NielsenR. ANGSD: analysis of next generation sequencing data. BMC Bioinformatics. 2014:15 (1 ):1–13. 10.1186/s12859-014-0356-4 24383880
Laidò G , ManginiG, TarantoF, GadaletaA, BlancoA, CattivelliL, MaroneD, MastrangeloAM, PapaR, De VitaP. Genetic diversity and population structure of tetraploid wheats (Triticum turgidum L.) estimated by SSR, DArT and pedigree data. PLoS One. 2013:8 (6 ):e67280. 10.1371/journal.pone.0067280 23826256
Laidò G , MaroneD, RussoMA, ColecchiaSA, MastrangeloAM, De VitaP, PapaR. Linkage disequilibrium and genome-wide association mapping in tetraploid wheat (Triticum turgidum L.). PLoS One. 2014:9 (4 ):e95211. 10.1371/journal.pone.0095211 24759998
Laitinen RAE , NikoloskiZ. Genetic basis of plasticity in plants. J Exp Bot. 2019:70 (3 ):795–804. 10.1093/jxb/ery404 30462241
Leinonen T , McCairnsRJS, O’HaraRB, MeriläJ. QST-FST comparisons: evolutionary and ecological insights from genomic heterogeneity. Nat Rev Genet. 2013:14 (3 ):179–190. 10.1038/nrg3395 23381120
Lemmon ZH , BukowskiR, SunQ, DoebleyJF. The role of cis regulatory evolution in maize domestication. PLOS Genet. 2014:10 (11 ):e1004745. 10.1371/journal.pgen.1004745 25375861
Levy AA , FeldmanM. Evolution and origin of bread wheat. Plant Cell. 2022:34 (7 ):2549–2567. 10.1093/plcell/koac130 35512194
Liu W , ChenL, ZhangS, HuF, WangZ, LyuJ, WangB, XiangH, ZhaoR, TianZ, et al Decrease of gene expression diversity during domestication of animals and plants. BMC Evol Biol. 2019:19 (1 ):1–11. 10.1186/s12862-018-1340-9 30606099
Long X , LiuQ, ChanM, WangQ, SunSSM. Metabolic engineering and profiling of rice with increased lysine. Plant Biotechnol J. 2013:11 (4 ):490–501. 10.1111/pbi.12037 23279104
Love MI , HuberW, AndersS. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014:15 (12 ):550. 10.1186/s13059-014-0550-8 25516281
Ludemann CI , GruereA, HefferP, DobermannA. Global data on fertilizer use by crop and by country. Sci Data. 2022:9 (1 ):1–8. 10.1038/s41597-022-01592-z 35013360
Luo MC , YangZL, YouFM, KawaharaT, WainesJG, DvorakJ. The structure of wild and domesticated emmer wheat populations, gene flow between them, and the site of emmer domestication. Theor Appl Genet. 2007:114 (6 ):947–959. 10.1007/s00122-006-0474-0 17318496
Maccaferri M , HarrisNS, TwardziokSO, PasamRK, GundlachH, SpannaglM, OrmanbekovaD, LuxT, PradeVM, MilnerSG. Durum wheat genome highlights past domestication signatures and future improvement targets. Nat Genet. 2019:51 :885–895. 10.1038/s41588-019-0381-3 30962619
Martin M . Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011:17 (1 ):10. 10.14806/ej.17.1.200
Martín-Robles N , LehmannA, SecoE, ArocaR, RilligMC, MillaR. Impacts of domestication on the arbuscular mycorrhizal symbiosis of 27 crop species. New Phytol. 2018:218 (1 ):322–334. 10.1111/nph.14962 29281758
Michaletti A , NaghaviMR, ToorchiM, ZollaL, RinalducciS. Metabolomics and proteomics reveal drought-stress responses of leaf tissues from spring-wheat. Sci Rep. 2018:8 :1–18. 10.1038/s41598-018-24012-y 29311619
Miflin BJ , HabashDZ. The role of glutamine synthetase and glutamate dehydrogenase in nitrogen assimilation and possibilities for improvement in the nitrogen utilization of crops. J Exp Bot. 2002:53 (370 ):979–987. 10.1093/jexbot/53.370.979 11912240
Oleksiak MF , ChurchillGA, CrawfordDL. Variation in gene expression within and among natural populations. Nat Genet. 2002:32 (2 ):261–266. 10.1038/ng983 12219088
Oliveira HR , JacocksL, CzajkowskaBI, KennedySL, BrownTA. Multiregional origins of the domesticated tetraploid wheats. PLoS One. 2020:15 (1 ):e0227148. 10.1371/journal.pone.0227148 31968001
Ondov BD , TreangenTJ, MelstedP, MalloneeAB, BergmanNH, KorenS, PhillippyAM. Mash: fast genome and metagenome distance estimation using MinHash. Genome Biol. 2016:17 (1 ):132. 10.1186/s13059-016-0997-x 27323842
Oostra V , SaastamoinenM, ZwaanBJ, WheatCW. Strong phenotypic plasticity limits potential for evolutionary responses to climate change. Nat Commun. 2018:9 (1 ):1005. 10.1038/s41467-018-03384-9 29520061
Page A , GibsonJ, MeyerRS, ChapmanMA. Eggplant domestication: pervasive gene flow, feralization, and transcriptomic divergence. Mol Biol Evol. 2019:36 (7 ):1359–1372. 10.1093/molbev/msz062 31039581
Patro R , DuggalG, LoveMI, IrizarryRA, KingsfordC. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017:14 (4 ):417–419. 10.1038/nmeth.4197 28263959
Plett DC , HolthamLR, OkamotoM, GarnettTP. Nitrate uptake and its regulation in relation to improving nitrogen use efficiency in cereals. Semin Cell Dev Biol. 2018:74 :97–104. 10.1016/j.semcdb.2017.08.027 28843981
Porter SS , SachsJL. Agriculture and the disruption of plant–microbial symbiosis. Trends Ecol Evol. 2020:35 (5 ):426–439. 10.1016/j.tree.2020.01.006 32294424
Purcell S , NealeB, Todd-BrownK, ThomasL, FerreiraMAR, BenderD, MallerJ, SklarP, De BakkerPIW, DalyMJ, et al PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007:81 (3 ):559. 10.1086/519795 17701901
Rapaport F , KhaninR, LiangY, PirunM, KrekA, ZumboP, MasonCE, SocciND, BetelD. Comprehensive evaluation of differential gene expression analysis methods for RNA-seq data. Genome Biol. 2013:14 (9 ):R95. 10.1186/gb-2013-14-9-r95 24020486
R Core Team . R: a language and environment for statistical computing. Vienna, Austria: R Foundation Statistical Computing; 2022. https://www.R-project.org/
Roberge C , GuderleyH, BernatchezL. Genomewide identification of genes under directional selection: gene transcription QST scan in diverging Atlantic salmon subpopulations. Genetics. 2007:177 (2 ):1011–1022. 10.1534/genetics.107.073759 17720934
Robinson MD , McCarthyDJ, SmythGK. Edger: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010:26 (1 ):139–140. 10.1093/bioinformatics/btp616 19910308
Rockström J , SteffenW, NooneK, PerssonÅ, ChapinFS, LambinEF, LentonTM, SchefferM, FolkeC, SchellnhuberHJ. A safe operating space for humanity. Nature. 2009:461 (7263 ):472–475. 10.1038/461472a 19779433
Sauvage C , RauA, AichholzC, ChadoeufJ, SarahG, RuizM, SantoniS, CausseM, DavidJ, GléminS. Domestication rewired gene expression and nucleotide diversity patterns in tomato. Plant J. 2017:91 (4 ):631–645. 10.1111/tpj.13592 28488328
Savary S , WillocquetL, PethybridgeSJ, EskerP, McRobertsN, NelsonA. The global burden of pathogens and pests on major food crops. Nat Ecol Evol. 2019:3 (3 ):430–439. 10.1038/s41559-018-0793-y 30718852
Smyth GK . Limma: linear models for microarray data. In: GentlemanR, CareyVJ, HuberW, IrizarryRA, DudoitS, editors, Bioinformatics and computational biology solutions using R and bioconductor. New York (NY): Springer; 2005. p. 397–420.
Sultana N , IslamS, JuhaszA, YangR, SheM, AlhabbarZ, ZhangJ, MaW. Transcriptomic study for identification of major nitrogen stress responsive genes in Australian bread wheat cultivars. Front Genet. 2020:11 :583785. 10.3389/fgene.2020.583785 33193713
Swanson-Wagner R , BriskineR, SchaeferR, HuffordMB, Ross-IbarraJ, MyersCL, TiffinP, SpringerNM. Reshaping of the maize transcriptome by domestication. Proc Natl Acad Sci U S A. 2012:109 (29 ):11878–11883. 10.1073/pnas.1201961109 22753482
Tajima F . Evolutionary relationship of DNA sequences in finite populations. Genetics. 1983:105 (2 ):437–460. 10.1093/genetics/105.2.437 6628982
Taranto F , D’AgostinoN, RodriguezM, PavanS, MinerviniAP, PecchioniN, PapaR, De VitaP. Whole genome scan reveals molecular signatures of divergence and selection related to important traits in durum wheat germplasm. Front Genet. 2020:11 :217. 10.3389/fgene.2020.00217 32373150
Tian T , LiuY, YanH, YouQ, YiX, DuZ, XuW, SuZ. agriGO v2.0: a GO analysis toolkit for the agricultural community, 2017 update. Nucleic Acids Res. 2017:45 (W1 ):W122–W129. 10.1093/nar/gkx382 28472432
Valente J , GerinF, MiniA, RichardR, Le GouisJ, Prigent-CombaretC, Moënne-LoccozY. Symbiotic variations among wheat genotypes and detection of quantitative trait loci for molecular interaction with auxin-producing Azospirillum PGPR. Microorg. 2023:11 (6 ):1615. 10.3390/microorganisms11061615
Vigouroux Y , McMullenM, HittingerCT, HouchinsK, SchulzL, KresovichS, MatsuokaY, DoebleyJ. Identifying genes of agronomic importance in maize by screening microsatellites for evidence of selection during domestication. Proc Natl Acad Sci U S A. 2002:99 (15 ):9650–9655. 10.1073/pnas.112324299 12105270
Watterson GA . On the number of segregating sites in genetical models without recombination. Theor Popul Biol. 1975:7 (2 ):256–276. 10.1016/0040-5809(75)90020-9 1145509
Whitehead A , CrawfordDL. Neutral and adaptive variation in gene expression. Proc Natl Acad Sci U S A. 2006a:103 (14 ):5425–5430. 10.1073/pnas.0507648103 16567645
Whitehead A , CrawfordDL. Variation within and among species in gene expression: raw material for evolution. Mol Ecol. 2006b:15 (5 ):1197–1211. 10.1111/j.1365-294X.2006.02868.x 16626448
Whitlock MC , GuillaumeF. Testing for spatially divergent selection: comparing QST to FST. Genetics. 2009:183 (3 ):1055–1063. 10.1534/genetics.108.099812 19687138
Wood SW , CowieA. A review of greenhouse gas emission factors for fertiliser production 2004. For IEA Bioenergy Task 38. Available online: https://www.sciencetheearth.com/uploads/2/4/6/5/24658156/2004_wood_a_review_of_greenhouse_gas_emission_factors.pdf
Zhang ZH , JhaveriDJ, MarshallVM, BauerDC, EdsonJ, NarayananRK, RobinsonGJ, LundbergAE, BartlettPF, WrayNR, et al A comparative study of techniques for differential expression analysis on RNA-Seq data. PLoS One. 2014:9 (8 ):e103207. 10.1371/journal.pone.0103207 25119138
Zhu T , WangL, RodriguezJC, DealKR, AvniR, DistelfeldA, McGuirePE, DvorakJ, LuoM-C. Improved genome sequence of wild emmer wheat Zavitan with the aid of optical maps. G3 Genes, Genomes, Genet. 2019:9 (3 ):619–624. 10.1534/g3.118.200902
