Based on the stress results of atrazine in the germination stage by Wang et al. (2021), 20 different resistant millet varieties, widely planted in China, were selected for field experiments to further verify their responses to atrazine (Table 1). The experiment was conducted in 2019 at the 5th management zone of Yinlang Farm (124° 44′ 13.89″ E, 56° 29′ 53.60″ N) in Daqing City, Heilongjiang Province, China. This site had an average annual temperature of 4.6°C, an active accumulated temperature of >2800°C, annual average sunshine of 2782.5 h, and a frost-free period of 136 days.

The experiment was conducted using a completely randomized design with three replications. The plot size was 9.75 m², with 5-m rows, 0.65-m row spacing, and 3 lines per repeat. The sowing density was 600000 plants/hm². Following the experimental screening results of atrazine stress in the germination stage by Wang et al. (2021), a mixture of approximately 13.55 ml of 38% fully mixed atrazine suspension and 800 ml of water (approximately 11 μmol/L) was sprayed evenly on each plot after sowing for 1 day. Water was used as the control. Additionally, 37.5 kg/hm² urea, 225 kg/hm² diammonium phosphate, and 120 kg/hm² potassium sulfate were applied as base fertilizers, after which approximately 120 kg/hm² urea was applied at the heading stage. Jinan Tianbang Chemical Co., Ltd., Jinan, Shandong province, China, provided the 38% atrazine suspension.

At the 3-leaf seedling stage, fresh leaf tissues were used to determine the physiological and biochemical indexes, including the chlorophyll content through ethanol extraction and colorimetric methods described by Liu et al. (2020), malondialdehyde (MDA) content using the thiobarbituric acid method as described by Akbulut and Yigit (2010), and soluble protein content using the coomassie brilliant blue G-250-binding method as described by Bradford (1976), with some modifications. Briefly, the whole plant at the 3-leaf seedling stage was transferred into paper bags, heated at 105°C for 30 min, then dried at 80°C to constant weight, and weighed using a 1/10000 gram balance with three biological repetitions. Subsequently, the average was taken after each determination.

At the 7–8-leaf seedling stage, fully expanded young leaves of 3 plants with the same growth pattern were selected for each treatment. Subsequently, we measured net photosynthetic rates, transpiration rates, and stomatal conductance of their inverter clovers using a LI-COR 6400 portable photosynthetic instrument (LI-COR Inc., Lincoln, NE, United States) on sunny days (9:00 to 11:00 h) at 1500 μmol m^–2 s^–1 PPFD. Ambient water vapor pressure, CO₂ concentration, and leaf temperature were maintained at 1.2 ± 0.1 kPa, 400 μmol m^–2 s^–1, and 25 ± 1.0°C, respectively. The average values were taken for data analysis.

At maturity, five panicles were randomly selected from each treatment to determine the grain weight per panicle. The average value was taken for data analysis.

All measurements conducted were repeated at least in triplicate. Relative values between the treatment and control were used in the data analysis through repetitions to avoid differences in genotypes. Data analysis was conducted using SPSS-26.0 (SPSS, Chicago, IL, United States) to determine the weight value of each indicator through principal component analysis (PCA), as described by Dong et al. (2016). The standardized principal component (PC) scores were extracted from the correlation matrix, and those PCs with eigenvalues of ≥1.00 were selected. It should be emphasized that to make the trend of each index consistent, the relative value of MDA should be converted to its reciprocal. Ultimately, the subordinate function method was used to comprehensively evaluate the resistance of different millet genotypes to stress.

According to the fold changes of each DEG and DEM, the Pearson’s correlation coefficient was calculated. Before calculating the correlation, the data were preprocessed using the Z-value transformation method. Correlations corresponding to a coefficient of R² > 0.8 and p-value of <0.05 were selected. The relationship between the metabolome and transcriptome was illustrated using BMKCloud (see text footnote 1).

The leaves of the millet plant were used for real-time quantitative reverse transcription PCR (qRT-PCR). The method used followed the procedures described in a previous report (Sun et al., 2018). The primers were designed and synthesized according to the gene sequences in the millet sequence database (see text footnote 2) and the 25S RNA sequence was used as the endogenous control (Supplementary Table 1). At the end of the PCR for each sample, based on the fluorescence logarithmic graph, the appropriate threshold was chosen and the Ct values were obtained to analyze the transcript levels of each gene using the 2^–ΔΔCt method (Livak and Schmittgen, 2001). The Bio-Rad CFX96 Real-Time System (Bio-Rad, Hercules, CA, United States) and CFX Manager System software version 2.0 were used in following the manufacturer’s instructions in this study.

