Genome-wide identification, evolution, and expression of the SNARE gene family in wheat resistance to powdery mildew

Identification of TaSNARE genes

The wheat genomes and annotations used were from the newest IWGSC (International Wheat Genome Sequencing Consortium) v1.0 (https://wheat-urgi.versailles.inra.fr/Seq-Repository). The hidden Markov models (HMMs) of the SNARE (PF05739), Syntaxin (PF00804), longin (PF13774), Synaptobrevin (PF00957), SEC20 (PF03908), V-SNARE-C (PF12352), V-SNARE (PF05008) and USE1 (PF09753) motifs were downloaded from the Pfam database (http://pfam.sanger.ac.uk/). The wheat SNARE protein sequences were analyzed with HMMER 3.0 (http://hmmer.janelia.org/) as the query and the default parameters (E < 0.01). All presumptive SNARE genes were retained and confirmed using the Pfam database and the NCBI conserved domain database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). The molecular weights (MW) and protein isoelectric point (pI) of the TaSNARE genes were obtained using the tools from the ExPASy website (https://www.expasy.org/).

Multiple alignments of SNARE proteins were performed using ClustalW (Larkin et al., 2007) in MEGA 7.0 (http://www.megasoftware.net/). Phylogenetic analyses were performed using the NJ (neighbor-joining) method in MEGA 7.0 (Kumar, Stecher & Tamura, 2016) with 1,000 bootstrap resampling, the Jones-Taylor-Thornton (JJT) model (Jones, Taylor & Thornton, 1992), and the pairwise deletion option. Gene Ontology (GO) enrichment analysis of SNAREs was implemented using the clusterProfier R package (R version=4.0.3).

Exon/intron structure analysis, conserved motif identification and cis-acting elements analysis

The gene structure provides important information, including disaggregated and evolutionary relationships among gene families. The SNARE genomic sequences and CDS sequences extracted from the plant database were compared in the gene structure display server program to determine the exon/intron organization of SNARE genes. The default parameters were used in the Multiple Em for Motif Elicitation (MEME) (http://meme-suite.org/) program for the identification of conserved protein motifs and a maximum number of 15 motifs. Promoter of the SNARE gene was used to analyze the cis-acting elements as previously described by Sun et al. (2017).

Chromosomal locations and gene collinearity analysis

The physical chromosome locations of all SNAREs were obtained from the gff3 files of wheat databases. TBtools (https://github.com/CJ-Chen/TBtools) software was adopted to visually map the chromosomal location. Gene duplication events were analyzed using the Multiple Collinearity Scan toolkit (MCScanX: http://chibba.pgml.uga.edu/mcscan2/). To exhibit segmentally duplicated pairs and orthologous pairs of SNARE genes, we used Dual Systeny Plotter software (https://github.com/CJ-Chen/TBtools) to drawn collinearity maps.

Fungus and wheat materials

The wheat–Ae. geniculata disomic addition line NA0973-5-4-1-2-9-1 (CS-SY159 DA 7M^g, (CS)/Ae. geniculata SY159//CS)) was used (Wang et al., 2016). ‘Shaanyou 225’ was the powdery mildew susceptible control variety. The wheat–Ae. geniculata disomic addition line TA7661 (CS-AEGEN DA 7M^g) was kindly provided by Dr. Friebe BR and Dr. Jon Raupp of the Department of Plant Pathology, Throckmorton Plant Sciences Center, Kansas State University, Manhattan, KS, USA. The powdery mildew isolates E09 was maintained on the susceptible wheat ‘Shaanyou 225’. All plants were cultured in a growth chamber with soil at 18 °C under a 16 h light/8 h dark photoperiod. The 14-day-old seedlings were inoculated with powdery mildew conidia from ‘Shaanyou 225’ seedlings infected 10 days previously using the dusting method. Wheat leaves were collected at 0 h, 6 h, 12 h, 24 h, and 48 h after powdery mildew infection, and quickly put into cryopreservation tubes and stored in liquid nitrogen. The leaves were used for the next step of RNA extraction and q-PCR experiments. The method by Wang et al., (2020) was used for CS (Chinese Spring), 7M CH (NA0973-5-4-1-2-9-1), 7M US (TA7661) and ‘Shaanyou 225’ to identify powdery mildew.

RNA-seq expression analysis of SNARE genes

To further understand the function of the SNARE gene, we investigated the reported RNA-seq data, including the developmental timecourse in five tissues (Choulet et al., 2014), grain layers (Pearce et al., 2015), grain layer developmental timecourse (Pfeifer et al., 2014), senescing leaves timecourse (Pearce et al., 2014), photomorphogenesis of DV92 and G3116 (Fox et al., 2014), and drought and heat effects (Liu et al., 2015). The data were analyzed using MeV (Multi Experiment Viewer) software. Data obtained from the RNA-seq expression atlas were normalized based on the mean expression value of each gene in all tissues/organs analyzed and clustered by the hierarchical clustering method.

The developmental time course in five tissues includes all of the wheat stage, as follows (Zadoks, Chang & Konzak, 1974): seeding (first leaf through coleoptile, Zadoks Scale 10, Z10), three leaves (three leaves unfolded, Z13), three tillers (Main shoot and 3 tillers, Z23), spike at one cm (pseudostem erection, Z30), two nodes (2nd detectable node, Z32), meiosis (flag leaf ligule and collar visible, Z39), anthesis (1/2 of flowering complete, Z65), 2 days after anthesis (DAA) (Kernel (caryopsis) watery ripe, Z71), 14 DAA(medium Milk, Z75), and 30 DAA (soft dough, Z85). The grain layers contained three parts at 12 days post-anthesis (DPA): the outer pericarp, inner pericarp, and endosperm. The grain layer developmental timecourse included the following seven-stages: 10 DPA whole endosperm, 20 DPA whole endosperm, 20 DPA starchy endosperm, 20 DPA transfer cells, 20 DPA aleurone, 30 DPA starchy endosperm, and 30 DPA aleurone plus endosperm. The senescing leaves timecourse contains three stages: heading date (HD), 12 DAA and 22 DAA. The photomorphogenesis of the wild winter wheat T. monococcum ssp. aegilopoides (accession G3116) and the domesticated spring wheat T. monococcum ssp. monococcum (accession DV92) was investigated. Drought and heat effect examinations included seven treatments,as follows: control, drought 1 h, drought 6 h, heat 1 h, heat 6 h, and drought plus heat 1 h, drought plus heat 6 h. Powdery mildew pathogen stress: included non-innoculation, powdery 24 h, powdery 48 h and powdery 72 h.

RNA extraction and real-time quantitative PCR

The total RNA was extracted from samples of fungal inoculated leaves using the optimized extraction procedure described by Zhang et al. (2014).

The SYBR Green Premix Ex Taq™ II quantitative PCR system (Takara, Dalian) was used for qPCR analysis. All experiments involving q-PCR were performed on a Q7 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The actin gene (GenBank: aK458277.1) was used as the reference gene. The PCR reaction and program were modified according to the manual. The PCR reaction (a total reaction volume of 10 µL) comprised 5 µL 2 × SYBR Green PCR Master Mix, 3 µL of the cDNA product, 1 µL of primer mix, and 1 µL of DNase/RNase-free water. The quantitative PCR thermal cycler program included 95 °C for 10 s, followed by 40 cycles at 95 °C for 5 s and 60 °C for 31 s. All primers for q-PCR were synthesized by the same company (AoKe, Yangling) (Table S4).

Identification of the SNARE protein in wheat

To identify SNARE proteins in wheat, the HMMER profile was implemented to identify the wheat genomes. The results showed that 173 hypothetical TaSNARE genes were characterized from wheat databases (Table S1). Qa, Qb, Qc, Qb+Qc and R SNARE proteins comprised 48(27.7%), 37(21.4%), 39(22.5%), 13(7.5%) and 36(20.8%) respectively. Among all 21 subfamilies, SYP1 contained a maximum of 33 proteins, and VAMP72 had the second most, at 15 proteins (Table S1). The encoded proteins comprised between 121 and 466 amino acid residues, the PIs ranged from 4.72 to 9.65, and the molecular weights were distributed from 13,687.37 to 51665.97 Da (Table S1).

All the sequences were divided into 64 groups in wheat (Table S1). Among these groups, 38 groups representing 114 genes contained three genes from each of the different subgenomes that were regarded as orthologous copies of a single SNARE gene named triplet. Five groups contained different homoeologous genes that were from the same homoeologous group (e.g., TaSYP43-4AL, TaSYP43-7AS, and TaSYP43-7DS). Eight groups contained two genes (e.g., TaSYP131-2BS and TaSYP131-2DS), and the remaining 8 groups consisted of only one gene (e.g., TaGOS12-6BS). Five groups had four genes, among which four groups had tandemly repeated genes (e.g., TaSNAP1-2A1, TaSNAP1-2A2, TaSNAP1-2B, and TaSNAP1-2D). We found that the Go term “vesicle-mediated transport” was most significantly enriched in the SNARE proteins (Table S2).

Chromosomal locations and gene collinearity analysis of SNARE gene family members in wheat

Most collinear gene pairs occur within the same chromosome group (Fig. 1). The chromosomal distribution of the SNARE gene family of T. aestivum was analyzed. Figure 2 revealed the chromosomal location of 173 SNARE genes. All 21 wheat chromosomes have several SNARE gene family members: the wheat 1 to 7 homoeologous groups had 14 (1A = 5,1B = 4,1D = 5), 22 (2A = 7, 2B = 8, 2D = 7), 31 (3A = 12, 3B = 9, 3D = 10), 25 (4A = 9, 4B = 9 4D = 7), 23 (5A = 8, 5B=9, 5D=6), 26 (6A = 9, 6B = 9, 6D = 8) and 27 (7A = 11,7B = 8,7D = 9) SNARE genes, and 4 had no chromosomal location. In addition to homoeologous group 1, the SNARE genes were evenly distributed in the wheat genome, and the number of genes on each chromosome is similar. The most striking result to emerge from Fig. 2 was that the triplets, which were from different subgenomes, were similar in terms of their relative position on chromosomes.

Phylogenetic, motif and structural analysis of the SNARE family genes

To further analyze the phylogeny, motif and structure of TaSNAREs, we selected one protein (genome group A chromosomal priority selection) from each of the 64 groups and obtained 64 SNAREs. The results showed that these proteins were primarily divided into 5 clades (Fig. 3). Most clades had three homoeologous proteins in the same branch, and these three homoeologous proteins were from three chromosomes in the same homoeologous group.

It is apparent from Fig. 3 that SNAREs in different subfamilies had different motifs. Qa had motifs 1, 2, 5, 8, 11 and 13. Qb had motifs 6, 7, 10, 12 and 13. Qc had motifs 5, 9 and 13. Qb+Qc had motifs 9 and 12. R had motifs 3, 4, 13 and 14. The results showed that motifs 6, 8, 10, 11, 13 and 15 were not predicted as being present in these SNAREs. Motif 1 was the SNARE domain; motifs 2 and 5 were syntaxin domains; motif 3 was the synaptobrevin domain, motifs 4 and 14 were longin domains; motif 7 was the SEC20 domain; and motif 12 was the V-SNARE-C domain. Qa, Qb, Qc, and R had motif 13, which was located in the C-terminus and associated SNAREs with lipid bilayers, and this motif was named the transmembrane (TM) domain (Lipka, Kwon & Panstruga, 2007).

Cis-acting elements of TaSNARE genes

Further analysis of the 2,000-bp promoter upstream of the 5′ end of the TaSNARE gene was performed. This promoter contains 9 types of resistance-related cis-acting elements (Table S3), including W-box (Cis-I), germ-related (Cis-II), MYB (Cis-III), SA-responsive (Cis-IV), Eth-responsive (Cis-V), EIRE (Cis-VI), G-box (Cis-VII), H-box (Cis-VIII) and IAA-responsive (Cis-IX) elements.

As shown in Fig. 4, Cis-I to Cis-IX were represented by 2230, 5054, 1647, 882, 170, 225, 152, 118 and 309 elements, respectively, in all 173 SNARE gene promoters. Cis-I, Cis-II and Cis-III made up 82.92% of all disease-related elements. Among the Cis-I elements, TaUSE12-7A was the most abundant (36). Among the Cis-II elements, TaSEC222-5B was the most abundant (87). Among the Cis-III elements, TaSEC222-5A was the most abundant (19).

In one triplet, for the promoters of the resistance-related elements, the numbers were similar. However, there were exceptions, as follows: TaSFT11-2A/B/D had 5∕21∕16 Cis-I elements; TaNPSN12-4A/B/D had 71∕17∕21 Cis-II elements, and TaSYP222-6A/B/D had 3∕11∕1 Cis-IV elements.

Expression analysis of TaSNARE genes by RNA-seq

To further understand the functions of the SNARE genes, we extracted gene expression information for 54 genes from six published RNA-seq databases (Fig. 5).

As shown in Fig. 5, in the growth period of wheat, the SNARE gene is expressed in roots, stems, leaves, seeds and spikes, showing low levels in seeds and leaves and high levels in roots, stems, and spikes. In the seeds of Z75, the expression levels of most genes (45) were very low, and in Z71-Z75-Z85, a high-low-high expression pattern was observed. Many genes (36) were most highly expressed in the 20 DPA aleurone layer during seed development. Most SNARE genes (49) are better expressed under light conditions than under dark conditions. Among these genes, SYP122-6A showed higher expression under light than in the dark in DV92, but in G3116, the opposite trend was observed. Compared with the control, the expression of 22 genes was upregulated 6 h after stress (drought 6 h, heat 6 h or drought plus heat 6h), and the expression patterns of 9 genes showed the opposite trend. In the process of leaf senescence, 42 genes showed the highest expression at 12 DAA. More than half of the genes (33) had the following expression distribution pattern in the grain layers: outer pericarp >inner pericarp >endosperm. In the powdery mildew pathogen stress, 3 genes were down-regulated by more than 0.5 times in 24 h, and one gene was up-regulated by more than 1 times. In total, 17 genes were up-regulated by more than 0.5 times, of which 8 genes were up-regulated by more than double. There were 6 genes down-regulated by more than 0.5 times at 48 h, and one gene was down-regulated by more than 1 times. There were 12 genes up-regulated by more than 0.5 times, of which five genes were up-regulated by more than 1 time. At 72 h, 5 genes were down-regulated by 0.5 times, of which 2 genes were down-regulated by more than 1 time. There were 13 genes up-regulated by more than 0.5 times. There were 6 genes that were continuously up-regulated between 24 h–72 h. TaSYP135 continued to be down-regulated between 24 h–72 h. These genes with the most dramatic changes in expression may have played a role in responding to powdery mildew infection.

Expression patterns of TaSNARE genes under powdery mildew treatment

We selected one gene from each of the 21 classes, and we obtained 21 TaSNAREgenes (TaYKT6 had no signal) specific to the designed primers (Table S3). As shown in Fig. 6, the expression patterns of different SNARE genes in the same sample and subfamily were similar. Most of the TaSNARE genes had similar expression pattern in 7M US and CS, but in 7M CH, a different expression pattern was observed. A majority of the TaSNARE genes in 7M CH had high expression levels at 6 h. TaSYP4, TaSYP8, TaMEMB and TaSEC22 had high expression levels at 6 h in ‘Shaanyou 225’ but not in the other wheat.

In the Qa subfamily, the expression of all genes changed little at each time point in CS. SYP121, SYP221, and SYP3 were upregulated in the 7M CH 6 h sample but downregulated in ‘Shaanyou 225’. SYP4 and SYP8 were upregulated in the ‘Shaanyou 225’ 6 h sample but not in the 7M CH sample. QaSNAREexpression was similar between 7M US and ‘Shaanyou 225’. This may indicate that the up-regulated expression of SYP4 and TaSYP8 played a negative role in the wheat response to powdery mildew infection.

In the Qb subfamily, GOS12 expression was upregulated at 6 h and then downregulated in four wheat varieties. MEMB expression was upregulated in ‘Shaanyou 225’ at 6 h, 24 h, and 48 h but not in CS. There was no significant difference in the expression at different time points in 7M CH and 7M US. The VTI12 expression patterns were similar to those of MEMB in ‘Shaanyou 225’ and 7M CH. NPSN11 expression in four wheat varieties were similar, showing downregulation at 6 h–48 h, except for the upregulation at 24 h–48 h observed in 7M CH. SEC203 was downregulated at 6h-48 h in 7M US and CS. This gene was downregulated at 12 h–24 h and upregulated at 48 h in ‘Shaanyou 225’. In 7M CH, SEC203 was downregulated at 24 h, and the other genes were upregulated. This shows that the NPSN11 gene has little effect in the early stage of wheat’s response to powdery mildew infection, while GOS12 and SEC203 play a certain role.

In the Qc subfamily, all QcSNAREs in the same variety were similar. In CS and 7M US, most genes were downregulated at 6–24 h. In 7M CH, certain genes were downregulated at 24 h, while the others were upregulated. In ‘Shaanyou 225’, the genes were downregulated at 12 h–24 h and upregulated at 6 h and 48 h. This means that after powdery mildew infects wheat for 48 h, BET1, SFT1, USE12, and SYP5 were resistant to 7M CH up-regulation and susceptible to 7M US down-regulation, which implies that the up-regulation of these genes has a positive effect in these two types of wheat.

In the Qb+Qc subfamily, SNAP1 was upregulated at 6 h in all varieties and downregulated at 12 h–48 h, except in ‘Shaanyou 225’ at 48 h, in which SNAP1 was upregulated. This may indicate that the SNAP1 gene may play a similar role in the four types of wheat.

In the R subfamily, in CS, VAMP712 was upregulated at 6 h and at the other time points, it showed no change. In ‘Shaanyou 225’, VAMP712 was downregulated at 6 h and 24 h, and at the other time points there was no change. In 7M CH, VAMP712 was upregulated at 6 h, 12 h, and 48 h and showed no change at 24 h. In 7M US, VAMP712was upregulated at 24 h and downregulated at the other time points. No signal for VAMP723 was detected in ‘Shaanyou 225’. In 7M US, VAMP723 was upregulated at 24 h and downregulated at other time points. In 7M CH, VAMP723 was upregulated at 12 h and downregulated at 24 h and 48 h. In CS, VAMP723was upregulated at 6 h and downregulated at 48 h. SEC222was upregulated at every time point in ‘Shaanyou 225’, and the other varieties showed no significant difference. This implies that the up-regulated expression of SEC222 plays a negative role in the powdery mildew infection of ‘Shaanyou 225’.

With the evolution of plants, wheat SNARE genes are constantly changing. SNARE genes have different expression patterns after being infected by powdery mildew, suggesting that these genes may be involved in the biological stress response to powdery mildew in different aspects.