Abstract

Background

   Alfalfa (Medicago sativa L.) is a vital forage crop with substantial
   economic and ecological significance in agriculture and animal
   husbandry. However, atrazine, a widely used herbicide, negatively
   impacts the growth and yield of alfalfa due to its residual presence in
   the environment. Transcriptomic analysis was performed to investigate
   the differences in tolerance and uncover the potential molecular
   regulatory mechanisms between the tolerant variety JN5010 and the
   sensitive variety WL363 when subjected to atrazine stress, using
   RNA-seq on pooled samples.

Results

   Based on the analysis of gene expression profiles, significant
   differences were observed between the tolerant variety JN5010 and the
   sensitive variety WL363 under atrazine stress: 2,297 upregulated and
   3,167 downregulated in the shoot parts, and 3,232 upregulated and 4,907
   downregulated in the roots of JN5010. In WL363, 2,937 genes were
   upregulated and 4,237 genes were downregulated in the shoot parts,
   while 5,316 genes were upregulated and 7,977 genes were downregulated
   in the roots. The DEGs in the shoot parts were mainly involved in
   biological regulation, metabolic processes, and cellular processes,
   including proline metabolic processes and S-adenosylmethionine cycle.
   The DEGs in the roots were predominantly associated with nitric oxide
   synthesis and metabolism, as well as processes related to cell wall
   biosynthesis and degradation. In the shoot parts of JN5010, six DEGs
   were mapped onto the proline metabolic pathway, including four
   upregulated genes involved in proline biosynthesis and two
   downregulated genes involved in proline catabolism. In the roots of
   WL363, eleven DEGs were mapped onto the phenylpropanoid biosynthesis
   pathway, including seven upregulated genes involved in flavonoid
   biosynthesis and four downregulated genes associated with lignin
   biosynthesis. These findings highlight the distinct genetic responses
   of the two alfalfa varieties to atrazine stress, with JN5010 exhibiting
   more consistent gene expression patterns compared to the sensitive
   variety WL363.

Conclusions

   The tolerant variety JN5010 shows improved tolerance to atrazine stress
   by maintaining stable gene expression and precise regulation in various
   pathways, such as antioxidant processes, signaling, photosynthesis, and
   toxin removal. This differential gene expression helps JN5010 maintain
   stability in its functions under stress, demonstrating better
   adaptability. These findings enhance our understanding of how alfalfa
   tolerates atrazine stress and provide important insights for developing
   atrazine-tolerant varieties.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12870-025-07129-x.

   Keywords: Alfalfa, Atrazine, Abiotic stress, Oxidative damage,
   Comparative transcriptomics

Introduction

   Alfalfa (Medicago sativa L.), often referred to as the “king of
   forages,” holds immense economic and ecological value in agriculture
   and animal husbandry. Renowned for its ability to provide high-quality
   feed for livestock, it also enhances soil fertility through nitrogen
   fixation [[44]1]. However, research has demonstrated that alfalfa
   monoculture systems cultivated for over six years display a notable
   decline in yield, accompanied by a degradation of soil microbial
   community structure and an increase in pathogens [[45]2]. Furthermore,
   alfalfa releases allelopathic compounds to inhibit subsequent alfalfa
   growth, resulting in continuous cropping challenges [[46]3, [47]4].
   Consequently, crop rotation with corn (Zea mays L.) has been widely
   adopted. Studies revealed that alfalfa-corn rotation improves soil pH,
   reduces salinity, and enhances organic matter content, thereby boosting
   soil fertility [[48]5]. The rotation system also significantly
   increased soil nitrogen utilization through heightening soil nitrogen
   content, enhancing phosphorus solubility, and stimulating of microbial
   activity. These benefits not only reduce the reliance on chemical
   fertilizers but also improve the yield and quality of corn. The
   effectiveness of this rotational system becomes more pronounced with
   the increasing duration of alfalfa cultivation, markedly enhancing land
   use efficiency and economic returns [[49]6]. Nonetheless, the
   widespread use of atrazine, a primary herbicide for controlling corn
   weeds, poses a challenge. Atrazine residues can hinder the growth and
   development of subsequent crops, ultimately leading to substantial
   yield losses.

   Atrazine is a commonly utilized triazine herbicide, primarily applied
   for weed control in crops such as corn and sugarcane. Despite its
   notable effectiveness in agricultural production, the environmental and
   health risks it poses are significant and warrant attention. Research
   has shown that atrazine accumulation in drinking water may pose
   reproductive and developmental risks to children [[50]7]. Moreover,
   with a half-life of up to 57 weeks in soil, atrazine leaves persistent
   residues that can harm subsequent crops and contaminate groundwater and
   surface water, endangering aquatic ecosystems [[51]8]. While the
   European Union banned atrazine in 2003, it continues to be widely used
   in agricultural powerhouse nations such as China and Brazil due to its
   high efficacy and versatile application methods [[52]9]. Studies have
   demonstrated that atrazine residues exhibit toxic effects on various
   crops, including cucumber, soybean, and peanut [[53]10]. Atrazine
   disrupts photosystem II (PS II), inhibits the electron transport chain,
   and reduces photosynthetic efficiency, which leads to the accumulation
   of reactive oxygen species (ROS) [[54]11]. This accumulation of ROS can
   cause oxidative stress and cellular damage. In recent years,
   advancements in transcriptomics and metabolomics have provided novel
   insights into the toxic mechanisms of atrazine on plants. Research has
   revealed that, under atrazine stress, genes involved in antioxidant
   defense (such as peroxidase (POD) and superoxide dismutase (SOD))are
   significantly upregulated in rice as a response to ROS accumulation
   [[55]12]. Metabolomic analyses have identified substantial alterations
   in amino acid metabolism and secondary metabolites, such as flavonoids
   and phenolic compounds, which appear to contribute to the plant’s
   detoxification and adaptive mechanisms [[56]13]. Further studies
   indicate that under atrazine stress, genes encoding cysteine synthase
   and spermidine synthase in alfalfa seedlings were significantly
   upregulated. This upregulation likely enhances the antioxidant capacity
   of alfalfa, aiding in its response to atrazine-induced oxidative stress
   [[57]14]. However, there are considerable variances in the tolerance of
   different alfalfa varieties to atrazine, and the specific molecular
   mechanisms underlying these differences remain poorly understood. Thus,
   elucidating the molecular response mechanisms in various alfalfa
   varieties under atrazine stress is crucial for a deeper understanding
   of atrazine detoxification and degradation strategies.

   In this study, we selected two alfalfa varieties with contrasting
   tolerance characteristics: JN5010, a tolerant variety, and WL363, a
   sensitive variety, both identified through previous field trials.
   Leveraging high-throughput RNA-Seq technology, we conducted a
   systematic comparison and analysis of their growth parameters and
   transcriptomic responses under atrazine stress. The objective is to
   integrate physiological indicators (including SOD, MDA, chlorophyll
   content, and soluble sugar levels) with transcriptomic sequencing and
   qRT-PCR validation to identify differentially expressed genes (DEGs)
   between the tolerant and sensitive varieties exposed to atrazine. The
   molecular response mechanisms of alfalfa tolerance to atrazine were
   elucidated through the integrated analysis of transcriptomic and
   phenotypic data. Furthermore, we sought to elucidate dynamic changes in
   key metabolic and signal transduction pathways, as well as to clarify
   the molecular response mechanisms distinguishing these tolerant
   varieties. This research provides valuable theoretical insights into
   the molecular mechanisms underlying alfalfa tolerance to atrazine,
   establishing a crucial foundation for molecular breeding practices
   aimed at developing tolerant varieties.

Materials and methods

Experimental materials and pretreatment

   Atrazine (suspension concentrate with a purity of 38%) was obtained
   from Qiaochang Modern Agriculture Company. The atrazine-tolerant
   alfalfa variety JN5010 was purchased from Zhengzhou Huafeng Grass
   Industry Technology Co., Ltd., and the sensitive variety WL363 was
   purchased from Beijing Zhengdao Seed Industry Co., Ltd. Seeds of both
   varieties were disinfected with 5% sodium hypochlorite and thoroughly
   rinsed with distilled water after 10 min of treatment. Petri dishes
   were lined with 10 cm × 10 cm filter papers, to which an adequate
   amount of sterile water was added to facilitate seed germination. The
   dishes were then placed in a light-incubated chamber. After five days,
   seedlings exhibiting uniform growth were selected, with six pots
   allocated per variety and 20 seedlings planted in each pot. The
   atrazine culture system used a modified Hoagland nutrient solution at
   half strength, which was replaced every three days. After three weeks,
   three pots from each of the six pots per variety were subjected to
   stress treatment with 2.0 mg/L atrazine (38%) for six days, during
   which the nutrient solution was not replaced. From each pot, 10
   uniformly growing alfalfa plants were carefully removed from the
   atrazine solution, quickly rinsed with sterilized DEPC water, and
   blotted dry with filter paper to remove excess moisture. Measurements
   were taken for plant height and the natural length of the primary
   roots. Subsequently, the shoot and root biomass of each plant were
   recorded. The stems and roots of the plants were then cut into segments
   roughly 1 cm in length, rapidly frozen using liquid nitrogen, and
   stored at -80 °C for subsequent total RNA extraction.

