Genome-wide identification and expression analysis of ADP-ribosylation factors associated with biotic and abiotic stress in wheat (Triticum aestivum L.)

Identification of ARF genes

Computer-based methods were used to identify members of the ARF gene family from the wheat reference genome IWGSC RefSeq v1.1 assembly (https://wheat-urgi.versailles.inra.fr). Known ARF protein sequences, including 19 ARFs from Arabidopsis thaliana (AtARFs) and 21 ARFs from Oryza sativa (OsARFs), were used as query sequences for BLASTp analysis with an e-value cutoff of <1 ×10⁻¹⁰ (Altschul et al., 1997). The obtained sequences were submitted to InterProScan (http://www.ebi.ac.uk/interpro/) to check the ARF domains Mulder & Apweiler (2007). The Pfam database (http://pfam.xfam.org/) was used to select sequences that contained the ARF-box domain (PF00025) (Finn et al., 2006). The hit sequences were further validated using SMART (http://smart.embl-heidelberg.de/smart) to remove unmatched proteins (Letunic et al., 2004).

Phylogenetic analysis and characterization of TaARFs

A total of 96 protein sequences (19 from Arabidopsis (Gebbie et al., 2005), 21 from rice (Muthamilarasan et al., 2016), and 56 from wheat) were compared by ClustalW2 software with default parameters (Oliver et al., 2005). An unrooted phylogenetic tree was created using the Neighbor-Joining (NJ) method with 1000 replicated-bootstraps in MEGA 7.0 software (Kumar, Stecher & Tamura, 2016). Then, the phylogenetic tree was edited by the Interactive Tree of Life (ITOL, Version 3.2.317, http://itol.embl.de/) to have the final illustration (Letunic, 2016). The identified TaARFs were used to perform protein characterization in ExPASy Server10 (https://prosite.expasy.org/) (Wilkins et al., 1999). The predicted protein features for each of the protein sequences were determined, including the length, molecular weight (MW), instability index, isoelectric point (pI), and aliphatic index. Plant-mPLoc (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/) was used to predict the cellular localization (Horton et al., 2006).

Exon-intron structure and motif analysis of TaARFs

Structural analysis was performed to identify the exon-intron structure of each TaARF gene using GSDS (http://gsds.cbi.pku.edu.cn/index.php/) (Hu et al., 2015). The software MEME (http://meme-suite.org/) was used to determine conserved TaARF motifs (Bailey et al., 2006). Motifs with E-values >0.001 were probably statistical artefacts rather than real motifs, and were excluded (Bailey et al., 2006). The parameters were employed as the following descriptions: the maximum number of motifs, 10; and the optimum width of each motif, between 6 and 50 residues. The motif prediction results were inserted into TBtools (https://github.com/CJ-Chen/TBtools/) to produce the illustrations (Chen et al., 2018). The locations of conserved domains were determined by SMART and visualized in MEME to reveal the diversification of TaARFs (Fang et al., 2020).

Chromosomal locations and synteny of TaARFs

The location information of TaARF genes were obtained from the reference genome in the IWGSC v1.1 database. Moreover, we generated chromosome locations using MapInspect version 1.0. The common tool “all against all BLAST search” was used to determine possible paralogous or orthologous sequences among wheat sequences with an E-value cutoff of 1e−10 and identity >80% (Gu et al., 2002). The R package “circlize” was employed to prepare a diagram showing the locations and homology relationships of TaARFs (Zuguang et al., 2014). The non-synonymous (Ka) and synonymous (Ks) substitution rates were calculated using DNA Sequence Polymorphism (DnaSP) 5.10 to analyze gene duplication events (Rozas, 2009). We calculated the Ka/Ks ratios for the replicated ARF gene pairs in wheat to explore whether Darwinian positive selection pressure has affected the functional proportion of replicated genes. We also compared the gene duplication events of ARF genes between wheat and Triticum dicoccoides, Aegilops tauschii, Arabidopsis, and rice. In general, when the ratio is greater than 1, the replicated gene is under positive selection, a ratio equal to 1 indicates genes under neutral evolution, and a ratio of less than 1 indicates genes under negative selection pressure (Zhang et al., 2006; Anton, Makova & Li, 2002).

Multiple conditional transcriptome analysis of TaARFs

Multiple RNA-seq data from different tissues, development stages, and treatments were downloaded from the NCBI Short Read Archive (SRA) database and mapped to the wheat genome using HISAT2. Cufflinks were used to perform gene assembly, expression level calculations, and identifications of differences in differentially expressed genes (Cole et al., 2012). The obtained transcripts per million (TPM) values reflecting the expression level of each gene were used to generate a heatmap of TaARFs using the R package “pheatmap” (Kolde, 2015). Triad expression analysis was carried out as described previously (Ramírez-González et al., 2018).

Plant materials and qRT-PCR

Hexaploid wheat (cultivar Emai 170; Triticum aestivum; AABBDD) was used to validate the expression patterns of selected candidate genes in all experiments. Seeds were surface sterilized with 1% hydrogen peroxide, germinated in an incubator at 28 °C for 2 d, and transferred to a greenhouse at 26 °C with a 16 h light/8 h dark cycle and relative humidity of 60–70%. Seven-day-old seedlings were treated with cold, ABA, and NaCl. Control seedlings continued to grow under standard conditions. Root tissue was collected 6 h and 12 h after treatment, including the set controls. All sample materials were quickly frozen in liquid nitrogen and stored at −80 °C before RNA extraction. Total RNA was isolated using the Trizol reagent (Invitrogen, Carlsbad, CA, USA). Quantitative real time (RT)-PCR was performed using the TaKaRa qRT-PCR system (RR047A). The gene-specific primers (Table S10) were designed using Primer 5.0 to amplify 80–350 bp fragments. The thermal cycling conditions were as follows: 95 °C for 30 s and 40 cycles of 95 °C for 5 s and 60 °C for 30 s. The relative quantity of target gene transcripts was calculated using the 2^−ΔΔct method with wheat β-actin and GAPDH as reference genes (Yin et al., 2018).

Cis-element analysis of putative promoter regions

The cis-elements in the promoter regions are related to gene expression patterns and functions (Anne-Laure, Adrien & Veitia, 2014). To investigate the cis-elements in the promoter regions of genes of interest, we downloaded 1.5 kb of the genomic DNA sequences upstream of the start codon corresponding to each gene from the hexaploid wheat database. The Plant CARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) were used to analysis the putative cis-elements in the promoter sequences (Magali et al., 2002).

Gene Ontology annotation and protein interaction network

A GO (Gene Ontology) database was used for functional annotation of the ARF genes using MAJORBIO CLOUD (https://cloud.majorbio.com/). GO annotations were mapped according to biological processes, molecular functions, and cellular components. TaARFs protein-protein interaction networks was assembled using the STRING tool (https://string-db.org/). All predicted TaARF proteins have been submitted to the STRING database. The minimum required interaction score was set to high confidence (0.900). The active interaction sources was setted come from ‘experiments’ and ‘databases’. The maximum number of interactors were no more than 5 on the first shell.

Identification and classification of ARF genes in wheat

In the present study, the entire hexaploid wheat genome was downloaded and used to construct a local database. After genomic retrieval, a total of 126 proteins that were similar to ARF were obtained from hexaploid wheat; however, only 77 wheat sequences were confirmed to be conserved in the ARF family domain (PF00025 and IPR006689) by Pfam and InterProScan. Sequences without complete conserved domains or with a length less than 170 aa were excluded. After elimination, 74 ARF proteins containing complete domains were obtained. Among the remaining ARF proteins, 18 were splice variants (Table S1). To investigate their evolutionary relationships, we constructed an NJ tree with MEGA7.0 using the amino acid sequences of putative ARF family members from Arabidopsis, rice, and wheat (Fig. 1, Table S2). The predicted TaARF proteins were classified into seven sub-groups: ARFA1, ARFB1, ARFC1, ARF3, ARL1, TTN5, and GB1 based on the phylogenetic tree and previous reports (Muthamilarasan et al., 2016). The ARL1 sub-group is the largest of these sub-groups with 17 members. The ARFA1 and the ARFB1 sub-groups have 14 and nine members, respectively. There are three proteins in the ARF3 and ARFC1 sub-groups, and the other sub-groups each contain five members. Interestingly, Arabidopsis had one ARFD1 protein, while rice and wheat did not have any. These results are consistent with the coevolutionary relationships between these species, indicating that the relationship between wheat and rice is closer than that with Arabidopsis. Each TaARF gene was named based on its phylogenetic relationship with AtARFs and OsARFs (Muthamilarasan et al., 2016). Genes corresponded equally across the three homoeologous subgenomes (A, B, and D) in wheat, which is referred to as a triad (Ramírez-González et al., 2018). We identified 15 TaARFs triads with reference to the results of Ramírez-González et al. (2018) (Table S3). TaARFs triads have identical gene names except for the sub-genome identifier (A, B or D).

Features of predicted TaARF proteins

Of the 56 putative TaARF proteins, the predicted MW was around 21 kDa on average, except for TaARFA1a-D and TaARFB1c-B. The number of exons and introns were similar within the same class but different between classes. The protein lengths ranged from 170 (TaGB1b-B) to 263 (TaARFA1a-D) amino acids (aa) and the predicted isoelectric point ranged from 5.67 to 9.35 (Table 1). The protein instability index shows that 82% TaARFs are stable proteins, but TaARFA1a-D, TaARFB1b-A, TaARFB1c-B, TaARFB1d-B, TaARFB1e-N, TaGB1a-A, TaGB1a-B, TaGB1a-D, TaGB1b-B, and TaGB1c-D were predicted to be unstable. The average hydropathicity (GRAVY) was less than 0, except in the case of TaARFC1-A/B/D, indicating that most of these proteins are hydrophilic (Table S1). Subcellular localization prediction showed that TaARFs are localized mainly in the cytoplasm, but are also found in the chloroplasts, peroxisome, nucleus, or mitochondria (Table S1).

Gene structure and motif analysis of the TaARF gene family

The structure of genes can be used to predict their expression and function. We found that TaARFs contained different exon-intron composition patterns by comparing the gene structures (Fig. 2). There were 38 TaARF sequences with both untranslated regions (UTRs), and of the remaining sequences, 14 did not have 5′- and 3′-UTRs, 2 (TaARFA1d-A and TaARFA1c-D) had only 5′-UTRs, and 2 (TaARFA1a-D and TaARFB1b-A) had only 3′-UTRs (Fig. 2B). The number of introns ranged from 1 to 7, and 43% contained 6 exons of varying lengths. TaARF genes in the same sub-group shared similar exon-intron structure. The majority of the ARFB1 sub-group contained 6 exons, with the exception TaARFB1d-B, which contains five. While the ARFA1 subgroup gene members have 5 exons, ARFB1, ARFC1, and ARF3 subgroups have 6, 2 and 8 exons, respectively.

The 10 most statistically significant motifs were chosen to describe the motif pattern in TaARFs, which were named motif 1-10 (Fig. 2C, Table S4). The lengths of the 10 motifs were between 8 (motif 8 and 9) and 49 (motif 1) aa residues. The number of motifs in each TaARF protein varied from 5 to 8. Notably, motif analysis indicated that most TaARFs have relatively conserved motif compositions. All TaARF proteins (except TaGB1c-D and TaARFB1d-B) contained motif 1 to 4, which contained the ARF-box domain. In general, many TaARFs of the same group encoded proteins with similar motif compositions, and therefore, these proteins may have similar functions. There were 4 motifs (motif 1, 3, 4, and 5) in the middle region of most TaARF proteins; however, different motifs were found at the N-terminal and C-terminal regions. For example, two specific motifs (motif 9 and motif 10) were only found in the C-terminus of the TaARL1 sub-group. Motif 7 appeared in both the N-terminus and the C-terminus of the TaARL1 subfamily, but only in the C-terminus of other subfamilies.

Chromosomal distribution and gene duplication events of TaARF genes

To further investigate the genetic differences in the TaARF gene family, we mapped their chromosomal locations. After positioning, 55 TaARFs were mapped on all 21 chromosomes, while TaARFB1e-N was not mapped on any chromosomes (Fig. 3A, Table 1). TaARFs were distributed roughly evenly across the three subgenomes (subgenome A, 18; subgenome B, 18; subgenome D, 19). Members from the same sub-group tended to be distributed at similar locations. However, the distribution of genes on chromosomes varied from one homoeologous group to another. The largest homoeologous group of 5 chromosomes (A, B, and D) had the most TaARF genes (12), followed by the smallest homoeologous groups of 2 (4 genes) and 7 (5 genes). By screening the sequence identity and position, 75 segmental duplication pairs were predicted, with no tandem duplication pairs being identified. Triads were not considered when predicting gene duplication events. The Ka/Ks ratio between TaARF gene pairs was less than 1 (average 0.08) in all cases (Fig. 3B, Table S5), suggesting that the TaARF gene family might have experienced strong purifying selective pressure (Song et al., 2019). Among the TaARF genes, 75 pairs of segmental duplication genes were found concentrated in the ARFA1 subgroup (Table S5).

Four comparative syntenic maps of wheat associated with four representative species, including Triticum dicoccoides, Arabidopsis, Aegilops tauschii, and rice, were constructed to further deduce the evolutionary origin and orthologous relationship of the wheat ARF family (Fig. 4). The numbers of orthologous pairs between the other four species (Triticum dicoccoides, Aegilops tauschii, Arabidopsis, and rice) were 61, 27, 67, and 88, respectively (Table S6). Six TaARF genes have both orthologous genes in four species. Among them, the ARFA1 subgroup accounted for 5 genes, which was a larger number than the other subfamilies.

Expression analysis of TaARF genes in different tissues

RNA-sequencing (RNA-seq) is a powerful tool to explore transcription patterns using high-throughput sequencing (Wang, Gerstein & Snyder, 2009). The RNA-seq data for five tissues (grain, spike, leaf, stem, and root) in wheat were used to characterize the expression of TaARF genes during growth and development. Out of the 56 full-length genes, 93% were expressed in at least one developmental stage, with a wide expression range between 1-852 TPM (TPM_max) (Fig. 5A; Table S1, Table S7). The remaining 7% of full-length genes had a very low level of expression with a TPM_max <  1, and were considered as “not expressed” (Fig. S1, Table S1). In general, the expression of TaARFs could be divided into three patterns: the first group contains members that are widely expressed in many tissues during multiple developmental stages, the second group contains those that are highly induced only at specific growth and development stages, and the last group contained members with no expression or low expression during all growth and developmental phases. Almost all ARF genes in the sub-group ARFA1 were highly expressed in multiple tissues, and may therefore be involved in the regulation of growth and development.

The differential abundance of homoeologs was analyzed using a previously described framework (Ramírez-González et al., 2018). Balanced expression of homoeolog triads was denoted when transcripts from every gene had a similar abundance. Suppressed and dominant categories were denoted if some transcripts were more abundant than others (Fig. 5B, Fig. S2). Expression data were obtained from a developmental time course of Chinese Spring wheat (Ramírez-González et al., 2018). The percentage of triads in the balanced category was between 60% and 73%, with an average of 68.3% (Table S8), which is consist with the values observed for transcripts from all wheat genes (Ramírez-González et al., 2018). BAR software was used to display the electron fluorescence diagram of TaARFA1a-A expression to better understand TaARF gene expression during growth and development. Our results suggest that some TaARFs may play an important role during wheat growth (Fig. 5C).

Expression analysis of TaARF genes under biotic and abiotic stresses

The TaARFA1 sub-group members that had higher expression during different developmental stages were selected for further analysis. The original RNA-seq data related to abiotic stress (drought and heat) and biotic stress (powdery mildew and stripe rust) from the NCBI database was used for expression profiling of TaARFA1 sub-group members (Table S9). The levels of TaARFA1 genes were up-regulated after 6 h of drought treatments compared with the control, and then down-regulated with prolonged exposure to stress (Fig. 6A). The pattern of expression was opposite under high temperature conditions, with down-regulation followed by up-regulation to near-control levels after 12 h. With drought and heat co-treatment, TaARFA1 genes showed similar expression patterns to under high temperature alone, although the magnitude of down-regulation was greater with multiple stressors. During biological stress, the expression levels of most TaARFA1 sub-group genes fluctuated. Powdery mildew infection caused the expression of TaARFA1 genes (except TaARFA1e-A/B/D) to be up-regulated compared to the control. The expression level of TaARFA1e-A/B/D were down-regulated after infection for 24 h, then up-regulated with prolonged treatment, exceeding that of the control after 72 h (Fig. 6B). The magnitude of changes in expression during stripe rust infection were much smaller than during powdery mildew infection.

Due to the lack of transcriptome data for ABA and osmotic stress in the roots, we selected TaARFA1 genes for qRT-PCR analysis under cold, salt, and ABA stress (Tables S10, S11). The qRT-PCR results revealed that TaARFA1 genes were responsive to all abiotic stress treatments, and their expression patterns varied based on the stress type. Under ABA treatment, the expression levels of most TaARFA1 genes were down-regulated compared with the control, except the expression levels of TaARFA1b-A, TaARFA1c-D, and TaARFA1e-A/B/D, which increased after 6 h of treatment but recovered to the control level after 12 h. The expression levels of TaARFA1c-D, TaARFA1d-A/B/D, and TaARFA1e-A/B were significantly up-regulated during NaCl and cold treatment, while TaARFA1f-A/B/D were down-regulated (Fig. 7).

Cis-acting elements in the promoter of TaARF

The distribution of different cis-acting elements in gene promoters may reflect differences in function and regulation. Identified cis-acting elements were divided into three major categories, namely those relating to growth and development, phytohormone response, and biotic/abiotic stress (Tables S12, S13). Among them, the CAAT-box and TATA-box were the most frequently observed, and related to growth and development. We predicted 3 cis-acting elements related to biotic/abiotic stress (ARE, G-box, and Sp1) and phytohormone response (ABRE, CGTCA motif, and TGACG motif) in TaARFA1 genes (Fig. 8A). The cis-acting elements related to growth and development were enriched in TaARFA1 genes (Fig. 8B).

Gene ontology annotation and protein interaction network of TaARF genes

Gene ontology (GO) annotation is currently one of the most important functional annotation methods. In the present dataset, GO terms related to: (1) biological processes (BP), including regulation of biological process (GO: 0050789), response to stimulus (GO: 0050896), single-organism process (GO: 0044699), cellular process (GO: 0009987), signaling (GO: 0023052), biological regulation (GO: 0065007); (2) cellular components (CC), including cell (GO:0005623), organelle (GO:0043226), cell part (GO:0044464); and (3) molecular function (MF), including binding(GO:0005488), were specifically enriched (Fig. 9A, Table S14). In this study, we also used string website to predict the protein interaction network in which the TaARF genes were involved (Fig. S3). Proteins that share similar functions or participate in the same pathway tend to show interaction networks, so gene clusters or modules were formed in the network of protein interactions. It’s obviously that TaARFA1 and PP2C proteins may have interaction relationship (Fig. 9B). PP2C is a key protein in the ABA signaling pathway; in the absence of ABA, PP2Cs were negatively regulated by repressors that suppress gene transcription (Nguyen, Jung & Cheong, 2019). TaARFA1 genes might interact with PP2C to respond to biotic and abiotic stress.