The wheat (Triticum aestivum) genome (wheat genome version IWGSCV1.1) and protein sequences were obtained from the Ensembl Plants website (http://plants.ensembl.org/index.html accessed on 11 March 2022). HMMER 3.0 (hidden Markov model, HMM) (http://hmmer.org/) software, BLASTP (https://blast.ncbi.nlm.nih.gov/) method and Pfam (http://pfam.xfam.org/) database searches were used to identify BES1 genes in Triticum aestivum, Aegilops tauschii, Triticum dicoccoides, Triticum turgidum, Triticum urartu, and Hordeum vulgare using reported Arabidopsis, rice, and maize BES1 genes as queries. BES1 proteins were verified using HMMER software (E value < 0.01) and the NCBI conserved domain database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). The HMM profile of the BES1 domain was downloaded from the Pfam (PF05687) and WheatOmics (http://wheatomics.sdau.edu.cn/tools/proteinfamily.html) databases. The physical and chemical properties of TaBES1 proteins from the longest transcripts were determined using the ExPASy website (http://web.ExPASy.org/translate/).

MEGA 10.0 software was used to reconstruct a phylogenetic tree of the TaBES1 family using the maximum-likelihood method with the full-length protein sequences. The BES1-like proteins from eight species were divided into subfamilies according to their clustering patterns. Each TaBES1 gene was mapped onto one of the 21 wheat chromosomes according to the IWGSC RefSeq v1.1 (‘Chinese Spring’) reference genome of the Triticeae Multi-omics Center (http://wheatomics.sdau.edu.cn/). The Ensembl database (http://plants.ensembl.org/biomart/martview/) was used to extract information about chromosome length, and a physical map was drawn using MG2C software v2.1 (http://mg2c.iask.in/mg2c-v2.1/).

Gene structures were visualized using the Gene Structure Display Server (http://gsds.gao-lab.org/). Conserved motifs of TaBES1 proteins were analyzed using MEME (http://meme-suite.org/index.html). TBtools software (https://github.com/CJ-Chen/TBtools/releases) was then used to analyze and visualize the TaBES1 gene structures and conserved motifs of TaBES1 proteins. A 2-kb upstream promoter region for each TaBES1 gene was retrieved from the Ensembl Plants website and then analyzed using PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/), finally the prediction of cis-regulatory elements was visualized using TBtools software. Alignment of BES1-domain sequences was performed using MEGA 10 software, and the sequence logo of this domain was created using the weblogo website (http://weblogo.berkeley.edu/).

Sequences of TaBES1 genes were uploaded to the psRNATarget website (https://www.zhaolab.org/psRNATarget/) to search for potential miRNA-mediated regulation based on known miRNAs and gene sequence information. The psRobot website (http://omicslab.genetics.ac.cn/) was used for the analysis of small RNA target predictions, and putative regulatory networks were visualized using Cytoscape software v3.8.2 (https://cytoscape.org/). Expression data based on transcriptome deep sequencing (RNA-seq) and microarray databases of T. aestivum ‘Chinese Spring’ from the WheatOmics database (http://wheatomics.sdau.edu.cn/). To analyze tissue-specific expression patterns, the coding sequence of each gene was submitted to the WheatOmics database to obtain expression data for roots, stems, leaves, spikes, anthers, and grain tissues. Transcripts per million (TPM) values for TaBES1 genes under drought, salt, cold, and heat stress conditions were also retrieved from the Triticeae Multi-omics Center database and visualized as heatmaps using TBtools software.

To investigate the tissue-specific expression profiles of TaBES1 genes, the following tissues were collected: root and shoot tissues at the seedling stage, leaves at the tilling stage, spikes and anthers at the full heading stage, and grain at the milking stage, from the wheat cultivar ‘Xiaobaimai’ (XBM), which serves as the original backbone parent wheat variety in northern China. As XBM is known for its sensitivity to abiotic stress at seedling stage (Li L. et al., 2018; Cao et al., 2021), its seeds were placed in Petri dishes containing two pieces of filter paper soaked with water at 23°C under a 16 h light/8 h dark photoperiod; germinated seeds were then transferred to plastic mesh grids at 28°C under a 16 h light/8-h dark photoperiod. Fourteen-day-old XBM seedlings were exposed to drought (25%, [w/v] polyethylene glycol 6000 [PEG-6000]), salt (250 mM NaCl), cold (4°C), or heat (40°C) for 0, 1, 2, 5, 10, or 24 h (Liu et al., 2015; Xu et al., 2020), and they were used for the gene expression analysis.

Seedlings from the cultivars ‘Taiyuan806’ (TY806), ‘Zhongyou9507’ (ZY9507), ‘Yumai8’ (YM8), ‘Xinong318’ (XN318), and ‘Zhongyin6’ (ZY6) were used to test drought resistance. TY806, XN318, and ZY6 were reported to be drought-tolerant, while XBM, ZY9507, and YM8 were considered drought-sensitive cultivars (Li L. et al., 2018; Cao et al., 2021). They were planted in pots containing a mixture of vermiculite and humus soil in equal volumes, and grown under appropriate watering conditions at 70% relative humidity and a temperature of 25°C, with a photoperiod of 16 hours of light and 8 hours of darkness in a climate incubator. When the seedlings reached the 2-leaf stage, watering was withheld for 2 weeks to induce drought stress. Additionally, PEG-induced drought stress was applied to these six wheat varieties. The seeds of the six wheat cultivars were germinated in a dark environment at a temperature of 23°C. After a two-day germination period, uniform seedlings were selected and transplanted into plastic containers filled with a hydroponic solution based on the Hoagland culture method. The growth conditions were maintained at a temperature of 28°C, with a photoperiod of 16 hours of light and 8 hours of darkness. Subsequently, ten-day-old seedlings with two leaves were subjected to PEG-induced drought stress (PEG-6000, 25% w/v) for 24 hours. To assess the impact of drought stress, the superoxide dismutase (SOD) activity, malondialdehyde (MDA) content, and peroxidase (POD) activity of the entire seedlings were measured. These measurements were conducted using the Total Superoxide Dismutase assay kit, Plant Malondialdehyde assay kit, and Peroxidase assay kit, respectively (Jiancheng Bioengineering Institute located in Nanjing, China). All assays were performed according to the manufacturer’s instructions. The expression levels of the other five wheat cultivars were assessed using 14-day-old seedlings under 0, 1, 2, 5, 10, and 24 hours of PEG treatment, which was the same as the treatment given to XBM. Finally, the seedlings of all six wheat cultivars under PEG treatment were used for gene expression analysis. Three independent experiments were performed for each line (Li L. et al., 2018; Cao et al., 2021)

Temperature was found to play a crucial role in controlling the fertility of TGMS lines, with moderate temperature promoting fertility and low temperature inducing sterility during the anther development (Liu et al., 2015; Xu et al., 2020). In this study, TGMS lines BS1453 and BS366 were grown in Beijing (BJ) (moderate temperature fertile environment) and Nanyang (NY) of Henan province (low temperature sterile environment). The common wheat cultivar ‘Jing411’ (J411) was used as a control and exhibited a fertile phenotype in both the BJ and NY environments. Anthers from BJ and NY environments were collected and frozen in liquid nitrogen immediately at each male reproductive development stage of central callose (S6), meiotic (S7), tetrad (S8), and young microspore (S9) (Browne et al., 2018; Xu et al., 2020). After confirming the anther development stages by microscope, anthers at the four developmental stages from BS1453, BS366 and J411 were performed for the following gene expression analysis by qPCR. The pollen starch was subjected to staining using a 10% (v/v) Lugol solution (Lee et al., 2016; Xu et al., 2020). The RNAs of 12 samples (anthers at S9 stage of BS1453, BS366 and J411 under fertile and sterile conditions, two biological replicates per stage of three cultivars) were subjected to 150 bp paried-end sequencing using Novaseq 6000 platform (Illumina), sequencing were performed according to the manufacturer’s standard protocol. Illumina RNA seq was constructed using the method of Külahoglu et al. (Browne et al., 2018; Xu et al., 2020). Each gene’s transcript profile was evaluated using fragments per kilobase of exon per million fragments mapped (FPKM).

Seedlings under various abiotic stress conditions and anthers from sterile and fertile conditions were collected and stored at -80°C. Total RNA was isolated from the samples using TRIzol reagent (Invitrogen) (Pandey et al., 2022). For qPCR, first-strand cDNA was synthesized using a Primer Script RT reagent Kit with gDNA Eraser (TaKaRa). The transcript levels were determined with ChamQ SYBR qPCR Master Mix (Vazyme, Nanjing, China), and the qPCR were conducted in a CFXTM qPCR detection system (Bio-Rad, California, USA). Reaction systems were prepared in 20 μL volumes as follows: 10μL of 2μL Taq Pro Universal SYBR qPCR Master Mix, 6 μL of RNase-free water, 2 μL of cDNA, 1 μL of forward primer, and 1 μL of reverse primer. The primers were designed using Primer Premier 5.0. (Table S2). The qPCR reaction conditions were 94°C for 30 s, followed by 40 cycles of 95°C for 5s, 60°C for 30 s. Then a melting curve of 95°C for 15 s, 60°C for 60 s, and 95°C for 15 s was used. TaACTIN gene (Gene ID: LOC542814) was performed as the reference gene and the relative gene expression levels were calculated using the 2^-ΔΔCt method. Three biological replicates were performed for each treatment, with three technical replicates per sample.

Utilizing bioinformatics tools such as HMMER analysis and searches with the BLASTP tool and the Pfam database, we identified 59 putative BES1-like genes from the genomes of common wheat and its relatives. Of these genes, 23 were in Triticum aestivum, four in red wild einkorn wheat (Triticum urartu L.), 13 in wild emmer wheat (Triticum dicoccoides L.), 6 in Tausch’s goatgrass (Aegilops tauschii L.), and thirteen in durum wheat (Triticum turgidum L.), (Table S1). According to the gene nomenclature, we named the 23 wheat BES1-like genes based on their chromosomal locations and homoeologous genes was named in the same number with different alphabet letters corresponding to their respective subgenomes (Table S1). TaBES1 genes were unevenly distributed across 17 of 21 wheat chromosomes (Figure S1). Each chromosome harbored one TaBES1 gene except chromosomes 3A, 3B, 3D, 6A, and 6B, which each carried two genes. No TaBES1 gene mapped to chromosomes 1B, 1D, 5A, or 5D. The 23 TaBES1 genes were almost evenly distributed across the A (8), B (8), and D (6) subgenomes. As listed in Table S1, the average amino acid length of the putative TaBES1 proteins was 365 amino acids (aa), ranging from 178 to 686 aa, and TaBES1-4D was the longest protein (686 aa). The predicted molecular weights of TaBES1 proteins ranged from 19.27 kDa (TaBES1-3D1) to 75.47 kDa (TaBES1-4D); all TaBES1 proteins were hydrophobic.

We constructed an unrooted phylogenetic tree using the full-length sequences for the 59 BES1 proteins from common wheat and its relatives identified in this study and those of reported BES1 proteins from Arabidopsis (eight), rice (six), and maize (eleven) (Table S3). The tree formed four major clades (I-IV) based on the phylogenetic relationship between each BES1-like protein and Arabidopsis BES1 members (Figure 1). TaBES1-4B, TaBES1-4D, TaBES1-6A1, TaBES1-6B1, and TaBES1-Un were closely related to AtBMY2 and AtBMY4, which clustered in clade I. TaBES1-3A2, TaBES1-3B2, and TaBES1-3D2 formed a tight clade with Arabidopsis BEH3 and BEH4 within clade II. Clade III contained TaBES1-6A2, TaBES1-6B2, TaBES1-6D, TaBES1-7A, TaBES1-7B, and TaBES1-7D, with no Arabidopsis member. Finally, TaBES1-1A, TaBES1-2A, TaBES1-2B, TaBES1-2D, TaBES1-3A1, TaBES1-3B1, TaBES1-3D1, TaBES1-4A, and TaBES1-5B were closely related to Arabidopsis tapetum development associated proteins BEH1, BEH2, BES1, and BZR1, which all belonged to the largest clade IV.

We analyzed the exon-intron structure of all wheat BES1-like genes by aligning their predicted coding sequences to the wheat genome. TaBES1 genes from the same subclade shared a common gene structure (Figure S2A). Most wheat BES1 genes contained one intron and two exons; TaBES1-4A had two introns and three exons, while TaBES1-4B, TaBES1-4D, TaBES1-6A1, TaBES1-6B1, and TaBES1-Un belonging to clade I had more than seven introns. The intron and exon number of clade I was far greater than that of the other clades.

We also examined cis-regulatory elements within 2-kb fragments of promoter sequence upstream of each TaBES1 gene. As shown in Figure S2B, we identified cis-acting elements related to phytohormone signaling, including methyl jasmonate, abscisic acid, salicylic acid, gibberellin, and auxin. We also detected several abiotic stress responsive cis-acting elements, such as those associated with wounding, drought, cold, and light.

Proteins belonging to the same subfamily and sharing similar motif compositions are likely to share similar functions, prompting us to characterize functional domains in each TaBES1 protein (Li et al., 2019b). Motif analysis highlighted the high level of conservation of TaBES1 proteins, with motif 1 being shared by all TaBES1 proteins except TaBES1-5B (Figure S2C). Moreover, we detected motif 2 in all proteins except TaBES1-6A1, while motif 5 was present in all proteins except TaBES1-3A1, TaBES1-3B1, and TaBES1-3D1. Closely related family members in the phylogenetic tree shared a common motif composition (Li et al., 2019b); we identified motif 10 in members from subclades I, II, and III, while motifs 3, 8, and 9 only existed within subclade I members (Figure S2C).

The BES1-type domain is conserved across plant species. An analysis of conserved amino acid residues among wheat BES1-type domains produced results similar to comparable analyses in Arabidopsis, maize, tomato, and other species, illustrating the high conservation of the N terminus, while the C terminus of TaBES1 members was less conserved (Figure S3A). BES1 proteins contain a putative nuclear localization sequence, a highly conserved N-terminal region, a BIN2 (BR-INSENSITIVE 2) phosphorylation domain, a PEST motif, and a C-terminal structure (Cao et al., 2020). Similar to other species, sequence logo analysis showed that the N-terminal BES1-type domain was highly conserved between Arabidopsis and wheat BES1 family members (Figure S3B).

To investigate the potential regulation of TaBES1 genes by small regulatory RNAs, we uploaded all 23 TaBES1 sequences to the wheat miRNA database; which predicted that eight wheat miRNAs have the capacity to regulate TaBES1 transcript levels, based on their high sequence complementarity (Figure 2). Most TaBES1 genes were regulated by only one miRNA, although TaBES1-2A, TaBES1-2D, TaBES1-3B2, TaBES1-6A2, TaBES1-6B2, and TaBES1-6D were each regulated by two miRNAs. Notably, Tae-miRNA9780 targeted the transcripts of five TaBES1 genes (TaBES1-2A, TaBES1-2D, TaBES1-6A2, TaBES1-6B2, and TaBES1-6D), suggesting its importance during BR signal transduction. The miR171, miR1119, miR9657, miR9776, and miR9780 families have been reported to play roles in plant stress responses, with the miR171, miR9657, and miR9677 families being also associated with anther development (Bai et al., 2017; Shi et al., 2018; Sun et al., 2018; Li et al., 2019a; Zeeshan et al., 2021; Zhang et al., 2022).

We obtained initial tissue expression data from the WheatOmics database. The expression levels of TaBES1 genes varied across tissues and developmental stages (Figure S4). Expression levels of nine of the 23 TaBES1 family members were extremely low during all developmental stages and in all tissues. By contrast, TaBES1-3A2, TaBES1-3B2, and TaBES1-3D2 on chromosome 3 showed higher expression levels during the whole life cycle (Figure S4). We confirmed the expression patterns of TaBES1 genes in roots, stems, leaves, spikes, anthers, and seed from XBM by qPCR, which revealed that expression levels in root tissues are relatively higher than those in other tissues (Figure S5). TaBES1-1A belongs to the same subgroup as SlBZR1, which also contains Arabidopsis BES1, BZR1, and BEH1 involved in anther development (Figure 1). According to the wheat-expression databases and qPCR results, TaBES1-3A2, TaBES1-3B2, TaBES1-3D2, TaBES1-4B, and TaBES1-4D exhibited relatively higher expression levels in anthers than in other tissues and were therefore considered anther-preferential genes (Figure S5).

To determine how TaBES1 transcription responds to abiotic stresses, we downloaded expression data for TaBES1 genes in plants exposed to drought, salt, cold, or heat from the WheatOmics database; most TaBES1 genes were substantially downregulated under these abiotic stress conditions (Figure S6). To validate the expression patterns of TaBES1 genes, we conducted qPCR using total RNA extracted from fourteen-day-old XBM seedlings subjected to drought, salt, cold, and heat abiotic stresses for 0, 1, 2, 5, 10, 24 h. Finally, we observed a high degree of congruence between qPCR and RNA-seq data, with the transcript levels of most genes being affected by stress treatments visualized as heatmaps using TBtools software with log scale base as 2 and log with as 1 (Figure 3). Drought stress considerably downregulated the expression levels of most TaBES1 genes, with the exception of TaBES1-2A, TaBES1-6D, and TaBES1-7A, which were upregulated under the same conditions (Figure 3A). In response to salt stress, the expression levels of all family members tended to decline, except TaBES1-3A2, which was considerably upregulated (Figure 3B). The expression levels of almost all TaBES1 genes were very low upon exposure to cold stress, while TaBES1-3D1 and TaBES1-6D exhibited an upregulation after low-temperature treatment (Figure 3C). The expression of TaBES1-3A2, TaBES1-6B1, TaBES1-6D, TaBES1-7D, and TaBES1-Un was induced at the beginning of heat treatment, while that of the other TaBES1 genes was substantially suppressed (Figure 3D).

To further elucidate the drought resistance level of each wheat cultivar at the seedling stage, we conducted both water-withhold treatment and PEG-induced drought stress on the six cultivars. The experimental design included a well-watered control group and a drought stress group where water was withheld, validating the previously reported levels of drought resistance. Specifically, TY806, XN318, and ZY6 were found to exhibit drought tolerance, while XBM, ZY9507, and YM8 were identified as drought-sensitive cultivars (Li L. et al., 2018; Cao et al., 2021) (Figure S8A). These findings were further confirmed by the application of PEG treatment, which resulted in significantly higher levels of SOD values in the three drought-resistant cultivars compared to the sensitive cultivars. Similarly, POD values exhibited a similar trend to that of the SOD values. Furthermore, the MDA values in the three drought-resistant cultivars were significantly lower than those observed in the sensitive cultivars (Figure S8B). Based on these significant findings, the aforementioned three cultivars (TY806, XN318, ZY6) were classified as drought-resistant, whereas the remaining three cultivars (YM8, ZY9507, XBM) were categorized as drought-sensitive.

We further investigated the roles of TaBES1 members in response to drought stress by measuring their expression profiles in drought-tolerant and drought-sensitive cultivars. Wheat cultivar XBM is stress-sensitive at the seedling stage (Li L. et al., 2018; Cao et al., 2021); we determined that most TaBES1 genes are considerably suppressed in response to drought stress (Figure 3A). By contrast, the expression patterns of TaBES1-2A, TaBES1-3B2, TaBES1-3D2, TaBES1-6B1, TaBES1-6D, and TaBES1-7D were different in the drought-tolerant cultivar TY806 from those in XBM (Figure S7). We observed similar expression profiles for these six genes in three drought-sensitive cultivars (XBM, ZY9507, and YM8) and three drought-tolerant cultivars (TY806, XN318, and ZY6) at the seedling stage (Li L. et al., 2018) (Figures 4A, B).In general, we detected no obvious commonality in the expression patterns of TaBES1-2A, TaBES1-6B1, or TaBES1-7A between drought-tolerant and -sensitive cultivars (Figure S8C). The relative expression levels of TaBES1-3B2 in drought-tolerant cultivars showed a substantial elevation compared to those measured in drought-sensitive cultivars (Figure 4C). TaBES1-3B2 is the homolog of the positive drought regulator gene TaBZR2 (TaBES1-3D2), which positively regulates the drought stress response(Cui et al., 2019). Thus, we hypothesize that TaBES1-3B2 might have a role in drought signaling, and its function should be investigated in more detail future experiments. The expression patterns of TaBES1-6D were different among drought-sensitive and drought-tolerant cultivars; TaBES1-6D expression was greatly upregulated in drought-sensitive cultivars but downregulated in drought-tolerant cultivars. We therefore proposed that TaBES1-6D is a negative regulator in the response to drought stress.

Previous studies have reported that BES1 members play roles in anther development (Chen et al., 2019a; Chen et al., 2019b). Therefore, besides the responses to different abiotic stresses, we also explored the potential functions of anther-preferential TaBES1 genes. Using data from the wheat expression database and qPCR results, we determined the expression levels of anther-preferential wheat BES1 genes (TaBES1-1A, TaBES1-3A2, TaBES1-3B2, TaBES1-3D2, TaBES1-4B, and TaBES1-4D) in sterile lines and a common wheat cultivar at the S6~S9 anther developmental stages (Browne et al., 2018), which are extremely important for regulation of anther fertility (Xu et al., 2020). In the case of anthers of TGMS lines under sterile conditions (NY), the sensitive period of anther development aligns with a low-temperature environment, leading to significant inhibition in the development of TGMS anthers and a low seed setting rate. Conversely, under fertile conditions (BJ), anther development is delayed compared to the sterile conditions, occurring in a more favorable moderate-temperature environment, resulting in normal seed setting (Li et al., 2006; Bai et al., 2017; Xu et al., 2020) (Figure S9A). Under moderate temperature fertile conditions, the pollen morphology of TGMS lines was similar to the common wheat cultivar, with almost full filled starch; under sterile (low temperature) conditions, however, the pollen of TGMS lines nearly no starch-filled, while the common wheat cultivar pollen was still comparable to that seen under fertile conditions (Figures 5A, S9A).

In general, no substantial changes were detected in expression levels between the TGMS lines and the common wheat cultivar under fertile conditions (Figure 5B); however, we detected substantial differences between the TGMS lines and the common wheat cultivar under sterile (low temperature) conditions (Figure 5B), especially for the expression patterns of TaBES1-3B2 and TaBES1-3D2 (Figure 5B). Indeed, TaBES1-3B2 expression showed no considerable difference in the common wheat cultivar between sterile and fertile conditions merely with a slight decline under sterile condition, whereas in the TGMS lines, TaBES1-3B2 expression was higher under sterile conditions as compared to fertile conditions. Notably, the expression pattern of TaBES1-3D2 was similar to that of TaBES1-3B2, showing a considerable rise in expression levels in the TGMS lines under low-temperature conditions. Then the transcriptome analysis was further performed to investigate the expression patterns of TaBES1s in TGMS lines and common wheat cultivar (Figure S10), and a high degree of congruence was observed between qPCR and RNA-seq results. Generally, in the common cultivar, TaBES1-3B2 and TaBES1-3D2 expression showed no considerable difference under fertile and sterile conditions; however, under sterile conditions, the transcript levels of TaBES1-3B2 and TaBES1-3D2 were substantially higher than the fertile conditions in TGMS lines.