To analyze the sensitivity of different millet varieties to atrazine, eight indicators in different growth stages were detected and the data were subjected to PC factor analysis. Kaiser–Meyer–Olkin (KMO) and Bartlett’s ball tests were conducted on the relative value data of each trait using SPSS (Supplementary Table 2). The results showed that the KMO test value was 0.819 > 0.5 and the probability of significance in Bartlett’s ball test was p < 0.05, indicating that the test data conformed to spherical distribution (Table 2). Therefore, the experimental data of this study were suitable for factor analysis.

In Table 3, the eigenvalues of the first three PCs were >1, and the cumulative contribution rate reached 76.722%, reflecting most of the information of the 8 indexes. Therefore, the first three PCs were selected to evaluate the eight indexes of atrazine tolerance in millet. The confidence level of the evaluation was 76.722%. Based on PCA, the weight value of each indicator was calculated. Results are shown in Table 3. Ultimately, the comprehensive values of atrazine tolerance in 20 foxtail millet varieties were calculated (Table 4). This result showed that the T value of Gongai 2 (GA2) was the highest (0.91), whereas that of Longgu 31 (LG31) was the lowest (0.09), indicating that GA2 was not sensitive to atrazine, whereas LG31 was sensitive to atrazine. Therefore, GA2 and LG31 varieties were used further analysis.

For a deeper understanding of the molecular mechanisms occurring in the millet plant under atrazine stress, its leaves (Figure 1) were used for RNA-Seq, and the raw data of RNA-Seq were uploaded to the NCBI website. With three biological replicates, a total of 87.16 Gb of clean data with 94.31% of bases scoring Q30 were obtained. Of the total clean reads, 93.54–94.91% were unique matches with the Setaria italic L. reference genome and 35,747 unique genes were expressed in the millet leaves, of which 1,163 were novel genes and only 856 genes were successfully annotated (Supplementary Table 3).

The expression levels of the transcripts for each sample were calculated to analyze the different resistant and sensitive millet varieties. A total of 3,981 DEGs were identified in the samples of both varieties: 2,208 (501 upregulated and 1,707 downregulated) genes in the Gongai variety and 1,773 (761 upregulated and 1,012 downregulated) genes in the Longgu variety (FC ≥ 2 and FDR < 0.01; Table 5). The detailed information of the genes is presented in Supplementary Table 4. The results also showed that the number of downregulated genes (1,707 in GA2 and 1,012 in LG31) was significantly higher than the number of upregulated genes (501 in GA2 and 761 in LG31) under atrazine stress. The results demonstrated that atrazine could induce changes in transcript levels in millet leaves and that the expression of most genes may be inhibited to adapt to atrazine stress.

To classify the functions of DEGs, GO, and KEGG analyses were performed to clarify the biological pathways and functions in millet leaves affected by atrazine. The results showed that the significantly enriched GO categories were “metabolic process,” “cellular process,” and “single-organism process” in the resistant and sensitive millet varieties (Supplementary Figure 1). However, the KEGG analysis results implied that the resistant and sensitive millet varieties had different metabolic pathways (Supplementary Table 4). Based on the results presented in Figure 2A, when comparing GT with GCK, the DEGs enriched in “phenylalanine biosynthesis,” “flavonoid biosynthesis,” “flavone and flavonol biosynthesis,” “ribosome biogenesis,” and “ribosome” showed remarkable downregulation, whereas those enriched in “lysine biosynthesis,” “arginine and proline metabolism,” “fructose and mannose metabolism,” “pentose phosphate pathway,” and “terpenoid backbone biosynthesis” were significantly upregulated. In the “LCK vs. LT” comparison, except for the “ribosome,” “ribosome biogenesis in eukaryotes,” and “vitamin B6 metabolism” pathways, the other metabolic pathways with DEGs were significantly upregulated (Figure 2B).

