Identification and characterization of the CRK gene family in the wheat genome and analysis of their expression profile in response to high temperature-induced male sterility

Wheat material, paraffin sections and characterization of phenotypes

The wheat seeds of cv. Zhoumai 36 were sown on October 28, 2022, in the experimental field of Zhoukou Normal University, Zhoukou, Henan Province, People’s Republic of China (33°C64′N, 114°C6′E). The wheat material was covered with plastic film in early April of the following year and all treatment procedures were carried out as described in our previously research (Liu et al., 2018a; Liu et al., 2022). Samples from normal anthers and HT-ms anthers were collected in six centrifuge tubes (1.5 mL), quick-frozen in liquid nitrogen and quickly stored at −80 °C. Three of the samples were used for qRT-PCR experiments within half a month and the remaining three were kept in reserve. In addition, anthers were treated with FAA fixative and placed in a 4 °C refrigerator and set aside for paraffin sectioning. We set the longitudinal section thickness to 12 µm and detected insoluble carbohydrates (especially polysaccharides and starch granules) with a periodic acid-Schiff (PAS) staining technique. The fertility or sterility of the samples was preliminarily based on the staining results. Images were collected with a light microscope. All experiments were completed by 20 October 2023.

Identification and analysis of wheat CRK family members

The predicted protein sequences were downloaded from Ensembl Plants to identify CRK genes in wheat (IWGSC RefSeq v1.1: http://plants.ensembl.org/index.html) (International Wheat Genome Sequencing Consortium , 2018). The hidden Markov model (HMM) profiles (PF01657; PF00069; PF07714) corresponding to the CRK gene family were downloaded from the Pfam database (Pfam 35.0: http://pfam.xfam.org/). The FASTA and GFF3 files containing candidate wheat sequences for the CRK domains were downloaded from the wheat genome (http://plants.ensembl.org). Subsequently, we rebuilt a new wheat-specific HMM file for each HMM file using the HMMER v3.0 suite hmmbuild (Eddy, 2009) and used it to identify all the wheat proteins. The identification results of the three hmm files were intersected to obtain a total of 132 CRK proteins with a standard E-value <1 × 10⁻²⁵. Then, Expasy Protparma (https://www.expasy.org/) was used to analyze the physicochemical properties of the 132 TaCRKs proteins, including the Compute pI/MW, CDS length, and number of amino acids. The number of transmembrane domains was predicted by the TMHMM server (https://services.healthtech.dtu.dk/services/TMHMM-2.0/) (Möller, Croning & Apweiler, 2001). The subcellular localization was predicted by the BUSCA web server (http://busca.biocomp.unibo.it/) (Savojardo et al., 2018). A total of 132 candidate CRKs were finally obtained from the wheat genome database using the wheat-specific CRK HMM profiles and further validated by the PFAM database (Paysan-Lafosse et al., 2023) and the SMART database (Letunic, Khedkar & Bork, 2021a), and these proteins were further analyzed.

Analyses of phylogeny, gene structure, conserved motifs and protein domain

The CRK protein sequences of tomato (Solanum lycopersicum), soybean (Glycine max), rice (Oryza sativa) and A. thaliana were retrieved and downloaded from the Ensembl Plants database (http://plants.ensembl.org/index.html) and TAIR database (https://www.arabidopsis.org). MEGA-X software (Kumar et al., 2018) and the neighbor-joining method (NJ) were used to construct the phylogenetic tree for amino acid sequences with the bootstrap repeated value set to 1,000 times, followed by exporting the newick file and embellishing the phylogenetic tree using iTOL (https://itol.embl.de/login.cgi) (Letunic & Bork, 2021b). The gene structure and conserved motif analysis were performed according to previously published methods (Liu et al., 2021c). The domain distribution was analyzed using the NCBI-CDD website (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). Finally, these files were visualized and optimized using TBtools (Chen et al., 2020).

Analysis of the CRK family promoters, chromosomal location, gene duplication, and synteny in wheat

Gene expression is associated with regulatory elements of promoter sequences. We extracted and obtained the upstream sequences of 95 CRK genes CDS of 2-kb length from the genome sequence and used them as identification sites for cis-regulatory elements in the promoter region. These sequences were then submitted to the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002) for cis-element prediction and visualized with TBtools. Two-way comparison of genomic sequences was performed by the blast tool in the TBtools software. MCScanX Wrapper (Wang et al., 2012) was used to identify TaCRKs duplications and the syntenic blocks contacting CRK genes of wheat, A. thaliana, foxtail millet, rice, and soybean. Subsequently, the gene duplication results and the collinearity comparison map were visualized with TBtools (Chen et al., 2020).

Prediction of putative miRNAs targeting TaCRK genes and enrichment analysis

All TaCRK CDSs were uploaded to the psRNATarget website (Dai, Zhuang & Zhao, 2018) with parameters set to default, which in turn predicted miRNA targets. Subsequently, the miRNA target network map of TaCRK was generated with the assistance of Cytoscape-v3.9.1 (Shannon et al., 2003). Gene Ontology enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment of the TaCRK genes were analyzed using the EGGNOG-Mapper and TBtools.

Expression profiles of TaCRK genes in different wheat tissues

The wheat RNA-Seq data (IWGSC Annotation v1.1) were downloaded from WheatOmics 1.0 (http://wheatomics.sdau.edu.cn/) to explore the expression profiles of TaCRK genes in five different tissues. TBtools software was used to draw the expression heatmap of wheat CRK genes.

Total RNA extraction and quantitative real-time PCR analysis

Total RNA was extracted from wheat anther samples at the mononuclear and trinuclear developmental stages using TRNzol Universal (Tiangen Biotech, Beijing, China) according to the manufacture’s requirements. The RevertAid first strand cDNA synthesis kit & DNase I (Thermo Fisher Scientific, Waltham, MA, USA) were used for reverse transcription according to the manufacturer’s instructions. The reverse-transcribed cDNA was used as a template, the wheat actin gene was used as an internal control, and each gene-specific primer was used as the primer. Fluorescence quantification was performed by BIO-RAD CFX Connect™ Fluorescent PCR Detection System according to the Power SYBR^® Green PCR Master Mix reaction system. The primer sequence designed with TBtools and Primer Premier 5.0 software are detailed in Table S1. The total reaction volume was 25-µL, which contained 1 µL diluted cDNA and 0.5 µL (10 mM) each primer. The amplification program consisted of an initial step at 95.0 °C for a duration of 3 min. Subsequently, the temperature was set to 95 °C for 10 s, followed by 55 °C for 20 s, 72 °C for 20 s, and 75 °C for 5 s. This cycle was repeated a total of 40 times. The wheat actin gene (AB181991.1) was used as a housekeeping gene for internal control. Three biological replicates of each treatment were performed, and the relative expression of genes was calculated using the 2^−ΔΔCT method. SPSS version 27.0 (IBM Corp, 2020, Armonk, NY, USA) was used to analyze the significance of differences. P < 0.05 indicated significant differences, while P < 0.01 indicated highly significant differences.

Identification of the CRK gene family in wheat

A total of 132 CRK proteins were identified from the wheat genome using Blast comparison and a HMMER 3.0 search with three PFAM files. The intersection was determined and (Fig. 1) analyzed by SMART, removing sequences without CRK structural domains, and combining different transcripts of the same gene (the longest transcript sequence represents the gene) to obtain a total of 95 genes, as shown in Table 1. The genes were named TaCRK1-1A to TaCRK95-U in sequence according to the order of the gene’s position on the chromosome. Analysis of the physicochemical properties of the protein revealed that the molecular weight of wheat CRK protein ranged from 54.35 kD (TaCRK72-3D) to 82.75 kD (TaCRK66-3B), and the theoretical isoelectric points (pIs) were distributed in the range of 5.21 (TaCRK44-2D) to 8.79 (TaCRK33-2B). The predicted transmembrane domain (TMDs) analysis of these 132 proteins resulted in 74 proteins containing one TMDs, 51 containing two TMDs, four proteins containing three TMDs, and others without TMDs. There were 117 hydrophilic coefficients less than zero as hydrophilic proteins, and 15 greater than zero as hydrophobic proteins. Furthermore, the predicted subcellular localization indicates that the vast majority of TaCRKs were localized at the plasma membrane, with only a few located in the extracellular space and endomembrane system (Table 1).

Phylogenetic tree construction, gene structure, motifs, and conserved domains of the CRK gene family in wheat

In order to reveal the phylogenetic relationships between wheat and other plant CRKs, an unrooted phylogenic tree of the CRK protein sequences of wheat, A. thaliana, rice, soybean, and tomato were constructed. A total of 330 total CRK protein sequences were derived from A. thaliana (44), wheat (132), rice (40), soybean (91) and tomato (23). According to the phylogenetic tree branching, the 132 proteins corresponding to the 95 TaCRK genes were classified into three major groups, Group 1, Group 2 and Group 3. Group 3 contained two subgroups, Group 3-1 and Group 3-2. A total of 14 wheat CRK proteins were present in Group 1 and only two in Group 2; Group 3-1 contained nine CRK proteins and Group 3-2 had 107 wheat CRK proteins. Group 3-2 contained the most wheat CRK and had 29 rice CRK proteins. All of the groups contained rice CRK proteins; for example, Group 1 had four rice proteins, Group 2 had three rice proteins, and subgroup Group3-1 had one rice protein. All of these rice proteins were present on the same branch as the wheat CRK proteins, indicating that the wheat CRK proteins were more closely related to the rice CRK proteins than the other species (Fig. 2).

Analysis of the gene structure, motifs, and conserved domains of the TaCRK gene family

The amino acid sequences of TaCRK family members were analyzed using the MEME online website to determine their functional regions, and 12 motifs were found. The evolutionary tree of 95 TaCRK family members is divided into three large subgroups, which are indicated by different background colors. Motif 1, 2, 5, 6, 9, 10 are motifs that are present in all protein members. Except for the shared motif, motifs 8 and 11 are found in the first major branch.

The second major branch has five TaCRK members, of which TACRK10-1D lacks motif 3. In addition, all four proteins lack this motif except for TaCRK5-1B, which contains motif 11. The third major branch has 15 TaCRK members, four of which (TaCRK42-2B, TaCRK55-3A, TaCRK63-3B, TaCRK91-7B) are lacking motif 1 (Fig. 3). The exons and introns of the TaCRK gene family were analyzed to further understand the gene structure of this family. There were two to seven exons found in the TaCRK gene (Fig. 3). Among these, three genes (TaCRK36-2B, TaCRK4-1A, TaCRK69-3B) contained seven exons, 13 genes contained five exons, 69 genes contained four exons, eight genes contained three exons, and two genes (TaCRK81-5B, TaCRK82-5B) contained two exons. There were one to nine TaCRK genes containing introns. Among them, 64 genes contained six introns and 13 genes contained seven introns. The genes containing nine introns were TaCRK69-3B and TaCRK61-3A, while those containing only one intron were TaCRK81-5B and TaCRK82-5B. Most of the genes had upstream and downstream UTR regions, but some genes had missing UTR regions, such as TaCRK70-3D and TaCRK37-2B. The predictions of the Conserved Domain Database (CCD) showed that most TaCRK genes contained the STKc-IRAK conserved domain and the stress-autifung conserved domain.

Distribution of wheat CRK genes on chromosomes and gene duplication

Wheat CRK genes are unevenly distributed on 18 chromosomes, and no CRK genes were found on chromosomes 4A, 4B, and 4D. Chromosome 2B contained 15 genes, followed by chromosomes 2A and 2D with 12 genes; chromosomes 6A, 6B, 7A, and 7B were the least numerous, containing only one gene each (Fig. S1). The distribution of TaCRK genes in the same subfamily or evolutionary branches on chromosomes also showed a tendency to be located within the same subgenome. For example, the gene clusters TaCRK56-3A to TaCRK62-3A are located on the same branch in the phylogenetic tree (Fig. S1; Fig. 3. The Ka/Ks ratios of 109 pairs of TaCRK paralogous homologs were calculated to further investigate the effect of evolutionary factors on the TaCRK family. The results showed that the Ka/Ks ratios of all TaCRK genes did not exceed 1, with the largest being 0.67 (TaCRK7-1B:TaCRK9-1B) and less than 0.5 for 99 pairs of genes. This Ka/Ks ratio results suggest that these paralogous homologous genes with ratios less than 1 may have undergone strong purifying selection during evolution (Fig. S2; Table S2). Segmental replication and tandem replication are thought to be the two major reasons for gene family expansion in plants. Seventy-two percent (68/95) of the wheat TaCRK members showed duplication events, with 15 tandem duplication events out of 68 duplication events (Fig. 4). For the rest, highly similar genes were found in different chromosomes with segmental duplication events. As shown in Fig. 4, duplication events occurred mainly on chromosomes 2A, 2B, and 2D, while no duplication events were observed on 4A, 4B, and 4D (Fig. 4). In addition, our cluster analysis of homologous genes revealed that almost all homologous genes were in the same evolutionary branch. For example, the segmental repeats TaCRK1-1A, TaCRK5-1B, and TaCRK12-1D have homologous loci in all three partial homologous chromosome groups (A, B, D) and are in the same branch of the evolutionary tree (Fig. 3).

Collinearity analysis of CRK genes between wheat and other representative plants

To investigate the origin and evolution of the TaCRK gene in wheat, collinearity analysis of the CRK gene family in wheat and A. thaliana, foxtail millet, rice, and B. distachyon was performed, and the results showed that wheat and A. thaliana, foxtail millet, B. distachyon, and rice have three pairs of orthologous genes, 12 pairs of orthologous genes, 15 pairs of orthologous genes, and 11 pairs of orthologous genes, respectively (Fig. 5). In detail, three pairs of orthologous genes (TaCRK36-2B and AT4G23180; TaCRK38-2B and AT4G23270; TaCRK44-2D and AT4G23230) were found in A. thaliana and wheat. The others are detailed in Table S3 and Fig. 5. The Ka/Ks values of these TaCRK genes were all less than 0.67, indicating that these CRK genes in wheat species have undergone strong purifying selection during evolution and are more functionally conserved. That is, the rate of nonsynonymous substitutions is less than the rate of synonymous substitutions when most nonsynonymous substitutions are deleterious and most nonsynonymous substitutions are purified (Fig. S2).

Analysis of cis-acting elements of the wheat CRK gene

To further understand the regulatory functions of the genes, the cis-elements in the promoter sequences were analyzed. A 2,000 bp promoter sequence upstream of the CRKs initiation codon was extracted from the wheat genome, and cis-element analysis was performed using the online tool PlantCARE. The TaCRKs promoter region was enriched with cis-acting elements that respond to phytohormone and stress in response to adversity. There are five types of cis-elements related to phytohormone response, namely the abscisic acid response element (ABRE), methyl jasmonate response element (CGTCA-motif; TGACG-motif), salicylic acid response element (TCA), gibberellin response element (GARE-motif, P-box and TCTC-box), and auxin hormone response element (TGA-element). The cis-elements associated with growth and development are the regulatory elements of zein protein metabolism. The cis-elements associated with stress response mainly include anaerobic induction response element (ARE), hypoxia response element (GC-motif), low temperature response element (LTR) and defense and stress response element (TC-rich repeats, MBS). This implies a potential role for the wheat CRK gene family in wheat growth and development and in a variety of hormones and stresses (Fig. 6; Table S4).

GO enrichment and KEGG enrichment analysis of TaCRKs

For the further exploration of the functions of TaCRKs proteins in wheat, we performed GO enrichment and analyzed the cellular composition, molecular functions, and biological pathway categories of these proteins. In the cellular composition, the main enrichment entries were “plasma membrane” and “membrane.” In the molecular function category, TaCRKs proteins were highly enriched in “kinase activity” followed by “transferase activity” and “catalitic activity”. Analysis of the biological pathways revealed that most TaCRKs proteins were mainly assigned to the “cellular protein modification process”, “protein metabolic process”, and “metabolic process” (Fig. 7A). Subsequently, KEGG enrichment analysis was performed for these TaCRKs proteins.“ Protein kinases” and “Protein families: metabolism” were the most enriched KEGG enrichment in TaCRK proteins (Fig. 7B).

The microRNA targets of TaCRK genes throughout the genome were identified

MicroRNAs (miRNAs) are a class of small non-coding RNAs that can regulate gene expression at the post-transcriptional level by inhibiting mRNA translation or promoting mRNA degradation. We determined 50 miRNAs targeting 78 genes to better understand how TaCRK genes are altered by miRNAs. As shown in the Sankey diagram, one miRNA was (tae-miR9657b-5p) targeted to 15 TaCRK genes, which was the largest number, followed by tae-miR164 and tae-miR9780. Similarly, a gene can be regulated by multiple miRNAs, such as TaCRK68-5B, TaCRK94-U, and TaCRK52-2D which target six miRNAs, five miRNAs, and five miRNAs, respectively (Fig. 8). To improve the visual representation, we mapped the network interactions of genes and miRNAs, and the results showed that five miRNAs were in the hub of the interactions map. These genes include tae-miR9657b-5p, tae-miR9780, tae-miR9676-5p, tae-miR164, and tae-miR531 (Fig. S3A) and further mapped the miRNA targeting sites of TaCRK85-5D and TaCRK52-2D (Figs. S3B, S3C). Table S5 provides all miRNA targeting sites/genes.

Expression patterns of wheat TaCRK genes in different tissues

To deeply investigate the expression pattern of TaCRK in different tissues, RNA-seq data of Chinese spring wheat varieties were obtained from WheatOmics 1.0. As shown in Fig. 9, TaCRK genes were more highly expressed in roots and stem compared with other tissues, with the highest expressed gene also occurring in roots (TaCRK54-2D). There are some genes that show high expression in various tissues including root, stem, leaf, spike, and grain, such as TaCRK36-2B and TaCRK49-2D (Table S6). These two genes are in the same branch of the evolutionary tree (Fig. 3) and are paralogous homologs with a Ka/Ks value of 0.51 (Fig. 4 and Fig. S2). Several other genes were highly expressed in roots, stems, spikelets, and grains, but showed low or no expression in leaves, including TaCRK16-2A, TaCRK43-2D and TaCRK44-2D. Figure 9 shows that there is a significant number of genes that exhibit high expression levels in one specific tissue, while displaying lower expression levels in all other tissues or exhibiting low expression levels across all tissues. Notable examples of such genes include TaCRK61-3A, TaCRK29-2B, TaCRK79-5A, and TaCRK11-1D.

The relationship between the wheat CRK family and anther sterility induced by high temperature was investigated using qRT-PCR

High-temperature stress has a negative impact on the growth and reproductive process of wheat plants. During development, HT disrupt the development of stamens, resulting in male sterility. HT-ms plants experience restricted growth, characterized by lower plant height (Figs. 10A, 10B). Additionally, HT stress can also affect the morphology and function of wheat anthers. Normal wheat anthers are slightly larger, and exhibit pronounced dehiscence, facilitating the release and dispersal of pollen. However, under HT conditions, anthers of male-sterile wheat may become smaller and less prone to dehiscence, thereby affecting the normal release of pollen (Figs. 10C, 10D). Longitudinal paraffin sections at the monocarpic stage did not exhibit significant differences. However, normal microspores were observed to be slightly larger than sterile microspores, and the epidermis of sterile anthers appeared thicker compared to that of normal anthers (Figs. 10E, 10F). For the trinucleate stage of anthers, normally-developing microspores exhibit abundant starch content in mature pollen grains, while HT-ms anthers show insufficient starch content in microspores, and some even lack starch accumulation. Additionally, there are also structural differences observed in the epidermis between HT-ms anthers and normal anthers (Figs. 10G, 10H).

To further validate the association of these trait changes with CRK genes at the gene level, we selected six genes that have collinearity with other species and analyzed their expression changes during the high-temperature male sterility process in anthers. For TaCRK17-2A, the expression level in the mononuclear stage of sterile anthers was 2.62-FC (fold change) higher than that in Normal anthers, which was significantly higher than that of Normal anthers; whereas the expression in HT-ms anthers at the trinuclear stage was also slightly higher than that of normal anthers but did not reach a significant level (Fig. 11A). The electronic fluorescent pictograph (eFP) demonstrated that this gene was relatively highly expressed in roots (0.92 TPM) and grains (1.42 TPM), while it was relatively low in the spikelets (0.19 TPM) (Fig. 11B). As shown in Figs. 11C–11J, the genes TaCRK38-2B, TaCRK44-2D, TaCRK77-5A, and TaCRK85-5D exhibited a significant increase in expression levels at the mononucleate stage of HT-ms anthers, compared to normal anthers, with fold changes of 2.26, 3.40, 4.18, and 3.97, respectively. However, at the trinucleate stage, they showed a decreasing trend but did not reach a significant level. The TaCRK38-2B gene exhibits high expression levels in the roots (1.55 TPM) according to the eFP, while its expression in other tissues is generally normalized. TaCRK44-2D, TaCRK77-5A, and TaCRK85-5D all show relatively high expression in grains (3.08 TPM, 0.01 TPM and 0.10TPM, respectively). Additionally, the TaCRK77-5A gene also exhibits high expression levels in spikelets and leaves, while the TaCRK85-5D gene shows a highly expressed state in roots as well. For the gene TaCRK93-7D, the expression level at the mononuclear stage of HT-ms anthers was 7.67-FC higher than normal anthers, showing a highly significant difference. Similarly, at the trinuclear stage, the expression level of HT-ms anthers was slightly higher than that of normal anthers, but the difference did not reach a significant level (Fig. 11K). The eFP results indicate that the gene TaCRK93-7D shows relatively higher expression levels in leaves and stems, while it exhibits relatively lower expression in other tissues (Fig. 11L). Detailed data from the qRT-PCR experiments are shown in Table S7. Overall, all genes at the mononucleate stage of sterile anthers exhibit a highly significant differential high expression compared to Normal anthers. Additionally, compared to normal anthers, there are four genes that show a downregulation trend in expression during the trinucleate stage, while two genes show an upregulation in expression, but the difference in expression levels is not significant compared to normal anthers.