The DUF506 Gene Family in Triticum aestivum: Genome-Wide Identification and Expression Profiling Under Salt Stress

Abstract

As a subfamily of the PD-(D/E)XK nuclease superfamily, DUF506 family shows great potential in abiotic stress responding of higher plant, yet its clues of structure, evolution and functions remain largely unexplored due to their distant phylogenetic relationship with other nuclease families, especially in Triticum aestivum. In this study, 26 T. aestivum DUF506 genes (TaDUF506) were identified from genome-wide level through bioinformatic techniques. Phylogenetic and structural analyses revealed that TaDUF506 genes exhibit conserved motif and gene structure patterns intra-phylogenetic clusters but display significant divergence inter-clusters. Gene duplication identification showed that whole-genome duplication event (WGD) was the primary driver of TaDUF506 family expansion, while Ka/Ks analysis indicated that whole TaDUF506 family experienced purifying selection generally. Gene ontology analysis and protein-protein interaction prediction suggested that DUF506 plays a potential role in transcription regulation and nucleotide-excision generally. Promoter analyses highlighted an enrichment of hormone-responsive elements linked to salt stress in TaDUF1.3-3D TaDUF5.1-3A, with expression analysis demonstrated their significant upregulation under salt stress, suggesting the potential roles in stress responses. Altogether, our study advances the understanding of DUF506 gene family in higher plant from structural, evolutional and functional aspects, and thereby provides a foundation for the development of salt-tolerant wheat varieties.

1. Introduction

Triticum aestivum (T. aestivum) is a critical cereal crop, which provides approximately 19% of the global dietary energy supply, widely used in animal feed and processed foods. However, abiotic stresses, particularly drought and salinity, have severely limited T. aestivum growth and yield, necessitating a deeper understanding of the molecular mechanisms underlying wheat’s stress responses [1].

“DUF” (Domains of Unknown Function) originally represents conserved but functionally unknown domains in Pfam database [2]. Currently, a plenty of studies have characterized DUF domain-containing proteins widely involved in participating abiotic stress response. For example, a DUF569 domain-containing protein in Arabidopsis is a positive regulator to drought, while some of DUF668 domain-containing proteins in rice have been shown to response salt and wound [3].

PD-(D/E)XK nucleases superfamily includes a series of DUF domain-containing protein family, such as DUF524, DUF506 and DUF1626 [4]. With great structural and functional diversity, PD-(D/E)XK nucleases mainly participating in DNA degradation, repairment and recombination [5]. As a subfamily of PD-(D/E)XK superfamily, members in DUF506 show low sequence similarity and distant homologous relationship with other PD-(D/E)XK nuclease proteins, causing difficulties to infer precise functions of single member [6]. Recent studies have identified part of DUF506s as potential abiotic stress response genes Arabidopsis and rice. For example, AtRXR3 is an Arabidopsis DUF506 (AtDUF506) protein that playing a role in root hair growth suppression under phosphate-limiting conditions [7], and most of DUF506 in Oryza sativa (OsDUF506) respond to ABA, JA and cold condition [8], underscoring the functional versatility of DUF506 family in higher plants. However, there remains a huge missing clue of DUF506 family in T. aestivum, which constrains our understanding of complete PD-(D/E)XK superfamily in higher plants, especially their potential role in T. aestivum.

In this study, 26 TaDUF506 genes were identified from T. aestivum genome through comprehensive bioinformatic analyses. We systematically investigated their evolutionary, structural, and functional properties of TaDUF506 members, as well as their behaviors under salt stress. Our findings not only filled the missing clue of DUF506 in T. aestivum. providing insights for understanding PD-(D/E)XK nuclear superfamily in higher plants, but also reveal potential value of TaDUF506 family in breeding salt-tolerance T. aestivum cultivars.

2. Materials and Methods

2.1. Genome-Wide Identification and Validation of the TaDUF506 Gene Family