For an overview of the metabolomic changes induced by atrazine treatment in millet, the widely targeted metabolomics was analyzed using LC–MS/MS. A total of 636 metabolites were obtained from all samples (Figure 3A and Supplementary Table 5) and divided into four groups according to PCA (Figure 3B). On the basis of FC of ≥1.5, the numbers of upregulated and downregulated metabolites in GCK and GT and LCK and LT were 82 and 95 and 110 and 120, respectively (Supplementary Table 5). Moreover, these DEMs were identified into 65 species and were mainly involved in phenylpropanoid biosynthesis (18), biosynthesis of secondary metabolites (59), metabolic pathways (82), glutathione metabolism (4), flavonoid biosynthesis (7), carbon metabolism (13), arginine and proline metabolism (8), glutathione metabolism (4), glyoxylate and dicarboxylate metabolism (7), purine metabolism (8), and others (Figure 3C and Supplementary Table 6). The DEMs related to proline metabolism, glutathione metabolism, and flavonoid metabolism were significantly accumulated in GCK and GT and LCK and LT. The results implied that atrazine can induce changes in antioxidant levels in millet leaves.

To analyze the relationship between DEGs and DEMs under atrazine stress in millet, a coexpression network analysis of DEGs and DEMs was performed (Pearson’s coefficient R² > 0.8 and p < 0.05; Supplementary Table 7). The results of the KEGG analysis between the DEGs and DEMs are presented in Supplementary Figure 2. The results implied that GA2 and LG31 millet varieties have different regulation mechanisms under atrazine stress. For example, the DEGs and DEMs in GCK vs. GT were mainly involved in flavonoid biosynthesis, glutathione metabolism, amino acid biosynthesis, and fructose/mannose metabolism (Supplementary Figure 2A). However, in LCK vs. LT, the galactose, starch and sucrose, and thiamine metabolisms were mainly enriched (Supplementary Figure 2B).

To further examine the relationship between DEGs and DEMs in millet under atrazine stress, a coexpression network analysis was conducted. Oxiglutatione, L-cysteine, and γ-glutamylcysteine were found to be related to glutathione metabolism (ko00480, Supplementary Figure 3 and Figures 4A,B) and more highly regulated in GCK vs. GT than in LCK vs. LT (Table 6). By contrast, except for glutathione synthesis, trans-5-o-(4-coumaroyl)shikimate and scopolin were highly enriched in LCK vs. LT and were correlated with phenylpropanoid biosynthesis (ko00940, Supplementary Figure 3 and Figures 4C,D). Moreover, the coexpression network of DEGs and DEMs in GCK vs. GT was mainly enriched in the biosynthesis of amino acids (ko01230, Supplementary Figure 3 and Figure 4E), such as L-cysteine, and N-acetyl-L-glutamic acid was more highly regulated in GCK vs. GT than in LCK vs. LT (Table 6).

To further analyze the effects of atrazine on genes and metabolites involved in glutathione biosynthesis in millet, the networks of DEGs and DEMs involved in glutathione metabolism pathway (ko00480) were analyzed (Figure 5A and Supplementary Table 8). The results showed that 12 and 4 genes were differentially expressed in GA2 and LG31 at transcriptional levels, respectively. Compared with the CK group, the expression of all genes, except Seita.9G06450 and Seita.4G241000, was downregulated in the treatment group. The functions of all these genes were related to glutathione S-transferase (EC: 2.5.1.18). In addition, other genes in this pathway, such as GCLC (encoding for glutamate cysteine ligase catalytic subunit, Seita.3G050000), GGT (encoding for γ-glutamyl transpeptidase, Seita.5G104100), and G6PDH (encoding for glucose-6-phosphate 1-dehydrogenase, Seita.9G424300), were downregulated, whereas GSS (encoding for glutathione synthase, Seita.3G348900) was upregulated in GA2 (Supplementary Table 9). DEMs such as oxiglutatione, L-cysteine, and γ-glutamylcysteine were related to glutathione metabolism and more highly regulated in GCK vs. GT than in LCK vs. LT. The results showed that the DEGs and DEMs related to glutathione metabolism in GA2 were involved in the response to atrazine stress.

