The hidden Markov model (HMM) profile of the RCC1 gene family (PF00415) in PFAM (http://pfam.xfam.org/) was downloaded and used to identify the RCC1 genes in the local protein database of bread wheat, emmer wheat (Triticum dicoccoides, 2n = 4x = 28, AABB) and Aegilops tauschii (2n = 2x = 14, DD) (downloaded from Ensembl Plants, http://plants.ensembl.org/index.html) with the hmmsearch tool of HMMER3.1 software (HMMER 3.1; http://hmmer.org/). To avoid missing RCC1 family members, an aligned file of a high-quality protein set (E value < 1 × 10⁻²⁰) in MEGA X software ( Kumar et al., 2018 ) was used to reconstruct the new HMM profile, which was used as the query to search all the RCC1 members (E value < 0.01) in all bread wheat, emmer wheat, and Aegilops tauschii proteins, respectively. All the detected protein sequences were submitted to the PFAM (http://pfam.xfam.org/), SMART domain search (http://smart.embl.de/smart/batch.pl) ( Letunic et al., 2021 ), and NCBI Batch CD-search database (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) to confirm the structural integrity of the RCC1 domain ( Marchler-Bauer et al., 2005 ). The non-redundant, verified genes encoding proteins with RCC1 domains were assigned as members of the RCC1 gene family.

Multiple alignments of the conserved RCC1 protein sequences of bread wheat, emmer wheat, and Aegilops tauschii were performed using Clustal Omega ( Sievers et al., 2011 ) using default parameters, and a phylogenetic tree was constructed using a maximum-likelihood method with 1,000 bootstrap replications in the RaxML_NG software ( Kozlov et al., 2019 ). Figtree 1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/) was used to visualize and optimize the phylogenetic tree.