Measurement of physiological indicators

   Chlorophyll content in leaves and soluble sugar levels of two alfalfa
   varieties were determined using Solarbio kits (BC0995 and BC0035;
   Solarbio, Beijing, China). Superoxide dismutase (SOD) activity and
   malondialdehyde (MDA) content were measured with Abbkine kits (KTB1030
   and KTB1050; Abbkine, Wuhan, China). Chlorophyll was quantified only in
   the shoot parts (stems and leaves), while all other parameters were
   assayed separately in both roots and aerial tissues (stems and leaves).
   To ensure accuracy and reliability, each physiological parameter was
   assessed in triplicate.

RNA extraction and quality control

   Total RNA was separately extracted from the shoot parts (stems and
   leaves) and roots of 20 alfalfa seedlings using TRIzol^® Reagent
   (Invitrogen, California, USA) following the manufacturer’s instructions
   (Pooled tissues), to reduce the impact of individual variation. Genomic
   DNA contamination was eliminated through treatment with DNase I
   (TaKaRa). The quality and concentration of the extracted RNA were
   evaluated using a 2100 Bioanalyzer (Agilent Technologies, California,
   USA) and a NanoDrop ND-2000 spectrophotometer (NanoDrop Technologies,
   Delaware, USA). Only RNA samples meeting stringent quality criteria
   (OD260/280 = 1.8–2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5, 28 S:18 S ≥ 1.0, and
   > 1 µg yield) were selected for library construction and sequencing.

Library Preparation and illumina sequencing

   RNA-seq libraries were prepared using the TruSeq™ RNA Sample
   Preparation Kit (Illumina, San Diego, CA) with 1 µg of total RNA.
   Messenger RNA was isolated via the poly(A) selection method with
   oligo(dT) beads and subsequently fragmented. Double-stranded cDNA was
   synthesized using the SuperScript Double-Stranded cDNA Synthesis Kit
   (Invitrogen, California, USA) with random hexamer primers (Illumina).
   The resulting cDNA underwent end-repair, phosphorylation, and ‘A’ base
   addition, followed by size selection for 300 bp fragments using a 2%
   Low Range Ultra Agarose gel. The libraries were then amplified through
   15 PCR cycles using Phusion DNA polymerase (NEB). Once quantified, the
   libraries were sequenced on the Illumina HiSeq X Ten/NovaSeq 6000
   platform, producing paired-end reads with a length of 150 bp. See Fig.
   [58]S1 for the sequencing data statistics table. The raw sequencing
   data were uploaded to the CNCB SRA database (Shoot: CRA024145; Root:
   CRA024142).