To investigate the effect of atrazine stress on the expressions of genes and metabolites involved in amino acid biosynthesis, the interaction of DEGs and DEMs involved in amino acid biosynthesis (ko01230) was analyzed (Figure 5B and Supplementary Table 10). We found that the expression of the DEMs N-acetyl-L-glutamic acid and L-proline was more highly regulated in GCK vs. GT than in LCK vs. LT. However, almost all DEGs related to this pathway were downregulated in GA2 (Supplementary Table 9), including the genes that may be regulated by enzymes related to the two abovementioned metabolites, such as NAGS (encoding for N-acetyl-L-glutamate synthetase, Seita.7G188800) and P5CS (encoding for delta-1-pyrroline-5-carboxylate synthetase, Seita.3G220500), except for ArgC (encoding for N-acetyl-gamma-glutamyl-phosphate reductase, Seita.9G167800). The results showed that the DEGs and DEMs related to amino acid biosynthesis in GA2 were involved in the response to atrazine stress.

To examine the effects of atrazine stress on the expression of genes and metabolites involved in phenylpropanoid biosynthesis, the interaction of DEGs and DEMs involved in amino acid biosynthesis was analyzed (Figure 5C and Supplementary Table 10). DEGs such as ToGT1 (encoding for scopoletin glucosyltransferase, Seita.5G239500) and C3′H [encoding for 5-o-(4-coumaroyl)-D-quinate 3′-monooxygenase, Seita.3G194300] were upregulated in LG31 (Supplementary Table 9). Moreover, the DEMs trans-5-o-(4-coumaroyl)shikimate and scopolin were highly enriched in LCK vs. LT and were correlated to phenylpropanoid biosynthesis (ko00940). The results showed that the DEGs and DEMs related to phenylpropanoid biosynthesis in LG31 were jointly involved in the response to atrazine stress.

To validate the accuracy of the RNA-Seq data, the expression levels of approximately 34 DEGs (Supplementary Table 11) that were related to glutathione, phenylpropanoid, and amino acid biosyntheses, such as GSTU6, HST, GLN1-3, and GSTF1, were measured using qRT-PCR. The results were similar to those of RNA-Seq (log₂FC), which verified the reproducibility and credibility of the RNA-Seq data.

several studies indicated that atrazine stress can lead to excessive accumulation of ROS in plant cells, resulting in the peroxidation of membrane lipids and oxidative damage to the cell membrane, proteins, and DNA (Bai et al., 2015; Baxter et al., 2015; Erinle et al., 2018; Singh et al., 2018). However, to avoid damage caused by ROS, plants have evolved a complex detoxification system (Yuan et al., 2007). For example, plants can use ROS scavenging enzymes (Superoxide Dismutase, Peroxidases, and Catalase) and non-enzymatic antioxidants (glutathione, glutathione reductase, and glutathione-Px) to scavenge excess ROS (Sharif et al., 2018). If the excess ROS is not eliminated by detoxification system, which could be marked by the accumulation of MDA (Bai et al., 2015). Sathish et al. (2022) found MDA content in Zea mays L. was significantly increased by water deficit stress. In the present study, the MDA content in GA2 was lower than that in LG31 under atrazine stress, indicating that the degree of peroxidation damage caused by atrazine in GA2 was less to maintain relatively better membrane stability. Moreover, the activities of ROS scavenging enzymes were not well detected in this study. These results indicate that ROS scavenging in millet plants mainly occurs by increasing the levels of non-enzymatic antioxidants, which corroborates the findings of a previous study on fulvic acid used for alleviating drought stress in tea plants (Sun J. H. et al., 2020).

Glutathione is a thiol tripeptide with a low molecular weight that widely exists in all kinds of cells and functions as a non-enzymatic antioxidant that scavenges for H₂O₂, ¹O₂, OH^•, and O2•- and protects different biomolecules to maintain ROS homeostasis (Das and Roychoudhury, 2014). Glutathione S-transferases (GSTs) can convert various herbicides using glutathione into non-toxic and water-soluble forms, making the plant resistant and reducing the damage to crops caused by herbicides (Shimabukuro et al., 1970). In this study, glutathione levels were not significantly different in the atrazine-resistant GA2 and atrazine-sensitive LG31 variety. Oxidized glutathione (GSSG) levels were increased in GA2, but showed no change in LG31. These results imply that the dynamic equilibrium between glutathione and GSSG may fluctuate because glutathione participates in resistance to environmental stress in GA2. Moreover, a previous study reported that maize and grain sorghum, as crops naturally tolerant to atrazine, showed rapid detoxification resulting from the high constitutive activity of GSTs (Timmerman, 1989). Similarly, studies on Lolium rigidum in Australia and Alopecurus myosuroides in the United Kingdom have also proved the detoxification mechanisms of GST (Cummins et al., 2013; Yu and Powles, 2014). In the present study, nine genes with functions related to GST were detected in GA2, but the expressions of the GST genes in the atrazine treatment group were downregulated at the transcriptional level. The reason for this may be that gene transcription occurs before substance metabolism. Therefore, we deduced that the glutathione pathway plays an important role in the resistance of GA2 to atrazine stress.