The RCC1 gene loci of wheat and its related genome donors were extracted from the corresponding annotated gff3 file (downloaded from Ensembl Plants, http://plants.ensembl.org/index.html) using a perl script. The Multiple Collinearity Scan toolkit (MCScanX) was used to analyze the gene collinearity among wheat, emmer wheat, and Aegilops tauschii with the default parameters ( Kozlov et al., 2019 ).

Homolog analysis of RCC1 genes among the A, B, and D genomes of wheat was performed based on the aligned result. The chromosomal distribution and collinearity of RCC1 genes among the wheat and its donors and of the homoeologous RCC1 genes among A, B, and D genomes were visualized by the circle package in R ( Gu et al., 2014 ).

Clustal Omega ( Sievers et al., 2011 ) was used to analyze RCC1 protein sequences of wheat, and RaxML_NG ( Kozlov et al., 2019 ) was used to construct a phylogenetic tree via a maximum-likelihood method with 1000 bootstrap replications. The domains of the RCC1 gene family in wheat were verified by the SMART domain search (http://smart.embl.de/smart/batch.pl) ( Letunic et al., 2021 ) and the NCBI Batch CD-search database (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) ( Marchler-Bauer et al., 2005 ). The conserved motifs of the RCC1 gene family in wheat were determined by the online Multiple Em for Motif Elicitation (MEME) suite program (http://meme-suite.org) ( Bailey et al., 2006 ). The software TBtools ( Chen et al., 2020 ) was used to visualize the gene structure, protein domains, and motifs of the RCC1 genes according to the annotated GFF files, the genome sequence of wheat, and the protein domain file from the SMART domain search database, as well as the motif result files from the MEME suite.

The promoter sequences (2-kb upstream) of the TaRCC1 genes were extracted from the wheat reference genome (IWGSC RefSeq v1.1, International Wheat Genome Sequencing, C 2018) using the GTF/GFF3 Sequences Extract function of TBtools ( Chen et al., 2020 ), and their potential cis-acting elements were predicted by submitting to PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

To analyze the tissue-specific and stress-induced expression of TaRCC1 genes, the RNA-seq expression data of five publicly available studies ( International Wheat Genome Sequencing, C 2014; Zhang et al., 2014; Li et al., 2015; Liu et al., 2015; Zhang et al., 2016 ) were obtained from expVIP Wheat Expression Browser (http://www.wheat-expression.com/) ( Ramirez-Gonzalez et al., 2018 ) and Triticeae Multi-omics Center (http://202.194.139.32/expression/index.html) and then visualized using the pheatmap package of R software.

A KnetMiner web application for wheat (https://knetminer.com/Triticum_aestivum/, HassaniPak and Keywan, 2017) was used to search gene-evidence networks extracted from the knowledge network and predict the transcription factors for the three TaRCC1 genes. The expression pattern of the predicted transcription factors was detected by qRT-PCR.

The full-length coding DNA sequences (CDS) of RCC1-3A (TraesCS3A02G362800), RCC1-3B (TraesCS3B02G395200), RCC1-3D (TraesCS3D02G356500), Myb-7B (TraesCS7B02G188000), and Myb-7D (TraesCS7D02G295400) were inserted into pCambia1300-35S-GFP, creating RCC1-3A::GFP, RCC1-3B::GFP, RCC1-3D::GFP, Myb-7B::GFP, and Myb-7D::GFP fusion vectors. The recombinant plasmids were mixed with the nuclear marker NLS-mCherry and transfected into wheat mesophyll protoplasts as previously described by Yoo et al. (2007). The transfection mixture was induced by PEG-Ca²⁺, and the protoplasts were cultured for 12 h at 25°C. The protoplasts were observed and photographed with a fluorescence microscope (Zeiss Imager A2, Germany).

The bread wheat cultivar, CS, was grown in a greenhouse with controlled conditions of 26°C/14 h light and 20°C/10 h dark. Three different treatments were applied, namely salt stress, cold, and drought stress induced by polyethylene glycol (PEG). During the two-leaf stage, seedlings were treated with Hoagland liquid medium containing 200 mM NaCl for 1, 3, and 6 h (salt stress), 4°C for 1, 3, and 6 h (cold stress), and 20% PEG4000 for 1, 3, and 6 h (drought stress). Seedlings grown in a normal environment without treatment were set as the control. Three biological replicates were set for all the trials.

The total RNA of plant materials was extracted using the RNAprep Pure Plant Kit (polysaccharide- and polyphenolic-rich) (TIANGEN, Beijing, China), following the manufacturer’s instructions. A NanoDrop One spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and agarose gel electrophoresis were used to assess RNA quantity and purity. Complementary DNA (cDNA) was synthesized using Reverse Transcriptase M-MLV (Takara, Beijing, China) according to the manufacturer’s instructions. qRT-PCR was performed using a 7300 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) according to the supplier’s instructions. A total of 6 L of DNase/RNase-free water, 11 μl of TB Green Real-Time PCR master mix, 2 μl of diluted cDNA product, and 1 μl of gene-specific primer was added to each reaction mixture. Three biological replicates were used for each tissue and three technical repeats for each biological replicate. The thermal cycle was set as follows: denaturing at 95°C for 30 s, then denaturing at 95°C for 15 s, and annealing and elongating at 58°C for 30 s with 45 cycles. The GAPDH gene was used as an internal reference for the normalization of the expression of the TaRCC1 genes. The relative expression levels were calculated using the 2^−ΔΔCt method.

To identify RCC1 genes in bread wheat, emmer wheat, and Aegilops tauschii, a genome-wide search was performed by local BLASTP using HMM profiles. In total, 149 RCC1 genes, comprising 76 TaRCC1s, 49 TdRCC1s, and 24 AetRCC1, were identified and verified by detecting the RCC1-conserved domain via the Pfam (http://pfam.xfam.org/), SMART domain search (http://smart.embl.de/smart/batch.pl), and NCBI Batch CD-search database (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi). The details of the identified RCC1 genes are listed in Supplementary Table S1 . The distribution of the RCC1 genes on chromosomes, different homoeologous groups, and sub-genomes was determined ( Figure 1 ). Emmer wheat (AABB) and Aegilops tauschii (DD) are two genome donors of bread wheat; thus, the RCC1 genes of emmer wheat and Aegilops tauschii were integrated together to compare with bread wheat. The number of RCC1 genes on chromosomes, homoeologous groups, and sub-genomes showed little difference between bread wheat and the combined data for emmer wheat and Aegilops tauschii. Most RCC1 genes were located in homoeologous groups 1, 2, and 3, while no RCC1s were detected in homoeologous group 4. Two RCC1 genes (TraesCSU02G009000LC and TraesCSU02G009100LC) on ChrUn of bread wheat were certificated to belong to Chr1B through phylogenetic and synteny analysis with the RCC1s of emmer wheat and Aegilops tauschii as follows.

To investigate the phylogenetic relationships and compare the evolutionary relationships of RCC1 genes among bread wheat, emmer wheat, and Aegilops tauschii, a maximum-likelihood phylogenetic tree was constructed using the protein sequences of RCC1s ( Figure 2 ). The best-fit model to construct the tree was LG+FC+G8m, and the RCC1s were classified into four subfamilies (sub. I–IV) and named RCC1 I–IV. The RCC1 I, II, III, and IV subfamilies contained 16 (nine for wheat, five for emmer wheat, and two for Aegilops tauschii), 55 (28 for wheat, 18 for emmer wheat and nine for Aegilops tauschii), 24 (12 for wheat, eight for emmer wheat and four for Aegilops tauschii), and 54 RCC1 genes (27, 18, and nine), respectively. Interestingly, in each subgroup, the number of RCC1 genes from bread wheat (AABBDD), emmer wheat (AABB), and Aegilops tauschii (DD) showed approximately a 3:2:1 ratio, which indicated that the RCC1 gene is evolutionarily conserved across the three species.

There are 76 wheat TaRCC1 genes mapped to 18 of the 21 wheat chromosomes, except 4A, 4B, and 4D, according to available annotation information of the wheat genome. Synteny analysis showed that most TdRCC1s (except TRIDC3AG069370.1) and all the AeRCC1s were highly collinear with the TaRCC1s, and phylogenetic analysis indicated the collinear RCC1 genes of the three species were clustered together ( Figures 2 , 3 ; Supplementary Table S2 ). Homologous gene analysis indicated that, except for TraesCS5D02G123500.1, most TaRCC1 genes were homoeologous to each other among the A, B, and D genomes, and clustered together ( Figures 2 , 3 ; Supplementary Table S3 ). It was noteworthy that two unanchored genes, TraesCSU02G009100LC.1 and TraesCSU02G009000LC.1 on ChrUn, should be anchored on chromosome 1B for their high collinearity to TraesCS1D02G073500.1 and TraesCS1A02G071100.1; therefore, we adjusted the positions of TraesCSU02G009100LC.1 and TraesCSU02G009000LC.1 on chromosome 1B for the visualization of collinearity analysis ( Figure 3 ). According to the descriptions, a chromosomal region within 200 kb containing two or more genes is defined as a tandem duplication event. The gene pairs TraesCS1A02G071100.1/TraesCS1A02G071200.1 and TraesCSU02G009000LC.1/TraesCSU02G009100LC.1 were each clustered into one tandem duplication event region on chromosomes 1A and 1B of bread wheat, respectively; moreover, no homoeologous gene on chromosome 1A was found for the two homoeologous genes TraesCS3B02G543500.1 and TraesCS3D02G489600.1. In general, the TaRCC1s, TdRCC1s, and AeRCC1s on corresponding chromosomes show high collinear with each other according to the synteny analysis ( Figure 3 ), and the number of RCC1 genes on corresponding chromosomes, homoeologous groups, and sub-genomes in bread wheat, emmer wheat, and Aegilops tauschii showed approximately a 1:1 ratio, respectively ( Figure 1 ), which indicated that RCC1 genes were highly conserved during the evolution of wheat, and the expansion of RCC1 gene family was mostly due to the genome polyploidization.

To further estimate the gene structure, protein-conserved domains, and motifs of wheat TaRCC1 genes, the full-length protein sequences of 76 TaRCC1s were aligned using Clustal Omega, and the phylogenetic tree was constructed using RaxML_NG ( Figure 4A ). The TaRCC1s in wheat were classified into four subfamilies named TaRCC1 I–IV. The TaRCC1 I, II, and III subfamilies contained 10, 27, and 12 genes, respectively, and carried only RCC1 domain repeats, while the remaining 27 TaRCC1s of the TaRCC1 IV subfamily contained multiple domains, including RCC1 domain repeats and PH or BRX domains ( Figure 4C ). The conserved motifs of the TaRCC1 genes were determined by the online MEME suite program (http://meme-suite.org): 20 conserved motifs with lengths from 11 to 41 amino acids were detected among the TaRCC1 genes ( Supplementary Table S4 ; Figure 4B ). TaRCC1s in the same cluster shared similar conserved motif compositions, which again indicated that there is high conservation of the RCC1 gene family sequence in wheat. Despite the similarity of motifs among closely related genes, the size of the gene fragments varied widely (390–23,462 bp), such that the TraesCS2D02G725400LC gene fragment was much smaller than TraesCS7B02G200200. The gene structure, including the size and number of intron–exon, varies a lot among different TaRCC1s (such as the number of exons is from 1 to 17) ( Figure 4D ). It is worth noting that the closely related members, especially homologous genes, showed similar exon–intron structure, and the difference among them was the exon–intron length. The homologous genes TraesCS7A02G284800, TraesCS7B02G200200, and TraesCS7D02G282600 were similar in motif, protein domains, and exon–intron structure, while their exon–intron lengths varied greatly (10,997 bp, 23,461bp and 14,257 bp, respectively) ( Figure 4 ; Supplementary Table S1 ).

Cis-acting elements in gene promoters are crucial regions for initiating transcription at transcription factor-binding sites, which play an important role in regulating gene expression. The potential cis-acting elements on the promoter regions (2 kb upstream) of TaRCC1s were analyzed by PlantCARE to further explore their possible biological functions (details in Supplementary Table S5 ). Various potential cis-acting regulatory elements in the promoter regions of TaRCC1 genes were predicted to be related to transcription, cell cycle, development, hormones, and response to stresses ( Figure 5 ; Supplementary Table S5 ). All of the TaRCC1 genes contained light-responsive elements. A total of 70 and 74TaRCC1s were detected with MeJA-responsive elements and ABA-responsive elements (ABRE), respectively. In addition, many elements were predicted to be involved in various abiotic stresses, such as drought, salt, cold, and light ( Figure 5 ; Supplementary Figure S1 ; Supplementary Table S5 ).

Using the available RNA-seq database of ( International Wheat Genome Sequencing, C 2014 ) obtained from the expVIP Wheat Expression Browser (http://www.wheat-expression.com/), the temporal and spatial expression patterns of 76 TaRCC1 genes in five different tissues (root, stem, leaf, spike, and grain) ( Supplementary Table S6 ) were visualized using the heatmap package of R software ( Supplementary Figure S2 ). The expression levels of TaRCC1s varied significantly among different tissues. Some TaRCC1s from the same group showed similar expression patterns, while others indicated diverse expression patterns. For example, TraesCS2B02G323100, TraesCS2A02G306400, and TraesCS2D02G304900 from TaRCC1 II were predominantly expressed in the root and stem, while the TraesCS3A02G362800, TraesCS3B02G395200, and TraesCS3D02G356500 from TaRCC1 II were most strongly expressed in the leaf, followed by the spike, root, and early stage of the stem. Similar expression patterns were observed for most homoeologous genes, although others presented diverse patterns. For instance, TraesCS6A02G153100 and TraesCS6D02G142500 were predominantly expressed in the root and stem, while the homoeologous gene TraesCS6B02G180900 presented very low expression in the five tissues. It was worth noting that, in the TaRCC1 I family, a total of four TaRCC1 genes exhibited no expression in the five tissues.

The available RNA-seq data from four studies ( Zhang et al., 2014; Li et al., 2015; Liu et al., 2015; Zhang et al., 2016 ), obtained from the expVIP Wheat Expression Browser (http://www.wheat-expression.com/) and Triticeae Multi-omics Center (http://202.194.139.32/expression/index.html), were used to study the expression of wheat RCC1s in response to salt, drought, heat, cold, and stripe rust stresses. The transcript-per-million-read (TPM) values of TaRCC1 genes are presented in Supplementary Table S7 ; values were transformed by log₂(x+1) and used for visualization with the pheatmap package of R software. The expression patterns of TaRCC1 genes varied a lot under different stresses ( Figure 6A ). Although some homoeologous genes presented diverse expression patterns, most of them exhibited similar expression patterns. The homoeologous genes TraesCS2A02G306400, TraesCS2B02G323100, and TraesCS2D02G304900 showed similar high expression trends under the four stresses and no significant differential expression under the different stresses, except for salt stress, indicating that these three genes might be induced by salt stress. The gene TraesCS2A02G456700 exhibited lower expression only under salt stress, while the homoeologous genes TraesCS2B02G478900 and TraesCS2D02G457000 presented differential expression, especially under cold stress, suggesting these two genes might participate in the cold tolerance of wheat. The homoeologous genes TraesCS3A02G362800, TraesCS3B02G395200, and TraesCS3D02G356500 showed similar and significant differential expression under the five different stresses. Under salt treatments at 6, 12, 24, and 48 h, the expression level of the three genes were higher than in the control; with drought and heat treatments, the three genes showed significantly higher expression under drought and/or heat treatments for 1 h, then decreasing after treatment for 6 h. Similar trends were detected under cold (4°C) and stripe rust pathogen stresses ( Figure 6B ). The results indicated that these three genes might be the key genes that participated in the early stress responses of wheat under stress, and might alleviate the stress injury of plants.

The gene-evidence networks extracted from the knowledge network of wheat through KnetMiner showed several transcription factors for the three TaRCC1 genes ( Figure 7A ). Five genes were identified as candidates for participating in the regulation of the three TaRCC1s, among which Myb-7B (TraesCS7B02G188000) and Myb-7D (TraesCS7D02G295400) were associated with the regulation of all three TaRCC1s (RCC1-3A, RCC1-3B, and RCC1-3D), suggesting that Myb-7B and Myb-7D might be the regulators of the three TaRCC1s. Myb-7B and Myb-7D encode two Myb-like transcription factors, which were related to the terms stripe rust response and drought tolerance in the KnetMiner knowledge network ( Figure 7A ). The expression patterns of Myb-7B and Myb-7D under multiple stresses, obtained from the expVIP Wheat Expression Browser (http://www.wheat-expression.com/) and the Triticeae Multi-omics Center (http://202.194.139.32/expression/index.html), showed similar and significant differential expression under different treatments of the five stresses. For instance, the expression level of the two Myb genes was higher than in the control for salt treatments at 6 h and 12 h and lower than in the control at 24 and 48 h. Under cold (4°C) and stripe rust pathogen stresses, the genes showed a similar gene expression pattern. Under drought and heat treatments, the Myb-7B/D showed reduced expression under dr_6h, heat_6h, and dr_heat_6h treatments, compared with the control ( Figure 7B ). Myb-binding sites on the promoters of TaRCC1-3A, TaRCC1-3B, and TaRCC1-3D were predicted using the PlantRegMap software (http://plantregmap.gao-lab.org/binding_site_prediction.php). A potential Myb-binding site was found in the promoters of TaRCC1-3A (−2119 to −2133) and TaRCC1-3B (−227 to −241), but no predicted Myb-binding sites were detected immediately upstream of the transcription start site of TaRCC1-3D, indicating that TaRCC1-3A/B are most likely directly regulated by the above Myb-like transcription factor, while TaRCC1-3D might be indirectly regulated ( Figure 7C ).

In Arabidopsis, RCC1 family proteins, such as UVR8 and TCF1, are located in the nucleus. Similarly, subcellular localization prediction suggested that all of the TaRCC1 family proteins are located in the nucleus. Our experiments to investigate the subcellular localization of the three TaRCC1 proteins (TaRCC1-3A (TraesCS3A02G362800), TaRCC1-3B (TraesCS3B02G395200) and TaRCC1-3D (TraesCS3D02G356500) in wheat protoplasts confirmed the results as predicted and showed that these proteins are located in the nucleus ( Figure 8 ). The subcellular localization of the Myb-7B/D proteins showed that these two Myb proteins are also located in the nucleus ( Figure 8 ). The primers used are listed in Supplementary Table S8 .

From the available RNA-seq data of several studies, the expression patterns of TaRCC1s in different tissues and multiple stresses had been analyzed, as described above ( Supplementary Figure S2 ; Figure 6 ). The expression levels of TaRCC1s varied significantly in different tissues and under multiple stresses; several TaRCC1s were induced by different stresses. Homoeologous genes TraesCS3A02G362800, TraesCS3B02G395200, and TraesCS3D02G356500 responded to five stresses (salt, drought, heat, cold, and stripe rust pathogen). We used qRT-PCR to verify the expression patterns of the three TaRCC1s in response to salt, cold, and drought stresses ( Figure 9 ). Overall, these three TaRCC1s were induced by almost all the treatments, showing similar patterns to the RNA-seq results mentioned above. The expression of Myb-7B and Myb-7D was analyzed by qRT-PCR, and they shared a similar gene expression pattern to RCC1-3A, RCC1-3B, and RCC1-3D under salt, drought, and cold treatment ( Figure 9 ), further suggesting that Myb-7B and Myb-7D might be regulators of the three TaRCC1s above. The qRT-PCR primers for RCC1-3A, RCC1-3B, RCC1-3D, Myb-7B, and Myb-7D are listed in Supplementary Table S9 .