The seeds of natural ecotypes of sandrice were collected from Dengkou (DK) County, Dulan (DL) County, and Aerxiang (AEX) Village of northern China in 2014 (Table 1), and consistently preserved in a seed cabinet to ensure their viability at a constant temperature of 20.00℃ and humidity of 5.00%. The three ecotypes have been growing in distinct environments for a long time, resulting in noticeable differences in their phenotypic traits (Yin et al., 2016a). The data on in situ environmental factors for the three ecotypes were downloaded from WorldClim (https://www.worldclim.org/), and displayed along with the phenotypic traits in Table 1. The annual mean precipitation of DK, DL, and AEX ecotypes is 137, 263, and 485 mm, respectively. In 2022, seeds of the three ecotypes were grinded with 1 mm diameter sand (the ratio of sand to seed is 2:1, v/v), and germinated overnight in a dark incubator at 25.00°C. The next day, the germinated seeds were sown in pots filled with sand and nutrient soil (9:1, v/v). After the cotyledons grew, pots were placed under artificially simulated climatic conditions: a photoperiod of 16 h/8 h (light/dark) with approximately 30.00%-50.00% relative humidity and a temperature of 25.00°C-30.00°C in the greenhouse at the Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences. There was only one healthy plant with the same growing consistency maintained in each pot by removing extra plants. A total of 36 pots were divided into six groups (6 pots for each group), including three control groups (CDK, CDL, and CAEX) and three drought treatment groups (DDK, DDL, and DAEX). All plants were timely watered to maintain a suitable soil moisture content until 35 d after germination (vegetative stage).

Water was continuously provided for the control group plants while withholding it from the drought treatment groups for 9 d until a drought stress phenotype with leaf wilting was observed. At this time, above-ground tissue samples were randomly pooled from 6 plants in each group, immediately frozen in liquid nitrogen, and stored in refrigerator at -80.00°C. A total of 18 samples with 3 biological replicates of each treatment were analyzed.

Soil moisture content, leaf relative water content (RWC; %), actual photochemical quantum yield (Fv°/Fm°, a universal indicator of plant stress severity), above-ground biomass (AGB; g), plant height increase (PHI; cm), leaf area (LA; mm²), stem diameter (STD; mm), and basal branch length (BBL; cm) were applied to evaluate drought resistance and physiological difference among the three ecotypes of sandrice. The soil moisture content of control and 3, 6, and 9 d after drought stress was monitored using a soil moisture sensor (Shun Koda TR-6, Beijing, China) at 9:00 am (LST). Fresh weight (FW; g), turgid weight (TW; g), and dry weight (DW; g) of four leaves from each sample were weighed immediately after being isolated from the plants, rehydrated in distilled water at 4.00°C for 24 h until fully turgid, and dried at 70.00°C for 72 h until constant mass, respectively. We calculated the leaf RWC according to the following formula: leaf RWC (%)=(FW-DW)/(TW-DW)×100.00% (Prasad et al., 2011). Fv°/Fm° was measured using MultispeQ v.2.0 software (PhotosynQ LLC Co. Ltd., East Lansing, Michigan, USA). LA and STD of the plants were measured using PhenoAI air (AgriBrain Co. Ltd., Nanjing, China) and a digital vernier caliper (DL91150, deli Co., Ltd., Yuyao, China), respectively. AGB was weighed using a 1/10,000 electronic analytical balance. PHI and BBL were measured with a ruler. Physiological data with 6 biological replicates of each treatment were statistically analyzed by one-way analysis of variance (ANOVA) and t-test at 0.050 and 0.010 levels using IBM SPSS v.17.0 software (SPSS Inc., Armonk, New York, USA). Bar chart was generated using the GraphPad Prism v.8.2.0 software.

Approximately 100 mg of samples were extracted with 2 mL of 70.00% v/v aqueous methanol using ultrasonic waves for 30 min at a time for a total of 2 h, and the mixed extracts were concentrated to nearly dry on a rotary evaporator at 35.00°C under reduced pressure. Before analysis, the residue was dissolved in 200 μL of 50.00% methanol/water solution and transferred to insert-equipped vials.

The sample extracts were analyzed using the ultra-performance liquid chromatography-mass spectroscopy (UPLC-MS) system including alkaloids, amino acids, catechins, flavones, flavone glycosides, phenolic acids, and theaflavins. The chromatographic separation was performed on a high strength silica (HSS) T3 (50.0 mm×2.1 mm, 1.8 μm, Thermo Fisher Scientific Waltham, USA) column using water plus 0.10% acetic acid solution and acetonitrile plus 0.10% acetic acid solution as solvent system at a flow rate of 0.3 mL/min. The gradient program was set as follows: 0.0 min, 90:10 v/v; 2.0 min, 90:10 v/v; 6.0 min, 40:60 v/v; 8.0 min, 40:60 v/v; 8.1 min, 90:10 v/v; and 12.0 min, 90:10 v/v. The column temperature 40.00°C, and the injection volume was 2 μL.

High resolution mass spectrometry (HRMS) data were recorded on a Q Exactive Hybrid Quadrupole-Orbitrap mass spectrometer (Thermo Fisher Scientific, Waltham, USA) equipped with a heated electrospray ionization source utilizing the single ion monitoring MS acquisition methods. The electrospray ionization source parameters were set as follows: spray voltage, -2.8 kV; sheath gas pressure, 40 arbitrary unit (arb); aux gas pressure, 10 arb; sweep gas pressure, 0 arb; capillary temperature, 320.00°C; and aux gas heater temperature, 350.00°C.

Data were acquired on the Q-Exactive using Xcalibur 4.1 (Thermo Scientific, Waltham, USA), and processed using TraceFinder™4.1 Clinical (Thermo Scientific, Waltham, USA). Quantified data were output in Excel software (Microsoft, Redmond, USA). Principal component analysis (PCA) and orthogonal partial least squares-discriminant analysis (OPLS-DA) were performed using the Metware Cloud, a free online platform for data analysis (https://cloud.metware.cn). The differentially accumulated flavonoids (DAFs) are filtered by a variable importance in the projection (VIP) ≥1.00 and P<0.050. Venn and hierarchical cluster heatmap analysis were performed using the OECloud tools at https://cloud.oebiotech.com. The stacking chart was performed using the Wekemo Bioinclod tools at https://bioincloud.tech.

We extracted total RNA from 18 sandrice above-ground tissue samples using the TRIzol reagent (Invitrogen, Carlsbad, USA) according to the manufacturer's protocol. RNA purity and quantification were evaluated using the NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, USA). RNA integrity was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, USA). Then we constructed the libraries using the VAHTS Universal V6 ribonucleic acid sequence (RNA-seq) Library Prep Kit from Vazyme Biotech Co., Ltd., Nanjing, China. The transcriptome sequencing was conducted by OE Biotech Co., Ltd., Shanghai, China.

The libraries were sequenced on a llumina Novaseq 6000 platform and 150 bp paired-end reads were generated. The raw sequencing files of transcriptomic data are now available in the national center for biotechnology information (NCBI) sequence read archive (SRA) database with the accession number PRJNA940574. Raw reads in FASTQ (FASTQ is a text-based format used to represent biological sequences with their quality scores) format were firstly processed using FASTQ and the low-quality reads were removed. The clean reads were mapped to the AEX genome using Hierarchical Indexing for spliced alignment of transcript (HISAT) software. The expression levels of each gene are represented by the fragments per kilobase million‌ (FPKM) value. The FPKM of each gene was calculated and the corresponding read counts were obtained by high-throughout sequencing (HTSeq)-count, a Phyton module for analyzing high-throughput sequencing data. Differential expression analysis was performed using the DESeq2, an R language package for analyzing differentially expressed genes (DEGs) in high-throughput RNA-seq data. For identifying DEGs, P-value<0.050 and |log₂(fold change)|≥1.000 were set as the threshold. The P-value, obtained through a t-test, is used to assess the statistical significance of differences. However, the q-value is the adjusted value of the P-value, and a more stringent statistical measure than the P-value. Fold change represents the degree of change between an initial value and a final value, calculated as the ratio of the final value to the initial value. Log₂ (fold change) is used to measure the difference in gene expression intensity between a sample and its corresponding control sample. Kyoto encyclopedia of genes and genomes (KEGG) enrichment analyses were employed to recognize the biological role and function of DEGs. KEGG enrichment analysis bubble map was produced using the OECloud tools at https://cloud.oebiotech.com.

Correlation analysis heatmap was produced using OriginPro v.2023 software (OriginLab Corporation, Northampton, USA). The network diagram was generated using the Metware Cloud, with a Pearson's correlation coefficient r≥0.800. And its significance of difference at P<0.050 was calculated and screened by SPSS v.17.0 software. We created the KEGG pathway heatmap based on the upstream and downstream relationships of flavonoid synthesis genes and metabolites in the KEGG pathway database.

RNA extracted from CDK, DDK, CDL, DDL, CAEX, and DAEX was used to validate the transcriptome data by qRT-PCR. Three technical replications were performed for each biological replication. Primers were designed using Primer Premier v.5.0 software (Premier Biosoft, Palo Alto, USA) (Table S1). The mean values of AsUBC22 and AsPP2A were used as the internal reference gene (Fang et al., 2023). Target genes were amplified using PerfectStart SYBR qPCR Supermix (TransGen Biotech, Beijing, China) and a Real-Time qPCR System (Mx3000P, Agilent Technologies Co., Ltd., Santa Clara, USA). The 2^−ΔΔCt method was used to analyze the normalized expression of each sample. It represents the fold change in the relative expression level of the target gene compared with the internal reference gene.

After a 9-d drought treatment, remarkable phenotypic characters were observed in the above-ground tissues of the three ecotypes (Fig. 1a), including wilting and curling of leaves, and plants got shorter and collapsed. In the whole experiment, the soil moisture content of the control group maintained 25.00%. The average soil moisture content of the treatment group gradually decreased to 10.70%, 5.20%, and 4.04% on 3, 6, and 9 d after the drought, respectively (Fig. 1b). Under normal culture conditions, AGB was similar in CDK and CDL (Fig. 1c). PHI, LA, STD, and BBL were the highest in CDL (Fig. 1d-g), compared with CAEX and CDK, indicating that DL has greater above-ground biomass. After the 9-d drought treatment, compared with the respective control groups, LA, STD, and BBL were almost not significantly different in the three ecotypes, except for the BBL of AEX. However, AGB significantly decreased by 75.83%, 51.42%, and 52.54%, and PHI decreased by 94.00%, 52.58%, and 72.93% in DDK, DDL, and DAEX ecotypes, respectively. Additionally, under normal moisture supply, Fv°/Fm°, a universal indicator of plant stress severity, was the highest in CDL (0.73) and the lowest in CAEX (0.65). However, compared with the respective control group, Fv°/Fm° significantly decreased by 51.46%, 49.21% and 34.44% in DDK, DDL, and DAEX, respectively (Fig. 1h). In addition, the leaf RWC was not significantly different among the control groups of the three ecotypes, while there was a significant decrease in the drought treatment group, compared with the respective control group, and leaf RWC was higher in DDL than in DAEX and DDK after 9 d treatment (Fig. 1i). In a word, DK and AEX were presented drought sensitive, while DL showed the strongest resistance to drought stress.

The composition of 15 flavonoid-targeted metabolites and their contents are shown in Table 2, with their contents ranging from 0.01 to 2859.34 ng/100 mg FW. The total amount of these compounds varied from 3560.70 to 5225.33 ng/100 mg FW. Among them, the total flavonoid content in AEX was the highest under both CAEX and DAEX. The PCA of the metabolome data explained 37.40% (principal component; PC1) and 21.90% (PC2) of the variation (Fig. S1a). The three replicates of each treatment group were all within the 95.00% confidence interval. However, drought treatment and control of each ecotype were separated, implying that sandrice responded to drought stress by regulating flavonoid accumulation. As shown in Figure 2a, rutin, isoquercitrin, and astragalin were identified as the most abundant flavonoids in all groups, and these metabolites accounted for more than 95.00% of the total 15 flavonoids. The content of other metabolites was relatively low (Fig. 2a). The OPLS-DA showed q>0.900 for each ecotype (Fig. S1b-d), indicating that these models were reliable and stable and could be used to further screen DAFs. There was a total of 12 DAFs, including 7 DAFs (rutin, quercetin, dihydrokaempferol, dihydroquercetin, naringenin chalcone, astragalin, and apigenin) specific to individual ecotype, 4 DAFs (vitexin, isorhamnetin, epicatechin, and naringenin) shared by the two ecotypes, and 1 DAF (luteolin) shared by the three ecotypes (Fig. 2b). Rutin, the most abundant compound, showed a gradual increase with precipitation gradient (Table 2), and significant increase in AEX after drought stress (Fig. 2c). Isoquercitrin and astragalin showed decrease in DK and DL, but increase in AEX under drought stress. In DL, astragalin showed a significant decrease (Fig. 2c). Those results demonstrated the flavonoids accumulation of the three ecotypes was affected differently, and rutin might play an important role under drought stress.

Total of 128.55 Gb raw reads from 18 samples were generated. After the removal of the low-quality reads, a total of 120.44 Gb clean data were retained for subsequent analyses. The total mapping rate was 95.68%-97.79% of each sample, and 92.19%-94.72% clean reads were uniquely mapped. Q30 (the percentage of bases with a mass value greater than or equal to 30) bases ranged from 90.54% to 91.92%, and the average GC (total proportion of guanine (G) and cytosine (C) bases in a DNA sequence) content was 44.69% (Table S2). The PCA of 18 samples was divided into the three ecotypes along PC 1, and control and drought treatments of each ecotype were separated along PC 2 (Fig. S2). CDK and CDL were closely grouped and were separated after drought treatment, indicating the response mechanisms to drought stress might be divergent.

DEGs analysis was performed between the control and drought treatments of the three ecotypes. A total of 3554, 1557, and 2374 up-regulated DEGs were identified in DK, DL, and AEX, respectively, while 4448, 2753, and 3594 down-regulated DEGs were identified (Fig. 3a). In the three ecotypes, the number of down-regulated DEGs was higher than that of up-regulated DEGs, and the total number of DEGs was the lowest in DL, indicating that DL was more stable than DK and AEX under drought stress. The Venn diagram showed 2584, 383, and 808 DEGs were unique to DK, DL, and AEX, respectively (Fig. 3b). Among the common 3177 DEGs of the three ecotypes, there were 1118 up-regulated and 2020 down-regulated DEGs (Fig. S3a-b). Clustering heatmap of the shared DEGs visually presented the difference between drought and control groups (Fig. 3c). KEGG enrichment analysis bubble map showed the top 20 metabolic pathways of these common DEGs (Fig. 3d), including peroxisome, plant hormone signal transduction, starch and sucrose metabolism, glyoxylate and dicarboxylate metabolism, carbon fixation in photosynthetic organisms, and photosynthesis, indicating that the three ecotypes respond to a drought stress environment by producing plant hormones, regulating cellular processes and metabolic pathways, and synthesizing metabolites. To be mentioned, we noticed that flavonoids biosynthesis pathway did not show a common trend to respond to drought treatment, however, the bubble map of KEGG in different ecotypes showed that flavonoids biosynthesis was triggered by drought in DL (Fig. S4).

To further study the difference in flavonoid synthesis between different ecotypes, we found 14 DEGs in the three ecotypes related to flavonoid biosynthesis (Fig. 4a), with 8 DEGs being down-regulated under drought stress. To our surprise, the expression of some genes was different among the three ecotypes. 3'GT (anthocyanin 3'-O-beta-glucosyltransferase) was up-regulated in AEX, but down-regulated in both DK and DL, while COMT5 (caffeic acid 3-O-methyltransferase) was down-regulated in AEX and up-regulated in DK and DL. These results demonstrated that 3'GT and COMT5 had expression specificity and complementary trends among different genetic lineages.

Correlation analysis between DEGs and DAFs under drought stress showed that astragalin was significantly decreased in DL, which was negatively correlated with CHI (chalcone isomerase), COMT5, and COMT6, while positively correlated with LDOX (leucoanthocyanidin dioxygenase), 3′GT, CYP75B1 (flavonoid 3'-monooxygenase), CHS (chalcone synthase), and ANR (anthocyanidin reductase) (Fig. 4c). Rutin had a significant accumulation in AEX, which was negatively correlated with AOMT2 (flavonoid 3′,5′-methyltransferase), CYP75B1, COMT5, and ANR, while positively correlated with 3′GT, FOMT (flavonoid 3'-O-methyltransferase), CHI, COMT1, and COMT6 (Fig. 4d). When the DAFs and DEGs in flavonoid biosynthesis were mapped to the KEGG metabolism pathways, we found that there were three different regulation models (DK, DL, and AEX) of flavonoid synthesis responding to drought stress (Fig. 5). Combining with Figure 4, 4 critical DEGs in flavonoid biosynthesis pathway shared among the three ecotypes, including CHI, FLS (flavonol synthase), CYP75B1, and ANR. To be mentioned, the flavonoid biosynthesis showed diverse among ecotypes under drought stress, except for the shared DAFs, i.e., vitexin. In DK, dihydrokaempferol and dihydroquercetin significantly increased, and the route dihydrokaempferol-dihydroquercetin activated by CYP75B1 might enhance the adaptation to drought stress (Fig. 5a). In addition, apigenin was significantly positive regulated by CHI (P<0.050) in DL (Fig. 5b). This ecotype might activate the metabolic flux route naringeninapigenin-luteolin to respond to drought stress. In AEX, we found rutin and quercetin were positively regulated by CHI (P<0.050), whose metabolic flux redirected from route quercetin-isorhamnetin to quercetin-isoquercitrin-rutin accompanied by transcriptional suppression on OMT (O-methyltransferase) (Fig. 5c). However, in DK and AEX, isorhamnetin was oppositely regulated under drought stress, which significantly increased in DK and decreased in AEX under drought stress.

To further confirm the reliability and repeatability of the transcriptome data, we subjected 8 structural genes involved in the biosynthetic pathway of flavonoid to qRT-PCR analysis. The results showed that the relative expression of 7 genes was consistent with the corresponding FPKM values in the transcriptomic data (Fig. 6).

The morphological differences among sandrice ecotypes displayed significant local adaptation to precipitation (Yin et al., 2016a, b). In our study, we observed that the seedlings of sandrice exhibited the strong drought tolerance (4.00% soil moisture content) with short stature, collapsed structures, yellowing leaves, and wilting and curling after drought stress (Fig. 1a), similar to phenotypic responses observed in other species (Hu et al., 2021, 2022; Huang et al., 2023). Meanwhile, according to the fluorescence parameters (Fv'/Fm', a universal indicator of plant stress severity) (Sharma et al., 2015), we found that the departure of sandrice from its original habitat and its cultivation under laboratory conditions can be considered a form of stress due to inadequate light exposure. However, throughout the experiment, all conditions, except for the water content, were kept consistent between control and drought treatment groups. This consistency ensured the integrity and validity of drought treatment. Notably, sandrice in DL had the highest in PHI (Fig. 1d), LA (Fig. 1e), Fv'/Fm' (Fig. 1h), and leaf RWC (Fig. 1i) under drought stress, compared with DK and AEX. Conversely, sandrice in DK exhibited the lowest leaf RWC under drought conditions. Leaf RWC is one of the optimum indicators to reveal the water retention capacity of plants under drought-related stress, which is positively correlated with the drought resistance of plants (Karimpour, 2019). Within a plant species, in addition, the leaf RWC of drought sensitive genotype was generally lower than that of resistant genotype (Rampino et al., 2006). Our physiological results indicated that DL displayed greater drought resistance, while DK was more sensitive when facing the same drought environment. DL, the highest-altitude ecotype from Qinghai-Xizang Plateau, due to long-term adaptation to multiple stresses such as cold, ultraviolet B (UV-B), and drought, may develop stronger self-protection mechanisms responding to drought stress. Investigating the precipitation of native habitats, we found the annual mean precipitation of DL (263 mm) is twice that of DK (137 mm). In terms of the precipitation gradient, DK, the arid land ecotype, should be the most drought-tolerant ecotype. In addition, previous data showed that under drought conditions, the arid land ecotype exhibited a greater capacity for adaptation compared with the semi-arid land ecotype, further supporting this finding (Fang et al., 2023). However, based on the phenotypic changes of DK under drought treatment, it appeared that DK could be a drought-sensitive ecotype. It is reasonable to consider the original habitat of DK, located near the Yellow River with abundant groundwater. The physiological traits of DK include short, broad leaf (Table 1), and a fragile photosynthetic system (Fig. 1h), which may be not detrimental to the plant's ability to adapt to drought stress, especially when compared with the sandrice ecotypes with above-ground characteristics such as smaller, thicker leaves, and more robust photoprotective systems (Fang and Xiong, 2015; Bhusal et al., 2021). Thus, to elucidate the molecular mechanism of sandrice adaptation to drought stress, it is essential to conduct more detailed studies on root development under drought treatment across more ecotypes.

The analysis of 15 flavonoid-targeting metabolites revealed that rutin, isoquercitrin, and astragalin were highly enriched in the above-ground tissues of sandrice, collectively accounting for more than 95.00% of the total 15 flavonoids (Fig. 2a). Numerous studies have demonstrated that rutin, isoquercitrin, and astragalin have a number of pharmacological effects, such as antimicrobial effects, anti-allergic activity, antioxidants, cardiovascular disorders, diabetes, cardioprotective properties, neuroprotective effects, anti-cancer, anti-inflammatory, and allergic reactions (Valentova et al., 2014; Ganeshpurkar and Saluja, 2017; Riaz et al., 2018; Resham et al., 2020; Negahdari et al., 2021). These results further supported the finding that above-ground tissues of sandrice contain large amounts of flavonoids with outstanding medicinal value under current culture conditions. Rutin, the most abundant compound, was 1231.57-2859.34 ng/100 mg FW (equivalent to 61.60-143.00 μg/g dry weight (DW)) (Table 2). This result contradicted with a previous study that found the content of isorhamnetin was the highest (557.50 μg/g DW) in sandrice from the common garden trial (Zhou et al., 2021a). The common garden trial was conducted outdoors on plots with relatively high sand content in Wuwei City, Gansu Province, China. Throughout the entire growth period, there was ample sunlight, and environmental factors such as temperature and humidity fluctuated greatly. However, the present study utilized pure sand from the field mixed with nutrient soil, and the experiment was conducted under artificially simulated climatic conditions. During the cultivation period, environmental factors like temperature and humidity remained stable. Although the light intensity was lower than natural sunlight, it was still sufficient to meet the requirements of drought treatment experiments. In summary, the common garden trial and the controlled experiments differed significantly in terms of planting soil, nutritional status, climate, and other factors. Therefore, we propose that the variations in cultivation conditions may be the primary factors contributing to the differences in flavonoid accumulation between the common garden trial and the ambient drought treatment experiment. In fact, this kind of phenomenon has been observed in other plants as well, where number and content of plant flavonoids are easily influenced by their growing environment (Jaakola and Hohtola, 2010; Shi et al., 2022).

In the present study, the content of rutin in the three ecotypes increased under drought stress, indicating its accumulation was conducive to the response to drought stress (Fig. 2a). A previous study has shown rutin accumulation in the cytosol scavenges hydroxyl radicals produced in response to salinity stress in quinoa (Chenopodium quinoa Willd.), which belongs to the same family as sandrice (Ismail et al., 2015). Furthermore, research on Fagopyrum tataricum (L.) Gaertn. also confirmed that the high expression of encoded rutin biosynthetic enzymes and the accumulation of rutin contributed to improving tolerance to abiotic stress (Zhang et al., 2017). Those results suggested that the accumulation of rutin might contribute to scavenging ROS and enhance sandrice adaptability to drought stress. In addition, the common garden experiments also demonstrated the high content of rutin (Zhou et al., 2021b; Fang et al., 2022), indicating that rutin might be involved into the adaptation to the extreme environment in sandrice. In future work, it will still need an in-depth experiment such as applying exogenous rutin to leaves to further confirm whether it enhances drought tolerance in sandrice.

The secondary metabolites accumulated in sandrice may also exhibit diversity in response to stressful environments. Multi-omics integration analysis of flavonoid synthesis pathways showed the diversity of flavonoid accumulation among the three ecotypes of sandrice when responding to drought stress. DK might activate the metabolic flux route dihydrokaempferol-dihydroquercetin in adapting to drought stress (Fig. 5a). DL might activate the metabolic flux route naringenin- apigenin-luteolin in responding to drought stress (Fig. 5b). AEX might redirect metabolic flux from route quercetin-isorhamnetin to route quercetin-isoquercitrin-rutin by transcriptional suppression on OMT (Fig. 5c). Metabolic flux redirection is a common response in plants (Dong et al., 2020), as observed in our previous study on sandrice when exposed to in situ environmental differentiation among ecotypes (Fang et al., 2022). The reason of redirection in sandrice remains unclear, but the unusually high accumulation of some flavonoids might readjust the homeostasis to cope with drought stress. Studies in Achillea pachycephala Rech.f. (Gharibi et al., 2019) and Chrysanthemum morifolium (L.) Ramat. (Hodaei et al., 2018) have shown that water stress can stimulate the synthesis of apigenin and luteolin, accompanied with the up-regulation of their biosynthesis-related genes (CHI, CHS, and F3H) to enhance the ability to adapt to drought stress. Some researches have reported that quercetin is a key intermediate metabolite of flavonoid biosynthesis pathway to enhance drought tolerance (Yeloojeh et al., 2020; Singh et al., 2021), and this kind of phenomenon was also found in the ecotype AEX of sandrice under drought treatment.

Additionally, the accumulation of flavonoids in sandrice also exhibited variability across different ecotypes under cold stress. Specifically, our recent research revealed that AEX and DL demonstrated a higher phytochemical diversity when coping with low temperature, compared with DK (Zhao et al., 2024). This, together with the drought treatment results, suggest that the diversity of flavonoid accumulation patterns is ubiquitous among sandrice ecotypes with diverse environmental conditions. Although the flavonoid accumulation observed in this research was a consequence of genotype and environment interactions (Hodaei et al., 2018), the significant genetic differentiation was found between central lineage (DK and DL) and eastern lineage (AEX) in sandrice, which were also geographically differed in annual mean precipitation (Qian et al., 2016, 2021). However, the upstream part of the biosynthesis pathway of flavonoids was triggered by drought in DL, differed to DK that belongs to the same genetic lineage. Considering to environmental characteristics of DL, which came from Qinghai-Xizang Plateau, plants from DL endures not only drought but also the other stressors such as low temperature and UV-B. Thus, we proposed that the downstream part of the flavonoid biosynthesis could be more or less involved into adaptation to UV-B radiation and/or cold temperature. Actually, the downstream flavonoids such as dihydroquercetin, isorhamnetin, kaempferol, and quercetin, were witnessed to be significantly accumulated under cold treatment in DL (Zhao et al., 2024). In this case, we would conclude that the phytochemistry diversity of flavonoids among sandrice ecotypes could be the consequences of local adaptation to the combination of multiple stresses across different deserts. In the future, it is essential to undertake a series of combined stress experiments to delve deeper into the molecular metabolic mechanisms of sandrice adaptation to multiple stresses (Jansen et al., 2021; Jan et al., 2022).