Proline, an osmotic regulator, is the main amino acid for the stabilization and protection of cell membranes and proteins under certain stress conditions, such as drought, salinity, and low and high temperatures (Rivero et al., 2014; Hasanuzzaman et al., 2019). The biosynthesis of proline mainly depends on glutamate, and the rate-limiting step in this pathway is catalyzed by pyrroline-5-carboxylate synthetase (P5CS; Rejeb et al., 2014). Previous studies reported that P5CS1 expression can be induced by drought, salt, abscisic acid (ABA), light, nitric oxide, and pathogens to improve the concentration of proline in Arabidopsis (Savouré et al., 1995; Strizhov et al., 1997; Hayashi et al., 2000; Fabro et al., 2004; Zhao et al., 2009). In general, tolerant genotypes of crop plants accumulated more osmoprotectants with high concentrations than sensitive genotypes under stress. For example, a near-isogenic line of spring barley genotypes (NIL143) that contains the wild allele pyrroline-5-carboxylate synthase 1-P5cs1 was observed higher proline concentrations in leaf and root and less severe symptoms of drought compared with any other genotypes under reduced water availability conditions (Frimpong et al., 2021).

In our study, we found the DEGs related to the proline biosynthesis pathway were all downregulated, including P5CS (Seita.3G220500). While proline was significantly accumulated in the GA2 group under atrazine stress. Many researchers also suggested that proline can not only improve osmotic tolerance in plants cell but also reduce ROS generation and scavenge excessive free ROS (Smirnoff and Cumbes, 1989; Alia et al., 1991; Shalata and Neumann, 2001). Recent studies also showed a positive correlation between proline and ROS levels under abiotic stress or exogenous treatment with H₂O₂ (Yang et al., 2009; Lopez-Delacalle et al., 2020). Previous studies showed that many metabolic mechanisms are involved in the response of plants to abiotic stress, including the production of signal molecules, activation and inactivation of gene expression, and synthesis of functional proteins (Zushi et al., 2014; Zandalinas et al., 2020). Among these, gene expression is rapid and occurs before protein synthesis. Thus, we only detected significantly increased levels of proline in the atrazine-tolerant GA2 variety, suggesting that proline may help in scavenging excessive ROS produced by atrazine stress and protecting cells from damage.

To protect plant themselves from various environmental stress factors, plant cells have developed an inducible defense mechanism that involves the production of toxic chemicals or antimicrobial compounds, such as secondary metabolites (Kliebenstein and Osbourn, 2012; Tiku, 2018). Scopoletin and its glucoside scopolin are coumarins, which are one of the secondary metabolite classes (Kai et al., 2006). The biosynthesis of scopoletin mostly occurs via the shikimic acid pathway (Bayoumi et al., 2008). Previous studies reported that treatment with CuCl₂, shortwave UV, and Triton X-100 induces scopoletin accumulation in sunflower leaves (Gutierrez et al., 1995). Scopoletin is also produced in leaves in response to pathogen attack (Barón et al., 2016). Moreover, its role in abiotic stress may be as an antioxidant for scavenging ROS to reduce oxidative stress in plant cells (Chong et al., 2002; Döll et al., 2018). In this study, DEMs and DEGs in the synthesis pathways of scopolin and metabolites [scopolin and trans-5-o-(4-coumaroyl)shikimate] as well as genes (ToGT1 and C3′H) were significantly upregulated in LG31 and downregulated in GA2. These results showed that scopolin may be accumulated in LG31 as an antioxidant to scavenge excessive ROS produced by atrazine stress.