Genome-Wide Mining of Wheat DUF966 Gene Family Provides New Insights Into Salt Stress Responses

Plant Material and Stress Treatment

Common wheat (Triticum aestivum L.) cultivar Zhengmai was used in this study. Plants were grown at 26°C under a 16 h day/8 h dark photoperiod. To conduct salt stress treatments, plants were treated with 150 mM sodium chloride (NaCl) at the stages of seedling with one and half leaves. Leaf and root samples were collected at 2, 4, 8, 12, 24, and 96 h after the salt treatment (Jiang et al., 2019). The experiment was performed with three independent biological replicates.

Identification of TaDUF966 Genes and Sequence Analysis

The TaDUF966 gene family members were identified by in silico analysis of the wheat reference genome assembly (RefSeq v.1.1; International Wheat Genome Sequencing Consortium [IWGSC]) (https://wheat-urgi.versailles.inra.fr/Seq-Repository/Assemblies) (He et al., 2020). To obtain the amino acid sequences of the putative TaDUF966 family proteins, BLASTp (Basic Local Alignment Search Tool protein) analysis was conducted using seven rice DUF966s (OsDSRs) (http://rice.plantbiology.msu.edu/index.shtml) (Luo and Tian, 2017) and five Arabidopsis (https://www.arabidopsis.org/index.jsp) DUF966s (AtUOFs) as query sequences against the wheat reference IWGSC v.1.1 (Jiang et al., 2020a). To ensure the reliability of TaDUF966 protein sequences, a cut-off value (E < 1 × 10⁻¹⁰) was used. The TaDUF966 family members were further screened and validated using the SMART (Letunic et al., 2004) and Pfam databases (version 32.0, PF06136) (http://pfam.xfam.org/).

Phylogenetic Analysis, Chromosomal Locations, and Duplication Patterns

Phylogenetic analysis of TaDUF966 proteins was conducted using the LG model (Le and Gascuel, 2008) based maximum likelihood (ML) method in MEGA 7.0 (Kumar et al., 2016). The test guidance was repeated 1,000 times (Felsenstein, 1985) to calculate the support for each node. A midpoint-rooted base tree was drawn using the Interactive Tree of Life (IToL, version 3.2.317, http://itol.embl.de) (Zhu et al., 2019b). The GFF3 gene annotation file was obtained from the wheat database IWGSC v.1.1, and TaDUF966 gene annotations were extracted. Chromosomal start and stop sites and the chromosomal location of each TaDUF966 gene were used to draw the physical map with the MapInspect software. All possible homolog of the DUF966 genes in each subgenome of common wheat were determined using “all-against-all” BLAST, with an E-value cut-off of 1 × 10⁻¹⁰ and identity > 75% (Letunic et al., 2004). The R package “circlize” was used to draw a graph showing its position and homology relationship (Fang et al., 2020). The non-synonymous (Ka) and synonymous (Ks) substitution rates were calculated using DnaSP 6.0 (Rozas et al., 2017) to diagnose the form of sequence evolution.

Gene Structure and Predicted Protein Feature Analyses

The predicted TaDUF966 protein sequences were analyzed using protein identification and analysis tools available online. The protein length, isoelectric point (pI), molecular weight (MW), instability, and index were predicted using the ProtParam (Gasteiger et al., 2005) tool (https://web.expasy.org/protparam/). The subcellular localization of these proteins was predicted using the Plant-mPLoc online tool (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/) (Chou and Shen, 2010).

The MEME (v5.0.3) (Bailey and Elkan, 1994) suite analysis tool and MAST (http://meme-suite.org/tools) motif search tool were used to identify conserved motifs in TaDUF966 proteins. The known DUF966 protein sequences, including AtDUF966 and OsDSR, were used for training the parameter settings. Then, the trained parameters were used to identify conserved TaDUF966 motifs as follows: each sequence may contain any number of non-overlapping repeats of each motif; number of different motifs (1–10); motif width = 6–50 amino acids (aa). The functions of predicted motifs were analyzed using InterPro (http://www.ebi.ac.uk/interpro) and SMART motif search programs (http://coot.embl-heidelberg.de/SMART). In addition, multiple sequence alignment analysis of TaDUF966 proteins was conducted using ClustalX2 (Kim and Joo, 2010).

Prediction of Cis-Acting Elements

To analyze putative cis-elements in TaDUF966 gene promoters, the 1,500 bp sequence upstream of the start codon of each gene was searched in the Plant-CARE database (Bailey and Elkan, 1994). The identified motifs were analyzed using the Fisher test, and those with an adjusted p-value < 0.05 were considered significantly enriched.

Analysis of Wheat Transcriptome Under Salt Stress

Transcriptome responses of two wheat cultivars to salt stress RNA-seq datasets generated previously under an abiotic stress treatments were downloaded from NCBI and mapped to the wheat reference genome using HISAT2. Then, genes were assembled using cufflinks to inspect the expression levels of TaDUF966 genes (normalized by TPM, transcripts per kilobase of exon model per million mapped reads). The R package “pheatmap” was used to draw the pheatmap of TaDUF966 genes (Zhu et al., 2019c).

RNA Extraction and Quantitative Real-Time PCR (qRT-PCR) Analysis

Total RNA was isolated from root samples from plants treated with or without NaCl using the TRIzol Reagent (Invitrogen, Carlsbad, CA, U.S.A). The isolated total RNA samples were quantified using the NanoDrop spectrophotometer (Thermo Fisher Scientific) and then treated with DNase I (Thermo Fisher Scientific) to remove any residual DNA contamination. The DNA-free RNA samples were converted to cDNA using a cDNA synthesis kit (Thermo Fisher Scientific). Gene-specific primers were designed using Primer Premier 5.0 (Table S1). Then, quantitative PCR (qPCR) was carried out using Maxima SYBR Green/ROX qPCR Master Mix (Zhang et al., 2019; Chen et al., 2020). The EF-1α gene was used as an internal reference. Three independent biological replicates, each containing three technical replicates, were conducted for each sample in both control and NaCl treatments. Gene expression was quantified using the 2^−△△Ct method (Yin J. L. et al., 2018).

Virus-Induced Gene Silencing (VIGS) Experiment in Wheat

Two cDNA fragments from different sites silenced TaDUF966-9B. Using BLAST analysis in NCBI, these fragments are specific (Zhu et al., 2018a). Primers for constructs in plant transformation were designed using Primer Premier 5.0 and are listed in Table S2. Barley stripe mosaic virus (BSMV) consists of α, β, and γ three RNA chains. According to the manufacturer’s instructions of RiboMAXTM, using the RiboMAXTM large-scale RNA production system-T7 (Petty et al., 1990) and Ribom7G Cap analogues (Promega, Madison, Wisconsin, USA) in vitro linearized plasmids containing the tripartite BSMV genome were sealed at the transcripts. The second leaf of the two-leaf wheat seedling was mechanically wiped with BSMV transcript and incubated at 25°C. Use BSMV: TaPDS (TaPDS, wheat octahydrolysin desaturase) as a positive control.

Control plants were treated with 1× Fes buffer (0.1M glycine, 0.06M K2HPO4, 1% w/v tetrasodium pyrophosphate, 1% w/v bentonite, and 1% w/v celite, pH 8.5) devoid of BSMV transcripts (Zhu et al., 2018b). At 10 dpi, 150 mM NaCl was poured and samples were taken at 1 d and 6 d for RNA isolation, followed by qRT-PCR analysis and phenotypic analysis. qRT-PCR primers were designed using Primer Premier 5.0 and are listed in Table S3. The experiment was repeated at least three times.

Identification and Classification of TaDUF966 Genes

A total of 28 TaDUF966 proteins were identified from wheat reference genome using seven rice and five Arabidopsis protein sequences as query (Table S4). These wheat proteins were named according to the following rules: (1) “TaDUF966” stands for the hexaploid wheat DUF966 gene family; (2) Arabic numerals represent the gene number; (3) “-A/B/D” represents the wheat subgenome to which the gene belongs [Table 1; (Song et al., 2019)]. We identified 7 triads TaDUFs (Table S5). To better analyze the phylogenetic relationships among wheat and other plant DUF966 proteins, we constructed a phylogenetic tree using DUF966 sequences from wheat, rice, and Arabidopsis (Figure 1A). All DUF966 proteins were divided into two major groups (Groups I and II) (Figure 2A). Each group contained 14 members, and TaDUF966 sequences were present in both groups. Group II contained DUF966 proteins from only wheat and rice, indicating that DUF966 genes in wheat were more closely related to those in rice (Figure 1A).

Chromosomal Location of TaDUF966

The identified 28 TaDUF966 genes were distributed on 13 chromosomes, and each chromosome contained 1–4 TaDUF966 genes (Figure 1B and Table 1). The 24 genes of the TaDUF966 gene were located on chromosome 3 (3A, 3B, and 3D, 11 genes), chromosome 4 (4A, 4B, and 4D, 9 genes) and chromosome 7 (7A, 7B, and 7D, 4 genes). The other 4 chromosomes (Chr1A, Chr2B, Chr5B, and Chr6B) contained only one TaDUF966 gene for each chromosome. Most members of Group I were distributed on 3 chromosomes (3A, 3B, and 3D), while members of Group II were evenly distributed on 13 chromosomes.

Homologous Gene Pairs and Synteny Analysis

The homologous pairs of DUF966 were identified in the common wheat ancestral species, Triticum urartu (AA), Triticum dicoccoides (AABB), and Aegilops tauschii (DD); 10, 15, and 9 homologous pairs of DUF966 were identified, respectively, named as TuDUF966, TdDUF966, and AeDUF966 based on the naming rules. The location and other details of these genes are summarized in Table 1. Comparisons between TuDUF966 vs. TdDUF966, TaDUF966 vs. TdDUF966, TaDUF966 vs. TuDUF966, TaDUF966 vs. AeDUF966, and TaDUF966 vs. TaDUF966 revealed 2, 21, 12, 6, and 21 pairs of putative paralogs, respectively (Figure 1C). To better understand the evolutionary factors affecting the DUF966 gene family, we calculated the Ka/Ks ratios among homologous pairs of TaDUF966, TuDUF966, TdDUF966, and AeDUGF966 genes (Table S6). Most of the DUF966 gene pairs showed Ka/Ks ratio < 1, indicating that these genes contained more synonymous than non-synonymous changes. Moreover, these genes were under negative selection pressure to protect the state of the ancestor. Three homologous pairs, including TaDUF966-11B/TuDUF966-9A, TaDUF966-11B/TuDUF966-8A, and TaDUF966-11B/TaDUF966-13B, showed Ka/Ks ratio > 1, indicating positive selection, accelerated evolution, and neofunctionalization. Analysis of duplication events of TaDUF966 genes revealed that gene duplication occurred mainly among Group II genes (Figure 1B, Table S6). Of the 28 TaDUF966 genes, 17 (60.7%) genes showed segmental duplication events. Furthermore, TuDUF966 genes showed segmental and tandem duplications, while TdDUF966 genes showed segmental duplications (Table S6). These data suggest that segmental duplication events played a key role in the expansion of the DUF966 gene family in wheat.

Analysis of Gene Structure and Motif Composition in TaDUF966 Genes

To investigate the structure of TaDUF966 genes, we extracted their detailed information from the GFF3 file of the reference genome (Jiang et al., 2020b). The results showed that all TaDUF966 genes contain introns. Comparative exon–intron structure analysis revealed the presence of 2–6 exons in Group I members and 4–8 exons in Group II members (Figure 2B). Motif usually refers to the basic structure that constitutes any kind of characteristic sequence. It is a subset of the structural domain, and its function is to reflect the various biological functions of the structural domain. The prediction of protein motifs is a useful protein analysis method. Analysis of 28 TaDUF966 proteins using MEME revealed 10 motifs (Figures 2C, D). Details of these 10 motifs are shown in Table S7. Each TaDUF966 protein contains 3–10 motifs, and most TaDUF966 proteins contain motifs 1, 2, and 3. In the phylogenetic tree, we found that members with relatively similar genetic relationships have more similar motifs, which indicates that DUF966 members gathered in the same subgroup may have more similar biological effects. In addition, all TaDUF966 proteins contain the DUF966 conserved Motif 1. Motif analysis using NCBI BLAST showed that seven of these 10 motifs (motifs 1–5, 9, and 10) belong to the DUF966 conserved domain, while the remaining three motifs were uncharacterized. Multiple sequence alignment analysis of TaDUF966 proteins revealed that the amino acid sequence of the DUF966 domain is highly conserved (Figure S1), suggesting that this domain plays a key role in protein function.

Analysis of TaDUF966 Protein Features

Protein feature analysis showed that TaDUF966 proteins contain an average of 422 aa (191 aa–644 aa) (Table 2). Other features of TaDUF966 proteins were as follows: average pI = 8.89 (range = 6.46–9.88); average MW = 46.1 kDa (range = 21.3–69.84 kDa); average instability index = 63.48 (range = 36.76–79.99) (Table S8). The instability index values indicate that TaDUF966 proteins are relatively unstable. The results of the subcellular localization prediction of all 28 TaDUF966 proteins showed that 27 proteins localized to the nucleus or chloroplast, whereas one protein localized to the cytoplasm.A

Prediction of Cis-Acting Regulatory Elements in TaDUF966 Gene Promoters

Cis-acting elements play important roles in the initiation of gene transcription. The TaDUF966 gene promoters were predicted to possess cis-acting regulatory elements, which are involved in three different physiological and biochemical processes. (Figure 3, Tables S9, S10). Among the biotic/abiotic stress-responsive cis-elements, ARE (antioxidant response element) and G-box were highly enriched in TaDUF966 gene promoters. Secondly, the CAAT-box and TATA-box (cis-elements related to growth and development) were also enriched in TaDUF966 gene promoters. Among the phytohormone-related cis-acting elements, the ABA-responsive element (ABRE) was the most frequently identified in TaDUF966 gene promoters. In addition, Cis-acting regulatory elements involved in the response to methyl jasmonate (MeJA), CGTCA-box and TGACG-box, were also identified in most of the TaDUF966 gene promoters. In general, TaDUF966 gene promoter may respond to various endogenous and exogenous stimuli.

Transcriptome Expression Profiling of TaDUF966 Genes Under Salt Stress

Gene expression acts as a link between the transmission and realization of genetic information in organisms, and genome-wide analysis of gene expression patterns helps to identify the gene expression pathway and its regulation. In this study, RNA-seq data downloaded from NCBI were mapped to the wheat reference genome using HISAT2, and expression patterns of TaDUF966 genes under salt stress were analyzed. The pheatmap, based on standardized TPM values, showed that salt stress strongly induced TaDUF966-5, -9, and -14 (A, B, and D) but did not affect the expression of other genes (Figure 4A, Table S11) (Appels et al., 2018). The majority of TaDUF966 genes showed low expression, whereas 9 genes showed high expression in two different wheat varieties. The expression of TaDUF966 genes increased with the duration of salt treatment, reaching a peak at 6 h, followed by a gradual decline.

Expression Analysis of TaDUF966 Genes

To further reveal the potential functions of TaDUF966 genes in association with abiotic stress (salt stress), we analyzed the expression patterns of several TaDUF966 genes by qRT-PCR. According to transcriptome analysis, expression of TaDUF966-5, -9, and -14 (A, B, and D) was relatively higher in salt-treated roots compared with control roots. The results of qRT-PCR analysis of TaDUF966-5D, TaDUF966-9B, and TaDUF966-14D were consistent with RNA-seq data. These three genes showed higher expression levels in the root and leaf tissues of wheat plants treated with NaCl compared with the control. Expression levels of these three genes increased in leaf tissues at 4 or 8 h after NaCl treatment but decreased in root tissues at 8 h (Figure 4B). In particular, the expression level of TaDUF966-9B increased in root tissues at 2 h after NaCl treatment, decreased in root tissues at 8 h, and was induced by salt stress 24 h later. Therefore, TaDUF966-9B may be better induced by salt stress. In subsequent trials we focused on this gene.

Silencing of TaDUF966-9B Reduces Salt Tolerance in Wheat

To further confirm the function of TaDUF966-9B, we performed a virus-induced gene silencing (VIGS) assay. Silencing of the TaDUF966-9B gene was achieved using two constructs (BSMV: TaDUF966-9B-1 and BSMV : TaDUF966-9B-2), each carrying an approximately 250 bp fragment (Figure S2) of the corresponding TaDUF966-9B gene, respectively, applied as independent treatments. All plants inoculated with BSMV showed mild chlorotic mosaic symptoms at 10 days post-inoculation (dpi), while leaves inoculated with BSMV : TaPDS showed photobleaching (Figure 5A), indicating that the BSMV-induced gene silencing system was functional. The expression of each homolog was also analyzed by qRT-PCR. At 10 days after friction-implanted virus, the leaf curling phenotype was significantly stronger in plants expressing either construct than in plants carrying the empty vector control (Figures 5B, C), and the level of BSMV : TaDUF966-9B-1/2 mRNA in NaCl-treated plants was significantly lower than that in the empty vector control (Figure 5D). These results indicate that the TaDUF966-9B gene is involved in salt stress tolerance in wheat.