The complete T. aestivum proteome was obtained from the Ensembl Plants database (https://plants.ensembl.org/index.html, accessed on 27 December 2024) [9], and the annotation file for Hidden Markov Model (HMM) search was obtained from the Pfam database (http://pfam.xfam.org/, accessed on 27 December 2024) (This database has been integrated into InterPro) [10,11]. DUF506 domain (Pfam ID: PF04720, ID: PDDEXK_6) was used for predicting candidate DUF506 proteins in T. aestivum through a Hidden Markov Model (HMM) search conducted by TBtools software (TBtools-ll (Toolbox for Biologists) v2.152), applying a threshold E-value of <0.01. Subsequently, InterPro database was used for checking the existence of DUF506 domain in searched proteins [11]. We also collected 13 DUF506 proteins in Arabidopsis by Ying et al. [7] for BLASTp alignment (E-value of <1 × 10⁻¹⁰), based on their conserved domain region, to cross-validate the results.

2.2. Chromosome Localization and Phylogenetic Analysis

The Triticum aestivum genome and corresponding GFF3 gene annotation files (Version 59, accessed by 27 November 2024) were downloaded from the Ensembl Plants database (https://plants.ensembl.org/index.html, accessed on 27 December 2024) [9]. Using TBtools [12], the chromosome length, gene density, and chromosomal coordinates for the 26 identified TaDUF506s were extracted and visualized. Phylogenetic analysis was conducted on the DUF506 members in Triticum aestivum and its ancestral species. The conserved domain sequences of the primary transcript-encoded proteins were aligned using MUSCLE [13]. The resulting alignment was used to construct a ML tree with a bootstrap value of 1000 iterations for reliability based on iqtree2 (iqtree2, version 2.3.6, http://www.iqtree.org/, accessed on 27 December 2024) [14]. The phylogenetic tree visualization was refined and finalized using the Interactive Tree of Life platform (iTOL, version 6.7.1, https://itol.embl.de, accessed on 27 December 2024) [15].

2.3. Gene Replication Events and Collinearity Analysis

Gene duplication events within the TaDUF506 family were analyzed using TBtools [12] and the Multiple Collinearity Scan Tool (MCScanX). Sequence alignments via DIAMOND, with parameters set to E <1 × 10⁻³, facilitated whole-genome comparisons across T. aestivum and its ancestral species to identify duplication events [16]. Following collinearity analysis, TBtools was employed for result simplification and visualization. To further investigate evolutionary patterns and estimate the divergence time of TaDUF506 genes, TBtools was used to compute the non-synonymous (Ka) and synonymous (Ks) substitution rates for duplicated TaDUF506 gene pairs. Selection pressures on the TaDUF506 family were inferred from Ka/Ks ratios, and divergence time was estimated using the formula T = [Ks/2 × (6.5 × 10⁻⁹)] × 10⁻⁶ (million years) [17,18].

2.4. Analysis of Protein Characteristics and Gene Structure

The ProtParam online tool (https://web.expasy.org/protparam/, accessed on 27 November 2024) [19] was utilized to predict the protein sequence length, instability index, isoelectric point, molecular weight, and hydrophobicity index. For conserved motif analysis, the MEME suite (Version 5.5.7, https://meme-suite.org/meme/tools/meme, accessed on 27 November 2024) [20] was applied with parameters set to identify motifs 6–50 amino acids in length, allowing up to 12 motifs. Gene structure visualization was performed using TBtools [12], with the predominant transcripts of TaDUF506 family members selected and mapped based on GTF genome annotation files (Version 59, accessed by 27 November 2024).

2.5. Cis-Regulatory Element Prediction and Gene Expression Analysis

To investigate the promoter regions of the TaDUF506 gene family, TBtools was utilized to extract the upstream 2000 bp sequences of all 26 TaDUF506 gene promoters. These sequences were then submitted to the PlantCARE database (https://bioinformatics.psb.ugent.be/webtools/plantcare/html, accessed by 27 November 2024) for the prediction of cis-acting regulatory elements [21]. The results were visualized by generating a heatmap in TBtools. To further analyze the expression patterns of TaDUF506 family members under salt stress, expression data for the 26 genes in both Chinese Spring and QingMai 6 T. aestivum cultivars under salt stress conditions were obtained from the WheatOmics database (http://wheatomics.sdau.edu.cn/, accessed by 27 November 2024) [22]. This data was also visualized through heatmap generation in TBtools.

2.6. Plant Cultivation and Tissue Sample Collection

Chinese Spring wheat seeds were sown in 5 cm × 5 cm flower pots containing 10 g of peat and placed in a growth chamber with controlled conditions (22 °C, 60% relative humidity, 400 ppm CO₂, and a 12-h light/dark cycle). After one week, the plants were transferred to a 0.2 × Hoagland nutrient solution for hydroponic cultivation, with the solution replaced every two days. At the three-leaf stage, the T. aestivum was treated with 0.2 × Hoagland solution containing 150 mM NaCl, while a control group received the same solution without NaCl. Leaf and root samples were collected at 6, 12, and 24 h post-treatment, immediately frozen in liquid nitrogen, and stored at −80 °C for subsequent RNA extraction [18].

2.7. RNA Extraction and Real-Time Quantitative PCR (RT-qPCR) Analysis

Total RNA was extracted from T. aestivum roots and leaves using the Trizol reagent kit (Carlsbad, CA, USA, Invitrogen), according to the manufacturer’s instructions. RNA concentration and purity were measured using a NanoDrop spectrophotometer, and the RNA samples were subsequently reverse-transcribed into cDNA with a cDNA synthesis kit (Beijing, China, TIANGEN). The resulting cDNA was stored at −20 °C for further analysis. Gene-specific primers were designed using Primer Premier 6.0 software and are listed in Table S1 [23]. Primer specificity was verified through the Primer BLAST database (https://www.ncbi.nlm.nih.gov/tools/primer-blast/, accessed by 27 November 2024) [24]. RT-qPCR was then performed using the 2 × RealStar Universal SYBR qPCR Mix (Beijing, China, Beijing Kangrun Biotech), with β-actin as the internal control gene. Each sample was analyzed in triplicate, and data analysis was conducted using the 2^−ΔΔCt method. A paired t-test was made for analyzing significance of results, wherein p < 0.05 is considered as significant [18].

2.8. Subcellular Localization of DUF506 Protein in T. aestivum

The coding sequence (CDS) of the predominant TaDUF506 transcript was cloned into the pCAMBIA1302 vector and fused with the GFP gene for expression (Figure S1). After confirming the sequence of the recombinant vector, the construct was introduced into Agrobacterium tumefaciens strain GV3101 via electroporation and cultured at 30 °C for 48 h. The bacterial suspension was then resuspended in 10 mM MgCl₂ containing 120 μM acetosyringone and injected into the lower epidermal cells of healthy Nicotiana benthamiana plants. The plants were grown under a 14-h light/10-h dark photoperiod for one month, followed by maintenance under low-light conditions for 48 h post-transfection. For visualization, leaf segments were mounted on glass slides with water droplets and examined using a laser confocal microscope (Tokyo, Japan, Nikon, model C2-ER) [25].

3. Results

3.1. Identification and Phylogenetic Analysis of TaDUF506 Family Members

To identify TaDUF506 members, HMM search combined with BLASTp analysis were performed, and 26 TaDUF506s were identified (Table S2). All TaDUF506s distributed evenly across the three sub-genomes of T. aestivum, with the highest density observed on chromosome 3 (Figure 1A).

Based on full-length of protein sequences, phylogenetic relationship between TaDUF506 and DUF506 in other plants were analyzed using ML method. All TaDUF506 members were classified into five distinct clusters according to branch length and tree topology, as shown in Figure 1B. Interestingly, the correlation between sub-genome number and counts of DUF506 members presents as a linear relationship in T. aestivum and its ancestral species (Figure 1C), suggesting a possibility on how genome duplication events have influenced the amplification and diversification of TaDUF506 gene family.

The physicochemical properties of the dominant transcript-encoded proteins from identified TaDUF506 genes were further analyzed (Table S3). The length of coding sequences of the TaDUF506 proteins ranged from 723 to 1146 bp, with corresponding amino acid sequences ranged from 240 to 381 residues. The isoelectric points (pI) ranged from 6.13 to 9.32, indicating that most TaDUF506 proteins are neutral or alkaline. The instability index, which spans from 39.16 to 69.23, suggests that the majority of TaDUF506 proteins are classified as unstable. Hydrophilicity indices range from −0.103 to −0.581, showing that all TaDUF506 proteins exhibit hydrophobic properties to varying extents. Notably, several physicochemical attributes appear consistent among members within the same phylogenetic cluster, implying that phylogenetically close proteins tend to share similar physicochemical characteristics.

3.2. The Gene Structure and Conserved Motifs of TaDUF506 Family Members

The conserved motif prediction of the TaDUF506 protein family was performed using the MEME online platform, identifying a total of 10 conserved motifs (Figure 2A). Motif 1 through 5 were consistently present in all TaDUF506 protein members, with a conserved sequence arrangement of Motif 2-Motif 3-Motif 1 across all members. This pattern likely represents a key sequence feature of the DUF506 functional domain (Figure 2B). In contrast, Motif 6, Motif 7 and Motif 9 through 12 were restricted to specific phylogenetic clusters within the TaDUF506 family. For instance, Motif 10 appeared exclusively in group 1, while Motif 8 was found only in group 2, suggesting potential functional divergence across the phylogenetic types.

The exon-intron structure is a critical aspect of gene characteristics and evolutionary dynamics. Analysis based on the published GTF genome annotation file of Triticum aestivum revealed that the each of identified TaDUF506 contains between 1 to 3 exons and up to 2 introns (Figure 2C). Most TaDUF506 genes possess both 3′- and 5′-UTR regions, though a few only feature a UTR on one end or lack UTR regions altogether. Additionally, members within the same phylogenetic cluster display a similar exon-intron organization, suggesting that genes encoding proteins with similar conserved domains tend to exhibit comparable structural patterns (Figure 2C).

3.3. Gene Duplication Events and Collinearity Analysis of TaDUF506 Family

Collinearity analysis provides valuable insights into gene family expansion in higher plants [26,27]. Among the TaDUF506 family members, 34 gene duplication events were identified, all of which were classified as whole genome duplications (WGD) or segmental duplications (Figure 3A). The duplication events observed between different TaDUF506 members were showing a symmetrical pattern within three sub-genomes that consistent with their phylogenetic classification, indicating a hallmark of WGD event.

To assess the evolutional pressure on the TaDUF506 gene family, we analyzed the ratio of non-synonymous to synonymous substitution rates (Ka/Ks) for duplicated gene pairs (Table S4). All duplicated gene pairs had Ka/Ks ratios below 1, suggesting an occurrence of significant purifying selection. Additionally, the estimated divergence times for these gene pairs were predominantly within the last 10 million years (Mya). However, in the group 2 phylogenetic cluster, certain TaDUF506 members showed divergence times exceeding 25 Mya, implying a potentially more ancestral origin within the family.

To investigate gene duplication events among TaDUF506 family members and DUF506 families from other species, we analyzed the orthologous relationships between TaDUF506 and the DUF506 families in Triticum dicoccoides, Triticum turgidum, and Triticum urartu. A total of 3, 35, and 25 gene duplication events were identified respectively (Figure 3B). Notably, the collinearity relationships do not exhibit a simple one-to-one pattern between the TaDUF506 and the DUF506 families in these species. Instead, multiple TaDUF506 members show redundant collinear relationships with the same orthologous genes in these other species. This observation aligns with the common presence of symmetrical collinear relationships among the three sub-genomes of Triticum aestivum, which likely resulted from gene duplication events within the TaDUF506 family.

3.4. Functional Prediction, Expression Profiling, and Cis-Regulatory Element Analysis of the TaDUF506 Family

GO enrichment analysis of the TaDUF506 family using AgriGO database revealed significant associations with several key molecular functions (Figure 4A) [28]. The RNA polymerase activity (GO:0034062) and DNA-directed RNA polymerase activity (GO:0003899) indicate that TaDUF506s may play a crucial role in transcriptional regulation, particularly in gene expression. Additionally, nucleotidyl transferase activity (GO:0016779) suggests involvement in RNA synthesis or modification processes, potentially influencing gene expression. The enrichment of transferase activity transferring phosphorus-containing groups (GO:0016772) points to the role of these genes in phosphorylation-based signaling pathways, which are vital for stress responses and cellular regulation.

RNA-seq analysis revealed distinct expression patterns of TaDUF506 genes under salt stress. The raw RNA-Seq expression data of DUF506 genes in wheat under control and salt stress conditions are presented in Table S5 and visualized through a heatmap for a clearer comparison of expression patterns in Figure 4B. We observed a significant upregulation of some TaDUF506s in phylogenetic cluster 1 (TaDUF1.3-3D, 1.3-3B & 1.8-6D) and in phylogenetic cluster 5 (TaDUF5.1-3A, 5.2-3B & 5.3-3D), indicating specialized roles of these TaDUF506s in salt stress response.

Cis-regulatory elements are crucial for the transcriptional regulation of genes. We identified a total of 27 cis-acting elements within the TaDUF506 family (Figure 4C). These elements were categorized into the following groups: Common elements, Phytohormone-responsive elements, Environment-responsive elements, Transcription Factor (TF)-associated elements, and Tissue-specific elements. The high abundance of key cis-regulatory elements, including A-box, ABRE, TATA-box, and CAAT-box, within the TaDUF506 family suggests a complex regulatory network driving the expression of these genes, particularly in response to environmental stimuli.

3.5. Protein Interaction Prediction of TaDUF506 Members

To explore potential regulatory patterns via protein-protein interaction (PPI) of TaDUF506 family members, we conducted PPI predictions using the STRING database with default parameters [29] (Figure 5A), with GO enrichment analysis performed to investigate their general functions (Figure 5B). These results indicating that proteins interacting with TaDUF506 may play a potential role in environmental stress response and nucleotide excision repair.

3.6. qRT-PCR Validation and Subcellular Localization

Under 150 mM NaCl treatments, 2 members of the TaDUF506 family exhibited distinct response patterns in different tissues (Figure 6A, Table S6). Specifically, TaDUF1.3-3D showed a mild negative response in the leaves, although not significant, while it exhibited a clear positive response in the roots. After 6 h of salt stress, its expression increased significantly, more than 10 times its initial level. TaDUF5.1-3A displayed a positive response, with significant expression changes in both roots and leaves, except at 6 h.

Subcellular localization results revealed that the coding protein of three TaDUF506 with duplication event (TaDUF5.1-3A, TaDUF5.2-3B, and TaDUF5.3-3D) were predominantly localized to plastids, where the green fluorescence signal was particularly strong. Although there was also some fluorescence in the cytoplasm, it was less intense compared to that in the plastids. These results suggest that these TaDUF506 proteins localize to the same subcellular compartment.

4. Discussion

The PD-(D/E)XK nucleases, originally identified within type II restriction endonucleases, have since been recognized as a large, highly diverse superfamily [4]. While the amino acid sequences of this superfamily often lack detectable similarity beyond their conserved active sites, these nucleases exhibit remarkable structural and functional diversit. Structurally, the PD-(D/E)XK superfamily features a conserved core with a four-stranded β-sheet flanked by two α-helices in an αβββαβ topology [4]. However, non-conserved peripheral elements often dominate the protein composition, complicating sequence-based predictions of function, which renders the PD-(D/E)XK nucleases one of the most complex families for protein sequence analysis and structural prediction. Within PD-(D/E)XK superfamily, DUF506 family represents a distantly homologous subgroup, causing difficulties in inferring molecular function of single DUF506 member. Recently, a systematically study focusing on AtDUF506 gene family has suggested their potential value in abiotic stress responses and other biological functions in plants [7], comparable analyses in T. aestivum remain scarce. In this study, we identified 26 high-confidence TaDUF506 genes at whole-genomic level, with a comprehensive analysis conducted to investigate their structures, evolutionary relationships and functions.

Except for chromosome 6A, TaDUF506 members with highly similar structures are evenly distributed among three sub-genomes, indicating potential occurrences of whole-genome duplication (WGD) event on TaDUF506 gene family (Figure 3A,B). Ka/Ks ratio analysis revealed significant selective purifying selection within the DUF506 gene family (Table S4), which likely reflects the family’s crucial role in the adaptive evolution of T. aestivum.

Protein-protein interaction (PPI) analysis, utilizing the STRING database, suggested that TaDUF506 family members interact with PD-(D/E)XK-like family proteins, implicating their involvement in essential biological processes, especially those associated with stress response (Figure 5A). Furthermore, PPI analysis revealed that TaDUF506 interacts with proteins involved in cellular stress response, implying a role in maintaining cellular homeostasis by protecting cells from oxidative stress and DNA damage. Additionally, interactions with proteins involved in nucleotide excision repair point to the family’s possible involvement in maintaining genome stability.

GO enrichment analysis indicated that these genes are closely associated with key stress response pathways (Figure 4A), while in TaDUF1.3-3D, 1.3-3B and 1.3-3A from phylogenetic cluster 1, we observed a significant enrichment of ABRE and G-box elements. Furthermore, RNA-seq data and RT-qPCR experiments also confirmed their positive responses to salt stress. In other higher plants, the expression of the Arabidopsis orthologous genes (At3g25240 and At2g39650) corresponding to these two groups of TaDUF506 members was also strongly induced by salt stress (log₂|FC| ≥ 2), while the rice ortholog (LOC_Os01g68650) of the three cluster 1 TaDUF506s positively responded to both abscisic acid signaling and drought stress [7,8]. Notably, a recent related study indicated that At3g25240 (AtRXR1) responds to phosphate starvation via a PHR1/PHL1-dependent pathway, whereas the segmentally duplicated gene At3g07350 of At3g25240 is involved in an auxin-mediated phosphate signaling pathway [30,31]. Therefore, even though the three cluster 1 TaDUF506 members that positively respond to salt stress share gene duplication events, they may function in a nonredundant manner.

Future studies should focus on potential pathway and mechanisms involved by salt responsive TaDUF506s, as well as how DUF506s functionally work based to their subcellular location in cell. Using CRISPR/Cas9 or RNA interference (RNAi) approaches can be conducted to assess the phenotypic impact of specific TaDUF506 genes under salt stress conditions, while co-immunoprecipitation (Co-IP) or yeast two-hybrid (Y2H) assays could be employed to confirm more precise protein-protein interactions of TaDUF506s. These approaches will provide a comprehensive understanding of the regulatory networks involving TaDUF506 genes and facilitate their application in T. aestivum breeding programs. Genetic engineering method could be applied for overexpressing salt responsive TaDUF506s by expanding their copy number, or interfere their downstream regulate genes after claiming more detail mechanisms.

5. Conclusions

In summary, this study systematically analyzed the TaDUF506 gene family in T. aestivum, providing valuable insights into their evolutionary, structural, and functional characteristics. Notably, our results highlight the significant roles of TaDUF506 genes in transcriptional regulation and stress adaptation, particularly under salt stress conditions. These findings offer practical implications for T. aestivum resilience and productivity, as the identified stress-responsive genes could serve as molecular targets for the development of salt-tolerant T. aestivum varieties. Future research focusing on functional validation and transgenic studies will further elucidate the potential of the TaDUF506 gene family in improving abiotic stress tolerance of T. aestivum, sustaining agricultural productivity in challenging environments.