Read mapping and transcript assembly

   The raw paired-end reads were trimmed and quality-controlled using
   SeqPrep ([59]https://github.com/jstjohn/SeqPrep) and Sickle
   ([60]https://github.com/najoshi/sickle) with default parameters. The
   clean reads were then aligned to the reference genome using the HISAT2
   software ([61]http://ccb.jhu.edu/software/hisat2/index.shtml) in an
   orientation-based mode. The mapped reads of each sample were assembled
   using StringTie
   ([62]https://ccb.jhu.edu/software/stringtie/index.shtml?t=example) in a
   reference-based approach.

Differential expression analysis and functional enrichment

   The expression levels of individual transcripts were quantified using
   the fragments per kilobase of transcript per million mapped reads
   (FPKM), calculated via RSEM
   ([63]http://deweylab.biostat.wisc.edu/rsem/). Differential expression
   analysis was conducted using DESeq2, DEGseq, or EdgeR, with a
   q-value ≤ 0.05 and|log[2]FC| >1 considered thresholds for identifying
   significantly DEGs. Functional enrichment analyses, including Gene
   Ontology (GO) classification and Kyoto Encyclopedia of Genes and
   Genomes (KEGG) pathway evaluation, were performed using Goatools
   ([64]https://github.com/tanghaibao/Goatools) and KOBAS
   ([65]http://kobas.cbi.pku.edu.cn/home.do), respectively, applying a
   Bonferroni-corrected p-value threshold of ≤ 0.05.

Weighted gene Co-expression network analysis and Protein-Protein interaction
network analysis

   Weighted gene co-expression network analysis (WGCNA) was applied to
   build networks of co-expressed genes and to delineate modules whose
   members share highly similar expression profiles. Modules exhibiting
   significant correlations with phenotypic traits were then selected for
   in-depth investigation of how core module genes relate to the observed
   phenotypes. For full implementation details see the [WGCNA package]
   ([66]https://cran.r-project.org/web/packages/WGCNA/index.html).
   Separately, proteins encoded by the set of DEGs were queried in the
   STRING database (v11.5; [67]https://string-db.org) to assemble a
   protein–protein interaction map, thereby uncovering putative functional
   links among these genes.

qRT-PCR analysis

   To validate the accuracy of our RNA-seq data, eight DEGs were selected
   for Real-Time Quantitative PCR (qRT-PCR) analysis. These genes were
   chosen from our dataset and have also been previously reported to be
   involved in stress responses. For the analysis, total RNA was first
   treated with gDNA Eraser to remove genomic DNA contamination and
   subsequently reverse-transcribed into cDNA using the PrimeScript™ RT
   Reagent Kit (TaKaRa, Cat# RR047A, Tokyo, Japan). The qRT-PCR was then
   performed using the resulting cDNA as a template on a CFX96 Real-Time
   PCR Detection System (Bio-Rad, USA). Each gene was assayed in three
   independent biological replicates, with MsACTIN2 used as the endogenous
   control for normalization. Gene-specific primers (Table. S2) were
   designed using Primer-BLAST (NCBI) and Oligo 7 software, with the
   following parameters: 18–24 nt length, 40–60% GC content, 58–62 °C
   melting temperature (Tm), and producing an 80–200 bp amplicon. The
   specificity of each primer pair was confirmed by melting curve analysis
   and subsequent Sanger sequencing of the PCR products.

Data analysis

   This study employed a triplicate experimental design, with each
   reported result representing the average of three replicate treatments,
   each consisting of 20 seedlings. Statistical analyses were conducted to
   identify DEGs across the libraries. Gene expression levels, or
   transcript abundance, were normalized to fragments per kilobase of
   transcript per million mapped reads (FPKM). DEGs between the libraries
   were identified using a Student’s t-test, with p-values adjusted for
   multiple testing to control the false discovery rate (FDR). The
   correlation among detected counts in parallel libraries was evaluated
   using the Pearson correlation coefficient. Biochemical and molecular
   response assessments involved determining significant differences
   between treatments through standard deviation and analysis of variance
   (ANOVA). ANOVA was also used for statistical comparisons across
   treatment groups, and when results showed significance at P < 0.05, a
   least significant difference (LSD) test was performed. All statistical
   analyses were carried out using SPSS version 27.0.

Result

Effects of atrazine stress on alfalfa growth and physiology

   In our previous field experiments, we identified the atrazine-tolerant
   variety ‘JN5010’ and atrazine-sensitive variety ‘WL363’ from 20 alfalfa
   varieties (Fig. [68]S1). In this study, the atrazine-tolerant variety
   ‘JN5010’ and atrazine-sensitive variety ‘WL363’ were exposed to
   atrazine stress, with subsequent quantification of critical growth
   parameters. The experimental results demonstrate that atrazine stress
   significantly impacts the growth and development of alfalfa seedlings
   (Fig. [69]1A). The atrazine-tolerant variety JN5010 and the sensitive
   variety WL363 exhibited markedly different responses to atrazine
   exposure. No significant change in plant height was observed in either
   variety (Fig. [70]1B, P > 0.05). However, root length was significantly
   inhibited in both varieties under atrazine stress (Fig. [71]1C). The
   reduction was more pronounced in WL363 (19.5%) compared to JN5010
   (12.1%) (P < 0.05). A clear distinction in biomass accumulation was
   observed between the two varieties (Fig. [72]1D, E). For the tolerant
   variety JN5010, atrazine exposure led to a significant reduction in
   biomass, with a 42.5% decrease in shoot weight and a 50.0% decrease in
   root weight (P < 0.05). In contrast, the biomass of the sensitive
   variety WL363 remained statistically unchanged, showing no significant
   alterations in either shoot or root weight (P > 0.05).

Fig. 1.

   [73]Fig. 1
   [74]Open in a new tab

   Effects of atrazine on the growth and physiological characteristics of
   two alfalfa varieties. Shown are the overall phenotype (A), plant
   height (B), root length (C), shoot weight (D), root weight (E),
   malondialdehyde (MDA) content (F), superoxide dismutase (SOD) activity
   (G), soluble sugar content (H), and chlorophyll content (I) of JN5010
   (tolerant) and WL363 (sensitive) seedlings. Data are presented as the
   mean ± standard deviation from three biological replicates. Different
   letters denote significant differences at P < 0.05

   Regarding antioxidant defense and oxidative stress, the two varieties
   showed different response patterns (Fig. [75]1F, G). Under atrazine
   stress, malondialdehyde (MDA) content increased significantly
   (P < 0.05) in both the shoots and roots of both varieties. However, the
   sensitive variety WL363 suffered more severe oxidative damage. After
   stress treatment, its shoot MDA content reached 84.76 nmol/g, a level
   significantly higher (P < 0.05) than the 57.79 nmol/g observed in the
   tolerant variety JN5010. Relative to their respective control levels,
   the MDA content increased by 20.0% in WL363 and 24.8% in JN5010. For
   superoxide dismutase (SOD) activity, a key distinction was observed in
   the roots: JN5010 root SOD activity was significantly induced by 40.5%,
   from 142.85 to 200.72 U/g (P < 0.05), whereas the root SOD activity in
   WL363 showed no significant change (P > 0.05). No significant
   alterations in shoot SOD activity were observed for either variety
   (P > 0.05). In terms of soluble sugar content (Fig. [76]1H), both
   varieties exhibited a slight but significant decrease in their shoots
   (P < 0.05). In the roots, a significant increase of 47.4% was recorded
   for JN5010 (P < 0.05), while no significant change was observed for
   WL363 (P > 0.05). Photosynthetic pigment content also differed
   substantially (Fig. [77]1I). In the tolerant variety JN5010,
   chlorophyll a content remained stable (P > 0.05), while chlorophyll b
   and total chlorophyll content significantly decreased by 12.1% and
   12.6%, respectively (P < 0.05). In contrast, the photosynthetic system
   of the sensitive variety WL363 was severely compromised. Its
   chlorophyll a content decreased sharply by 36.4% (from 0.77 to
   0.49 mg/g), and its chlorophyll b content decreased by 45.5% (from 0.55
   to 0.30 mg/g), leading to a 37.6% reduction in total chlorophyll
   content (P < 0.05).

Gene segregation in tolerant alfalfa under atrazine stress

   A principal component analysis (PCA) was conducted on the four gene
   libraries derived from the shoot tissues of two varieties, JN5010 and
   WL363, under both atrazine-treated and control conditions. The analysis
   revealed a distinct separation between the atrazine treatment group and
   the control group (Fig. [78]2A). Notably, in the PCA, the gene
   libraries subjected to atrazine stress-JN5010_ATZ and
   WL363_ATZ-clustered more closely together compared to the control
   groups, JN5010_CK and WL363_CK. This pattern suggests that under
   atrazine stress, a significant number of shared genes in the shoot
   tissues of both varieties exhibited differential expression, resulting
   in greater similarity in gene expression profiles under stress
   conditions. Venn diagram analysis further revealed a total of 71,338
   genes identified across the four gene libraries (Fig. [79]2C). Among
   these, 16,430 genes were commonly expressed in the treatment and
   control groups of JN5010, while16,759 genes were uniquely expressed in
   both treatment and control groups of WL363. In contrast, 13,607
   uniquely expressed genes were identified in JN5010_CK and WL363_CK,
   compared to 11,688 uniquely expressed genes in JN5010_ATZ and
   WL363_ATZ. These findings indicate that the number of specifically
   expressed genes between the two treatment groups exceeds the number of
   uniquely expressed genes between the two varieties. Additionally, under
   atrazine treatment, the number of uniquely expressed genes between the
   two varieties was lower than that observed under control conditions.

Fig. 2.

   [80]Fig. 2
   [81]Open in a new tab

   Comprehensive analysis of gene expression patterns in alfalfa under
   treatment conditions. (A) Principal component analysis (PCA) of shoot
   and root transcriptomes showing distinct clustering by treatment. (B)
   Root-specific PCA demonstrating treatment group separation. (C) Venn
   diagram of shoot genes with shared and unique responses. (D) Venn
   diagram of root genes highlighting treatment-specific signatures

   Principal component analysis revealed a clear separation between the
   atrazine treatment group and the control group in the root samples
   (Fig. [82]2B). Unlike the shoot tissues, the separation between the
   treatment groups (JN5010_ATZ and WL363_ATZ) on the principal coordinate
   plot was more pronounced than that of their respective control groups
   (JN5010_CK and WL363_CK). Additionally, the distance between the
   treatment and control groups of the WL363 variety was greater than that
   observed in the JN5010 variety. The Venn diagram analysis identified a
   total of 70,989 genes commonly detected across the four root
   transcriptome libraries (Fig. [83]2D). Among these, 15,034 genes were
   specifically expressed in the treatment and control groups of the
   JN5010 variety, while 18,163 specific genes were detected in WL363.
   Under control conditions, 11,657 genes were specifically expressed
   between the two varieties, whereas this number increased to 12,322
   under atrazine treatment. These findings demonstrate that the number of
   specifically expressed genes between the treatment groups of the two
   varieties exceeded that observed between their respective control
   groups. Moreover, under atrazine treatment, the number of specific
   genes expressed between the two varieties was higher than under control
   conditions. This evidence further suggests that under atrazine-induced
   stress, a significant number of genes in the roots of both alfalfa
   varieties exhibited specific expression, with WL363 exhibiting a
   stronger induction of gene expression compared to JN5010.

   Under atrazine stress, notable shifts in gene expression were observed
   in the shoot and root tissues of two alfalfa varieties, JN5010 and
   WL363. In the shoot tissues of JN5010, 2,297 genes were upregulated,
   while 3,167 genes were downregulated. Meanwhile, in WL363, 2,937 genes
   were upregulated, with 4,237 genes downregulated (Fig. [84]3A, B). In
   the root tissues, JN5010 showcased 3,232 upregulated genes and 4,907
   downregulated genes, whereas WL363 exhibited 5,316 upregulated genes
   alongside 7,977 downregulated genes (Fig. [85]3C, D). These findings
   suggest that gene expression in both shoot and root tissues of the
   alfalfa varieties underwent significant alterations under atrazine
   stress. Plants respond to abiotic stress by modulating gene expression
   to mitigate adverse conditions. PCA and Venn diagram analyses revealed
   that atrazine stress influenced the transcriptomic profiles of both
   alfalfa varieties in their roots and shoot tissues. Intriguingly, gene
   expression in the shoot parts became more uniform under atrazine
   stress, while the root tissues exhibited a noticeable degree of
   differentiation. This indicates that, despite changes in gene
   expression patterns in both shoot and root tissues, the mechanisms of
   response varied between these regions. Furthermore, a substantial
   number of common genes were affected in the shoot tissues, in contrast
   to the roots, which demonstrated distinct gene expression changes
   between JN5010 and WL363 under atrazine stress.

Fig. 3.

   [86]Fig. 3
   [87]Open in a new tab

   Heatmaps and bar charts of DEGs in shoot and root tissues between two
   alfalfa varieties under treatment and control conditions. (A) Heatmap
   of shoot tissues; (B) Bar chart of shoot tissues showing the number of
   up- and down-regulated DEGs; (C) Heatmap of root tissues; (D) Bar chart
   of root tissues showing the number of up- and down-regulated DEGs

GO annotation and enrichment analysis of DEGs under atrazine stress

   GO annotation analysis was conducted on the DEGs identified from the
   comparisons between shoot and root libraries, with the results
   presented in Fig. [88]4A. Under atrazine stress, the DEGs in the shoot
   and root tissues of alfalfa were predominantly categorized into three
   main groups: Biological Process (BP), Cellular Component (CC),
   Molecular Function (MF). Significant biological processes included
   biological regulation, metabolic processes, and cellular processes. The
   affected cellular components encompassed protein-containing complexes,
   organelle parts, organelles, membrane parts, and cell parts. Key
   molecular functions primarily involved catalytic activity and binding
   functions.

Fig. 4.

   [89]Fig. 4
   [90]Open in a new tab

   GO annotation and functional enrichment of DEGs between two alfalfa
   cultivars (JN5010 and WL363). (A) GO classification of DEGs into
   Biological Process, Cellular Component, and Molecular Function
   categories, showing results for root (R) and shoot (S) tissues. (B)
   Functional enrichment analysis of DEGs in shoot tissues. (C) Functional
   enrichment analysis of DEGs in root tissues

   A functional enrichment analysis was conducted on the DEGs in the shoot
   parts and roots of two alfalfa varieties under atrazine stress. The top
   10 enriched functions were compared, as illustrated in Fig. [91]4B, C.
   In the shoot parts, the proline metabolic process was significantly
   enriched in both JN5010 and WL363. Notably, JN5010 exhibited unique
   enrichments, including the S-adenosylmethionine cycle, meristem
   development, establishment of plastid localization, chloroplast
   localization, chloroplast migration, glycine biosynthesis, plastid
   localization, glycine biosynthesis via the glyoxylic acid transaminase
   pathway, and iron ion transport. Conversely, WL363 demonstrated unique
   enrichments among its top 10 significantly enriched functions, such as
   pre-replicative complex assembly involved in nuclear DNA replication,
   pre-replicative complex assembly involved in cell cycle DNA
   replication, pre-replicative complex components, initiation of DNA
   replication during mitosis, branched-chain amino acid catabolism,
   initiation of nuclear DNA replication, initiation of cell cycle DNA
   replication, and the proline catabolic process converting proline to
   glutamate.

   In the roots, the synthesis and metabolism of nitric oxide were
   significantly enriched in both JN5010 and WL363 within the DEGs. Among
   the top 10 enriched DEGs, JN5010 exhibited root-specific functions such
   as cysteine biosynthesis from cystathionine, the conversion of proline
   to glutamate, proline catabolism, seed development, β-glucan metabolic
   processes in the cell wall, and the biosynthesis and metabolism of
   cellulose in the plant cell wall. In contrast, WL363 displayed distinct
   functions among the top 10 enriched DEGs, including meiosis II, meiotic
   cell cycle processes in male meiosis II, saponin biosynthesis,
   glucoside biosynthesis, as well as pre-replicative complex formation
   and assembly essential for cell cycle DNA replication.

KEGG pathway enrichment analysis of DEGs under atrazine stress

   Employing KEGG database for functional annotation analysis facilitates
   the classification of DEGs based on their roles in metabolic pathways
   and biological functions. As shown in Fig. [92]S2, the DEGs from the
   shoot parts and roots of two alfalfa varieties subjected to atrazine
   stress are primarily grouped into the following categories: Organismal
   Systems, Genetic Information Processing, Cellular Processes,
   Environmental Information Processing, and Metabolism. Within the
   Genetic Information Processing category, the DEGs are involved in
   pathways linked to Environmental Adaptation, Transcription, Replication
   and Repair, Folding, Sorting and Degradation, and Translation. In the
   Cellular Processes category, the DEGs are associated with Transport and
   Catabolism. For the Environmental Information Processing category, they
   are connected to pathways such as Membrane Transport and Signal
   Transduction. The Metabolism category spans several pathways, including
   Carbohydrate Metabolism, Amino Acid Metabolism, Biosynthesis of
   Secondary Metabolites, Energy Metabolism, Lipid Metabolism, Terpenoid
   Metabolism, Polyketide Metabolism, Other Amino Acid Metabolism,
   Cofactor and Vitamin Metabolism, Glycine Biosynthesis and Metabolism,
   as well as Nucleotide Metabolism.

   The functional enrichment analysis of DEGs in the shoot and root parts
   of two alfalfa varieties revealed both shared features and distinct
   differences among the top 20 significantly enriched pathways
   (Fig. [93]5). Pathways significantly enriched in both varieties include
   β-alanine metabolism, amino sugar and nucleotide sugar metabolism,
   cysteine and methionine metabolism, galactose metabolism,
   alanine-aspartate-glutamate metabolism, propionate metabolism, starch
   and sucrose metabolism, glycine-serine-threonine metabolism,
   photosynthesis-antenna proteins, glycolysis/gluconeogenesis,
   selenocompound metabolism, valine-leucine-isoleucine degradation, and
   glyoxylic acid and dicarboxylic acid metabolism. Distinctly, the JN5010
   variety demonstrates unique enrichment in pathways such as pyruvate
   metabolism, carbon fixation in photosynthetic organisms, peroxisome
   metabolism, tyrosine metabolism, phenylpropanoid biosynthesis, fructose
   and mannose metabolism, and fatty acid degradation. Conversely, the
   WL363 variety shows unique pathway enrichment related to ribosome
   function, ascorbate and aldarate metabolism, phenylpropanoid
   biosynthesis, pyrimidine metabolism, cyanogenic amino acid metabolism,
   phagosome formation, and histidine metabolism.

Fig. 5.

   [94]Fig. 5
   [95]Open in a new tab

   Presents the KEGG enrichment analysis of DEGs in alfalfa. (A) shows the
   distribution of enriched metabolic pathways in shoot tissues, combining
   data from two alfalfa varieties. (B) shows the distribution of enriched
   metabolic pathways in root tissues, combining data from two alfalfa
   varieties. In the enrichment plots, dot size represents the number of
   genes associated with each pathway, and dot color indicates the level
   of significance (P-value)

Gene co-expression network

   Under atrazine stress, alfalfa demonstrates a multifaceted genetic
   response. Transcriptome analysis of various atrazine-tolerant alfalfa
   varieties identified over 60,000 DEGs. By employing Weighted Gene
   Co-expression Network Analysis (WGCNA), correlations between DEGs and
   atrazine stress-related physiological traits were examined
   (Fig. [96]6A). Twelve distinct gene co-expression modules were
   identified, and module-trait association analyses revealed significant
   correlations between gene expression levels in specific modules and
   particular traits (Fig. [97]6B). For instance, the brown module
   exhibited a strong positive correlation with plant height, with
   correlation coefficients ranging from 0.825. Additionally, it showed
   significant negative correlations with root length and root weight,
   with coefficients ranging from − 0.822 to -0.910. In contrast, the
   turquoise module demonstrated significant negative correlations with
   plant height, shoot biomass, and chlorophyll B content, with
   correlation coefficients ranging from − 0.808 to -0.825. Conversely, it
   exhibited a significant positive correlation with superoxide dismutase
   (SOD), with a correlation coefficient of 0.839 (Fig. [98]6B). These
   findings suggest that genes within the brown module play a crucial role
   in promoting plant height while suppressing root length and root weight
   under atrazine stress. Similarly, in various atrazine-tolerant alfalfa
   varieties, plant height, shoot biomass, and chlorophyll B content were
   strongly negatively correlated with gene expression in the turquoise
   module, while a significant positive correlation with SOD was observed.
   The brown module comprises 12,771 standard genes, whereas the turquoise
   module contains 20,419 standard genes. Among these, the top 30 hub
   genes, identified by connectivity, were recognized as critical genes
   representing the primary functions of these modules (Fig. [99]6C, D).
   Visualization of these hub genes using Cytoscape 3.9.1
   ([100]https://cytoscape.org/) [[101]15] revealed their involvement in
   numerous physiological processes, including photosynthesis,
   respiration, carbon metabolism, cytokinin signaling, stress responses,
   and root and embryo development (Fig. [102]6C, D). These findings
   underscore the vital roles played by these genes under atrazine stress
   conditions.

Fig. 6.

   [103]Fig. 6
   [104]Open in a new tab

   WGCNA reveals gene-physiology correlations in alfalfa under atrazine
   stress. (A) A module classification tree uses colors to represent
   different modules, with grey indicating genes that have not been
   assigned to a specific module. (B) Module-phenotype correlations
   heatmap depicting relationships between modules and physiological
   parameters. (C) Gene Co-Expression Network of DEGs in MEbrown. (D) Gene
   Co-Expression Network of DEGs in MEturquoise

Validation of RNA-seq identified DEGs via qRT-PCR in transcriptional
profiling

   To further validate the expression trends of DEGs, we selected eight
   genes for qRT-PCR experiments. These genes primarily include those
   involved in ethylene synthesis (MS.gene045407, Fig. [105]7A) and
   abscisic acid (ABA) synthesis (MS.gene39848, Fig. [106]7B), both of
   which play crucial roles in regulating plant stress responses.
   Additionally, the ASMT gene (MS.gene77494, Fig. [107]7C) mediates the
   synthesis of melatonin, which scavenges reactive oxygen species (ROS),
   while the APX6 gene (MS.gene67066, Fig. [108]7D) enhances oxidative
   defense through its antioxidant activity. The CHS gene (MS.gene02598,
   Fig. [109]7E) and F3H gene (MS.gene04220, Fig. [110]7F) are responsible
   for synthesizing flavonoids and anthocyanins through phenylpropanoid
   metabolism. The ABCB1 gene (MS.gene005430, Fig. [111]7G) utilizes the
   ATP-binding cassette (ABC) transport system to expel toxins or regulate
   hormone distribution, while the GA2ox gene (MS.gene70788, Fig. [112]7H)
   regulates the degradation of gibberellins. The relative expression
   levels indicate that the qRT-PCR results are consistent with the
   expression trends obtained from RNA sequencing, demonstrating the
   accuracy of our sequencing results and their applicability for
   subsequent transcriptomic analyses. To ensure the reliability of the
   transcriptome data, we validated it by calculating the linear
   regression coefficient between the RNA-seq data (log2(FC)) and the
   qRT-PCR results (ΔΔCt method). The results showed a strong positive
   correlation (R² >0.6) between the two methods, indicating the
   reliability of our transcriptome results. The correlation is
   illustrated in Fig. [113]S3.

Fig. 7.

   [114]Fig. 7
   [115]Open in a new tab

   The results of the qRT-PCR analysis involved eight key genes: (A) ACCO,
   (B) AAO3, (C) ASMT, (D) APX6, (E) CHS, (F) F3H, (G) ABCB1, and (H)
   GA2ox. To ensure reproducibility and reliability, qRT-PCR analyses were
   conducted on three independent biological replicates. The relative
   expression levels were calculated using the 2^−∆∆Ct method, and the
   values presented represent the means of the three biological
   replicates. Duncan’s test was employed to assess statistical
   significance, with different lowercase letters indicating significant
   differences between treatments (P < 0.05)

Discussion

   Atrazine has a soil half-life of up to 57 weeks, and its residues can
   affect alfalfa-maize rotation systems for as long as two years,
   resulting in yield losses that cannot be underestimated [[116]16,
   [117]17]. In our field trials, the root biomass of the sensitive
   variety WL363 was completely suppressed, whereas JN5010 retained 98% of
   its root biomass (Fig. [118]1E). Elucidating the molecular mechanisms
   of atrazine tolerance not only sheds light on the fundamental nature of
   plant herbicide tolerance but also offers practical guidance for
   improving productivity in rotational cropping systems. Exogenous
   stress-induced differential gene expression in plants highlights their
   role in stress response mechanisms. A scatter plot of DEGs demonstrates
   a substantial presence of DEGs in the shoots and roots of both plant
   varieties under atrazine stress. To validate the RNA-seq results and
   identify key pathways involved in atrazine tolerance, we performed GO
   enrichment analysis on the DEGs, validated by qRT-PCR. The root DEGs of
   JN5010 are primarily associated with cell wall functionality and the
   metabolism of proline and cysteine. The cell wall, as the first line of
   defense against external stressors, undergoes functional changes that
   play a critical role in plant adaptation to stress. Proline functions
   as an osmo-protectant, contributing to growth and development while
   scavenging reactive oxygen species; its levels fluctuate under various
   stress conditions [[119]18–[120]24]. Interestingly, the ALDH7A1 gene,
   involved in proline metabolism, showed significant downregulation in
   the roots of JN5010 under atrazine stress. This suggests a potential
   shift in proline metabolism within the roots of the tolerant variety
   under stress. Cysteine supports stress response mechanisms by
   catalyzing H[2]S production, which, as a gaseous signaling molecule,
   mitigates oxidative damage and enhances stress tolerance through
   synergistic effects [[121]25–[122]27]. In contrast, the root DEGs of
   WL363 are primarily enriched in cell cycle-related functions. Research
   indicates that heavy metal stress can trigger oxidative damage, leading
   to cell cycle arrest and subsequent inhibition of plant growth
   [[123]28–[124]30]. Notably, nitric oxide (NO), an essential signaling
   molecule, plays a pivotal role in seed germination, root development,
   and stress responses [[125]31, [126]32]. The biosynthesis and metabolic
   processes of NO are enriched in the root DEGs of both plant varieties.
   Functional analysis of the shoot DEGs in WL363 similarly highlights
   pathways related to the cell cycle and proline metabolism. However, the
   shoot DEGs of JN5010 are enriched in chloroplast- and plastid-related
   functions. Furthermore, distinct, significant metabolic pathways were
   identified in both varieties under atrazine stress. In JN5010, DEGs are
   predominantly enriched in pathways linked to pyruvate metabolism,
   photosynthetic carbon fixation, and peroxisome functions. In WL363,
   shoot DEGs are significantly enriched in pathways associated with
   ascorbate and aldehyde acid metabolism, as well as phenylpropanoid
   biosynthesis, while its root DEGs are enriched in the citric acid cycle
   (TCA cycle) and phosphoinositide metabolism pathways. For example, the
   ABCB1 gene, which encodes an ABC transporter, has been shown in other
   species to confer herbicide tolerance by pumping toxins out of the
   cell. These specific pathways, including ascorbate synthesis, plant
   hormone signal transduction, the MAPK signaling pathway, photosynthetic
   carbon fixation, the Calvin cycle, glycolysis, and ABC transporters,
   offer valuable insights into the variations in tolerance between the
   two varieties. Examining the response patterns of key genes within
   these pathways to atrazine stress offers valuable insights into the
   tolerance differences among various alfalfa varieties.

Ascorbate synthesis pathway

   Ascorbate (AsA) is a vital water-soluble antioxidant in plants, playing
   a central role in various physiological processes, such as
   photosynthesis, respiration, cell division, growth regulation, hormone
   signal integration, and senescence [[127]33–[128]35]. Its homeostasis
   is maintained through mechanisms including synthesis, recycling,
   degradation, and intracellular and intercellular transport. Among
   these, the galactose metabolic pathway constitutes a major route for
   AsA production [[129]36–[130]38]. In this biosynthetic pathway,
   GDP-D-mannose-3’,5’-epimerase (GME) acts as a key mediator in the
   interconversion between GDP-D-mannose, GDP-L-galactose, and
   GDP-L-glucose, while also participating in the L-galactose and L-gulose
   pathways [[131]39, [132]40]. Research shows that GME expression levels
   typically decrease under exogenous stress, and GME inactivation can
   result in reduced cell division and impaired development [[133]41]. In
   this experiment, the expression levels of three genes encoding GME in
   both the shoot parts and roots of WL363 were downregulated under
   atrazine stress (Fig. [134]8). In the roots of JN5010, a similar
   downward trend was observed, but the fold change in downregulation was
   less pronounced than in WL363. Meanwhile, in the shoot parts of JN5010,
   two of the three GME-encoding genes were differentially upregulated,
   with MS.gene44007 showing no significant change (Table. S3). The enzyme
   L-galactose-1-phosphate phosphatase (GPP), encoded by VTC4,
   dephosphorylates L-galactose-1-phosphate to L-galactose, making it an
   essential component of galactose metabolism for AsA synthesis
   [[135]42]. Studies have revealed that mutations in this gene may lead
   to a 50% reduction in ascorbate content [[136]43]. In this experiment,
   the expression of three VTC4-encoding genes in the shoot parts and
   roots of WL363 was downregulated under atrazine stress. In comparison,
   only one gene in the shoot parts and roots of JN5010 exhibited
   downregulation, with this decrease being less pronounced than in WL363.
   Similar expression patterns were observed in genes encoding L-galactose
   dehydrogenase (GalDH), gulonolactone oxidase (Gulo), and
   dehydroascorbate reductase (DHAR) in both varieties (Fig. [137]8).
   L-galactose dehydrogenase catalyzes L-galactonate-1,4-lactone,
   representing the final step in AsA synthesis through the galactose
   pathway. The transcription levels and enzyme activity of L-galactose
   dehydrogenase are positively correlated with L-ascorbate content
   [[138]44–[139]46]. Gulo catalyzes L-gulonolactone-1,4-lactone, a
   critical step in AsA synthesis through the gulose pathway; plants
   genetically transformed to express this enzyme exhibit enhanced
   tolerance to abiotic stress [[140]47]. Dehydroascorbate reductase,
   which is an essential antioxidant enzyme, plays a pivotal role in the
   ascorbate-glutathione (AsA-GSH) cycle, ensuring normal levels of AsA in
   plants. Extensive research indicates that dehydroascorbate reductase
   regulates reactive oxygen species levels, improves leaf photosynthesis,
   and boosts stress tolerance, thereby mitigating plant damage
   [[141]48–[142]51]. Plants’ synthesis and metabolic capacity of AsA
   significantly influence their antioxidant capacity and ability to
   withstand stress [[143]52–[144]54]. For example, research by Mirza
   Hasanuzzaman et al. highlighted that nitrogen supplementation in wheat
   notably increased the activity of AsA and related enzymes, such as SOD,
   CAT, GR, and DHAR, under salt stress, thereby bolstering the
   antioxidant defense system [[145]52]. Moreover, atrazine stress
   significantly induced root SOD activity in JN5010 (142.85 to 200.72
   U/g) but not in WL363 (Fig. [146]1G). While both varieties showed
   increased shoot MDA content under atrazine stress, the increase was
   less pronounced in JN5010 (46.31 to 57.79 nmol/g) compared to WL363
   (70.61 to 84.76 nmol/g) (Fig. [147]1F), which indicated that both
   varieties activate stress response pathways. However, the tolerant
   variety JN5010 exhibited greater root SOD induction and a smaller
   increase in shoot MDA, potentially contributing to its atrazine
   tolerance (Fig. [148]1F, G). In summary, atrazine stress significantly
   suppresses the expression of genes involved in the AsA synthesis
   pathways in the shoot parts and roots of WL363. By contrast, the highly
   tolerant variety JN5010 experiences less inhibition, with certain genes
   even exhibiting upregulated expression in response to atrazine stress.
   This more stable AsA synthesis in JN5010 under atrazine stress may
   serve as a key factor contributing to its robust atrazine tolerance
   (Fig. [149]8).

Fig. 8.

   [150]Fig. 8
   [151]Open in a new tab

   Molecular regulatory network for shoot and root in response to Atrazine
   stress in atrazine-tolerant and sensitive alfalfa varieties

Plant hormone signal transduction and the MAPK signaling pathway

   Ethylene serves as a pivotal regulatory factor in the processes of
   plant growth and development. It plays a crucial role in enhancing
   plant tolerance to various abiotic stresses by synthesizing secondary
   metabolites, interacting with other hormones, and regulating metabolic
   functions. Research has demonstrated that suppressing ethylene
   production and its signal transduction can adversely affect tolerance
   to abiotic stress [[152]55–[153]59]. The ethylene response factor (ERF)
   family, positioned downstream in the ethylene signaling pathway,
   contributes to abiotic stress responses by modulating the expression of
   various functional genes. Extensive studies have revealed that ERF
   genes are influenced by external stressors, with their expression
   closely linked to plant tolerance to conditions such as salinity,
   drought, and heat [[154]60–[155]62]. For instance, research by Zhang et
   al. [[156]63]. found that overexpression of the GmEFR3 gene in soybeans
   increased levels of free proline and soluble carbohydrates, thereby
   improving drought tolerance. Similarly, Wu et al. [[157]64]. showed
   that overexpression of ERF genes in tomatoes enhanced tolerance to
   salt, drought, and cold stresses, highlighting the pivotal role of ERF
   gene expression in stress resilience. In this experiment, two genes
   encoding ERF1 exhibited downregulated expression under atrazine stress,
   with the downregulation in WL363 being more pronounced than in JN5010
   (Fig. [158]8). EIN3, a key component of the ethylene signaling pathway,
   functions downstream of EIN2 and MAPK3/6 and is negatively regulated by
   the ubiquitination/proteasome pathway controlled by EBF1/2. This
   pathway transmits regulatory signals to ethylene response factors
   (ERFs) [[159]65, [160]66]. In the current experiment, EIN3 was
   upregulated in both varieties of alfalfa under atrazine stress;
   however, expression differences in the regulators MPK3 and EBF1/2
   exhibited varietal specificity under atrazine stress (Table. S3). Thus,
   variations in expression levels of the key ERF genes within the
   ethylene signaling pathway under atrazine stress may constitute one of
   the factors influencing atrazine tolerance in alfalfa.

   H₂O₂ is a reactive oxygen species with significant roles in diverse
   biological systems [[161]67]. It not only participates in essential
   developmental processes such as photosynthesis, stomatal movement, and
   leaf senescence but also functions as a pivotal signaling molecule in
   plant growth, development, and stress responses [[162]68–[163]71].
   Studies suggest that increasing endogenous H₂O₂ levels or modulating
   the expression of genes associated with H₂O₂ production can enable
   plants to better adapt to abiotic and biotic stresses [[164]72]. The
   activation of MPK4 by H₂O₂ requires MEKK1, highlighting its central
   role in maintaining H₂O₂ mediated reactive oxygen balance through
   interaction with its downstream target, MPK4. Research has revealed a
   correlation between reactive oxygen species and salicylic acid, with
   mekk1 and mpk4 mutants exhibiting elevated salicylic acid levels
   [[165]73, [166]74]. Moreover, MEKK1 is capable of activating two highly
   homologous MAPKs, MPK4 and MPK6, facilitating signal transduction
   related to reactive oxygen, as well as biotic and abiotic stresses via
   downstream MPK4 and MPK6. The gene expression profiles of MEKK1 and
   MPK4 show considerable overlap in common stress-response genes
   [[167]73, [168]75, [169]76]. Consistent with previous findings, the
   gene expression patterns of MEKK1 and MPK4 in the shoot parts and roots
   of both varieties under atrazine stress were generally similar in this
   experiment. Significant upregulation was observed in WL363, while in
   JN5010, only the expression levels of MPK4 in the roots were
   upregulated, albeit to a lesser extent than in WL363 (Fig. [170]8).
   These observations suggest that WL363 accumulated more H₂O₂ under
   atrazine stress, thereby activating the MEKK1 and MPK4 signaling
   pathways. At the downstream level, MPK4 may also suppress the
   expression of the GST-1 gene, which encodes a key detoxification enzyme
   crucial for eliminating harmful exogenous substances in plants. A
   kinase induced by oxidative signals, particularly hydrogen peroxide,
   OXI1 plays an essential role in transmitting hydrogen peroxide signals
   downstream [[171]77, [172]78]. Stress conditions have been shown to
   induce the gene encoding this kinase, potentially resulting in leaf
   cell death [[173]79, [174]80]. The WRKY transcription factor family
   plays critical roles in plant growth, development, metabolism, and
   responses to environmental stresses, with WRKY22 being a validated
   participant in tolerance mechanisms across multiple plant species
   [[175]81–[176]83]. For example, partial knockout of the CsWRKY22 gene
   increased citrus tolerance to canker disease, while overexpression of
   LlWRKY22 enhanced heat tolerance in lilies, and silencing it produced
   the opposite effect [[177]84, [178]85]. In this experiment, under
   atrazine stress, the expression of the H₂O₂-induced OXI1 gene in
   alfalfa showed an upregulated trend. Meanwhile, the WRKY22 gene at the
   end of the pathway exhibited clear varietal specificity. Of the four
   differential genes encoding WRKY22, two were downregulated in the shoot
   parts and roots of JN5010, whereas no expression differences were
   evident in the shoot parts of WL363, which exhibited marked
   upregulation in the roots (Table. S3). These findings indicate that
   atrazine stress stimulated the signaling pathway from OXI1 to WRKY22,
   with WL363 showing more pronounced changes. The atrazine induced
   increase in H₂O₂ in alfalfa plants activated related signaling
   pathways, including WRKY22-mediated cell death and the inhibition of
   GST-1 gene expression by MPK4. These mechanisms may explain why WL363
   displays greater toxicity under atrazine stress.

   Abscisic acid (ABA) is a vital plant hormone involved in various
   biological processes, playing a key role in plant growth and
   development [[179]86–[180]88]. Extensive research on the mechanisms of
   the ABA signaling pathway under stress conditions has demonstrated that
   ABA enhances plant stress tolerance by inducing the expression of
   specific stress response genes during adverse conditions
   [[181]89–[182]91]. At the core of this signaling pathway lies the
   PYL-PP2C-SnRK2 complex, which operates by interacting with PYL family
   proteins to inhibit the phosphatase activity of PP2Cs in the presence
   of ABA. Conversely, in the absence of ABA, PP2Cs interact with SnRK2s
   through dephosphorylation, leading to the inactivation of SnRK2s, which
   prevents them from phosphorylating downstream response factors. This
   study investigated the expression changes of ABA-related genes under
   atrazine stress, since the plants are more likely to produce ABA under
   stress. The findings revealed that in the shoot parts and roots of
   WL363, the ABA-activated PYL-encoding genes (e.g., MS.gene027877,
   MS.gene031922, MS.gene054122, MS.gene057108, MS.gene29692,
   MS.gene86009) were significantly upregulated, possibly leading to ABA
   insensitivity. In contrast, for JN5010, only MS.gene027877,
   MS.gene031922, and MS.gene29692 were significantly upregulated in the
   shoot parts, while other genes showed an upregulation trend but did not
   reach statistically significant levels (Table. S3). The lower
   upregulation of these PYL genes may result in greater ABA sensitivity.
   Furthermore, most PP2C-encoding genes exhibited a downregulation trend
   under atrazine stress, with the exception of MS.gene64771, which showed
   significant upregulation in WL363, suggesting a possible negative
   feedback mechanism to counteract the effects of ABA. This
   downregulation of PP2C allowed for the significant upregulation of the
   SnRK2 gene in the roots of WL363, resulting in a marked increase in the
   expression of associated genes, and a greater reaction to the ABA.
   Interestingly, the gene encoding MPK1_2 (MS.gene69463) was
   significantly downregulated in both varieties, with statistically
   significant reductions observed in the shoot parts and roots of WL363,
   which may mean that it plays an important role in stress responses.
   Additionally, the katE-encoding gene (MS.gene064838) displayed
   significant upregulation in both varieties across the shoot parts and
   roots (Table. S3), likely to combat stress. In summary, atrazine stress
   stimulates ABA production, with the upregulated PYL genes negatively
   regulating PP2C, thereby activating downstream response factors
   (Fig. [183]8). These findings align with previous studies and offer
   valuable insights into the regulatory mechanisms of the ABA signaling
   pathway. In the roots of the highly tolerant variety JN5010, the GST
   gene (MS.gene21351) was upregulated and acted synergistically with the
   ABCB1 transporter (e.g., MS.gene07361), which was upregulated 4.1-fold,
   enabling efficient extrusion of atrazine into the extracellular space
   or vacuolar compartments. Therefore, this combined detoxification
   mechanism offers two targets, the GSTC domain and the ABCB1 coding
   region, suitable for both marker assisted selection and CRISPR based
   improvement.

Photosynthetic carbon fixation and sugar metabolism

   Photosynthesis is the process whereby plants convert carbon dioxide and
   water into organic matter using light energy while releasing oxygen.
   This process is vital for the survival of living organisms and plays a
   key role in the carbon and oxygen cycles on Earth [[184]92, [185]93].
   Abiotic stress often disrupts photosynthesis by reducing photochemical
   electron transport rates, deactivating photosystem (PS) centers,
   degrading pigments and proteins, and interfering with other biological
   pathways, ultimately affecting normal plant growth [[186]94, [187]95].
   In this experiment, atrazine stress led to a decrease in the expression
   levels of most genes associated with photosynthesis (Table. S3).
   Specifically, for genes related to photosystem I (PS I), MS.gene018588
   and MS.gene32425 were downregulated in the shoot parts of WL363, while
   MS.gene044778 and MS.gene88194 were downregulated in JN5010. Regarding
   PS II genes, MS.gene039151, MS.gene036154, MS.gene040618, and
   MS.gene065561 in WL363 exhibited expression changes under atrazine
   stress, while no significant expression differences were observed in
   the shoot parts of JN5010. For electron transport-related genes,
   atrazine suppressed the expression of MS.gene66754, MS.gene003300,
   MS.gene02098, and MS.gene014378 in WL363, whereas only MS.gene002059
   was upregulated in JN5010 (Table.S3). These findings demonstrate that
   atrazine stress suppressed the expression of most
   photosynthesis-related genes, with expression variations in PS II and
   electron transport genes being more closely linked to atrazine
   tolerance. Gene expression changes were significant in WL363, while
   only a limited number of genes were affected in the highly tolerant
   variety JN5010 (Fig. [188]8). This indicates that tolerance differences
   among alfalfa varieties primarily lie within the PS II and electron
   transport processes, consistent with previous research on atrazine
   toxicity mechanisms [[189]8, [190]96].

   The Calvin cycle is a crucial component of photosynthesis, acting as an
   essential pathway for plants to fix carbon from the environment.
   However, its processes can be disrupted by factors such as reduced
   carbon dioxide concentration or fluctuations in light intensity. Under
   abiotic stress, plants modulate the expression of key enzymes and
   intermediates via multiple pathways, leading to alterations in the
   Calvin cycle. The inhibition of these enzymes is often responsible for
   the stunted growth observed in plants subjected to such stressful
   conditions [[191]97–[192]101]. Ribulose-1,5-bisphosphate
   carboxylase/oxygenase (Rubisco), a pivotal enzyme in the Calvin cycle,
   plays a central role in carbon fixation. Its catalytic activity depends
   on the cooperative function of large (L) and small (S) subunits, which
   collectively enhance its efficiency [[193]102, [194]103].
   Transcriptomic analyses have revealed that under atrazine-induced
   stress, the expression levels of rbcS-encoding genes in the shoot parts
   of both studied plant varieties were upregulated. Conversely, the
   expression of key enzymes-glyceraldehyde-3-phosphate dehydrogenase
   (GAPDH) and phosphoglycerate kinase (PGK)-responsible for the
   conversion of 3-phosphoglycerate into glyceraldehyde-3-phosphate, was
   suppressed (Fig. [195]8). The carbon fixation product,
   glyceraldehyde-3-phosphate, is a precursor for the synthesis of starch
   or fructose-1,6-bisphosphate through gluconeogenesis and glycolysis
   pathways. Notably, the genes encoding aldolase (ALDO) and triose
   phosphate isomerase (TPI), which are integral to these subsequent
   steps, also exhibited a downregulation trend. Numerous studies have
   demonstrated that the Calvin cycle serves as a source of organic matter
   and chemical energy, underpinning high crop yields and stress tolerance
   in plants. However, environmental changes can impair certain enzymes
   and reduce carbon fixation efficiency [[196]104–[197]106]. The findings
   from this experiment indicate that atrazine stress disrupted the carbon
   fixation process in alfalfa, ultimately influencing plant growth and
   development. Sugars, as central molecules in biological systems, play a
   vital role in gene expression, metabolism, the cell cycle, and the
   overall development of both eukaryotes and prokaryotes [[198]107].
   Glycolysis, an essential carbohydrate metabolism pathway, oxidizes
   sugars to produce ATP, pyruvate, and NADH in plants. Under abiotic
   stress, plants regulate enzymatic processes within glycolysis, which
   not only serve as signaling mechanisms but also function as key centers
   for managing energy during stress tolerance [[199]108, [200]109].
   Within the glycolytic pathway, hexokinase acts as a critical
   rate-limiting enzyme in glycolytic metabolism and plant respiration. It
   catalyzes the phosphorylation of fructose to generate
   fructose-6-phosphate, marking the initial step of glycolysis, while
   simultaneously influencing sugar metabolism and storage. Furthermore,
   hexokinase plays an integral role in sugar signaling pathways, with
   research showing its expression is induced under conditions of salt,
   osmotic, and cold stress [[201]109–[202]111]. Fructose-6-phosphate
   kinase and pyruvate kinase also play significant roles in glycolysis by
   consuming ATP to catalyze the conversion of fructose-6-phosphate into
   fructose-1,6-bisphosphate and phosphoenolpyruvate into pyruvate,
   respectively. These reactions are both irreversible and represent key
   regulatory nodes within the glycolytic pathway. Experimental results
   indicated that the expression levels of genes encoding these three
   enzymes were consistently downregulated under atrazine stress. However,
   notable differences were observed between the two varieties studied:
   while hexokinase and fructose-6-phosphate kinase showed reduced
   expression levels in WL363, their downregulation was less pronounced in
   JN5010. On the other hand, pyruvate kinase exhibited a sharper decline
   in expression in JN5010 compared to WL363 (Fig. [203]8). These findings
   highlight that atrazine stress leads to changes in the expression
   patterns of key enzyme-encoding genes in the glycolytic pathways of the
   shoot parts and roots of purple flowering alfalfa, with distinct
   variations between the two plant varieties.

   Trehalose (Tre) is a non-reducing disaccharide composed of two glucose
   molecules. It is typically found in low concentrations across most
   plants but plays a crucial role in various physiological processes,
   including carbohydrate storage, transport, and signaling regulation.
   Numerous studies have demonstrated that trehalose and its precursor,
   trehalose-6-phosphate, serve as vital signaling molecules, regulating
   plant growth and development while contributing to metabolic control
   and gene expression. Furthermore, trehalose is strongly associated with
   plant responses to abiotic stress [[204]112–[205]115], playing a key
   role in supporting normal growth and development under unfavorable
   environmental conditions. Research has shown that the exogenous
   application of trehalose can effectively alleviate high-temperature
   stress during the grain-filling stage in wheat, primarily by decreasing
   reactive oxygen species levels and enhancing antioxidant enzyme
   activity. Additionally, trehalose can mitigate drought stress in sweet
   sorghum, improving the plant’s seedling drought tolerance [[206]98,
   [207]115, [208]116]. Trehalose-6-phosphate synthase (TPS), the enzyme
   critical to trehalose synthesis, has been observed to exhibit
   upregulated gene expression under atrazine stress. Experimental results
   indicated that the gene encoding TPS showed higher levels of
   upregulation in JN5010 compared to WL363 in the aboveground parts of
   the plant (Fig. [209]8). Similarly, most TPS-related genes demonstrated
   more pronounced upregulation in JN5010 roots (Table. S3). These
   findings suggest that both plant varieties deploy trehalose synthesis
   pathways as a response to atrazine stress. However, the
   atrazine-tolerant variety JN5010 exhibited a greater degree of TPS gene
   expression upregulation, potentially driving increased trehalose
   synthesis and enhancing its tolerance to atrazine stress (Fig. [210]8).

Detoxification pathways mediated by ABC transporters and flavonoids

   The plant ABC transporter family comprises eight subtypes of membrane
   transport proteins essential for detoxification processes, facilitating
   the transport of diverse substrates [[211]117]. These proteins utilize
   ATP to expel harmful substances, such as heavy metal ions, organic
   pollutants, secondary metabolites, and xenobiotics, either to the
   extracellular space or vacuoles, thereby mitigating cellular toxicity
   [[212]28, [213]118, [214]119]. Research has revealed the functional
   diversity within the ABC transporter family. The ABCA subfamily is
   primarily responsible for lipid transport, maintenance of membrane
   structure, and stress response [[215]120–[216]122]. For instance, the
   soybean protein GmABCA7, localized in peroxisomes, promotes fatty acid
   β-oxidation and facilitates seed germination in Arabidopsis experiments
   [[217]123]. The ABCB subfamily predominantly regulates and transports
   hormonal substances. In tomatoes, SIABCB demonstrates high expression
   in roots, suggesting its involvement in ion and heavy metal transport
   [[218]123]. Likewise, the ABCC subfamily has been identified to
   transport metabolites, such as anthocyanins and phytic acid, alongside
   detoxification functions [[219]124–[220]126]. Studies in Arabidopsis
   have shown that AtABCC1 and AtABCC2 contribute to arsenic, cadmium, and
   mercury tolerance. Under cadmium stress, the expression of AtABCC6 is
   upregulated, further validating its role in heavy metal detoxification
   [[221]119, [222]127, [223]128]. The ABCG subfamily, the largest branch
   of ABC transporters, is integral to plant detoxification and defense
   mechanisms [[224]121, [225]129]. Research indicates that this subfamily
   enhances tolerance to biotic and abiotic stresses by mediating the
   secretion of secondary metabolites and facilitating the efflux of
   harmful substances. For example, transport proteins such as
   NpPDR1/NpABC1 and NtPDR1 have been shown to promote the secretion of
   antifungal terpenoids and antimicrobial metabolites, playing a pivotal
   role in plants’ chemical defense systems [[226]130, [227]131]. Further
   supporting the role of ABC transporters in stress response,
   transcriptome analysis suggested an upregulation of the GST gene
   (MS.gene21351) in the roots of the highly tolerant variety JN5010 under
   atrazine stress. Subsequent qRT-PCR analysis confirmed a significant
   upregulation of the ABCB1 transporter (MS.gene005430) under the same
   conditions (P < 0.05). The tolerant variety JN5010 appears to establish
   a robust atrazine tolerance mechanism through the coordinated action of
   the ABCA3 gene in the shoot parts and the ABCG2 gene in the underground
   parts, while the sensitive variety WL363 primarily depends on the ABCB1
   gene in the shoot parts.

   Notably, the CHS gene, responsible for synthesizing flavonoids and
   anthocyanins through phenylpropanoid metabolism, was upregulated by
   3.5-fold in both root and shoot tissues of WL363 under atrazine stress
   (Fig. [228]7E). However, its MDA content still increased significantly
   (from 70.61 to 84.76 nmol/g, Fig. [229]1F), suggesting that while
   flavonoid (especially anthocyanin) accumulation may partially alleviate
   membrane lipid peroxidation, it fails to completely inhibit oxidative
   damage. Given that MYBL2 acts as a negative regulator of flavonoid
   biosynthesis, its CRISPR-mediated knockout could further amplify this
   endogenous antioxidant response in WL363 to enhance stress tolerance.
   For the selection and improvement of highly atrazine-tolerant alfalfa
   lines like JN5010, we propose focusing on detoxification pathways.
   Although the GST gene (MS.gene21351) was significantly upregulated in
   JN5010 roots under atrazine stress, the ABCB1 transporter
   (MS.gene005430) showed a slight yet significant increase in expression
   (Fig. [230]7G). Thus, the ABCB1 locus could still be considered as a
   potential molecular marker for preliminary screening. Thereafter,
   CRISPR-Cas9 multiplex editing could be employed to simultaneously
   target ABCA3 and ABCG2 to enhance detoxification pathways in tolerant
   lines, and CHS and MYBL2 to boost antioxidant defenses specifically in
   sensitive lines like WL363.

Conclusion

   This study revealed the differences in gene expression between the
   tolerant variety JN5010 and the sensitive variety WL363 under atrazine
   stress, as well as their potential molecular regulatory mechanisms. The
   results showed that JN5010 exhibits significantly higher gene
   expression stability in both root and shoot tissues compared to WL363,
   indicating that JN5010 displayed stronger atrazine tolerance. Further
   analysis demonstrated that the DEGs in JN5010 were primarily enriched
   in proline metabolism and antioxidant defense pathways, whereas DEGs
   involved in cell cycle regulation and phenylpropanoid biosynthesis
   pathways were enriched in WL363. Despite photosynthesis and sugar
   metabolism were affected under atrazine stress in both alfalfa
   varieties, WL363 showed significantly greater suppression of the
   expression of genes related to photosystem II and electron transport
   compared to JN5010. Moreover, plant hormone signaling and MAPK
   signaling pathways played crucial roles in the response to atrazine
   stress, further supporting the advantages of JN5010 in detoxification
   and metabolism. Notably, the ABC transport proteins in JN5010 were
   significantly upregulated under stress, suggesting their key roles in
   detoxification and substance transport. Differential gene expression
   contributed to JN5010’s functional stability under stress,
   demonstrating better adaptability. These findings enhance our
   understanding of how alfalfa tolerates atrazine stress and provide
   important insights for developing atrazine-tolerant varieties.
   Nevertheless, our future research still needs to further validate
   specific pathways, especially through metabolomic analysis of the
   dynamic changes of proline and flavonoid substances, as well as through
   proteomic methods to verify the activity and content changes of key
   enzymes.

Supplementary Information

   Below is the link to the electronic supplementary material.
   [231]Supplementary Material 1^ (12.4KB, xlsx)
   [232]Supplementary Material 2^ (17.3KB, docx)
   [233]Supplementary Material 3^ (1,004.9KB, docx)
   [234]Supplementary Material 4^ (44.5KB, xls)
   [235]Supplementary Material 5^ (459.9KB, docx)
   [236]Supplementary Material 6^ (275.8KB, docx)

Acknowledgements