Genome-wide identification and expression analysis of the VQ gene family in soybean (Glycine max)

Identification of VQ genes

The Hidden Markov Model (HMM) profiles of the VQ motif PF05678 were downloaded from the Pfam database (Punta et al., 2012). HMM searched VQ motif (PF05678) from the G. max proteins database with the values (e-value) cut-off at 0.1 (Punta et al., 2012). The integrity of the VQ motif was determined using the online program SMART (http://smart.embl-heidelberg.de/) with an e-value < 0.1 (Letunic, Doerks & Bork, 2012). In addition, the three fields (length, molecular weight, and isoelectric point) of each VQ protein were predicted by the online ExPasy program (http://www.expasy.org/tools/) (Rueda et al., 2015).

Phylogenetic analysis

To investigate the phylogenetic relationship of the VQ gene families among A. thaliana, O. sativa, and G. max, AtVQ and OsVQ proteins were downloaded from phytozomes (http://www.phytozome.org) based on the previous studies (Cheng et al., 2012; Li et al., 2014; Goodstein et al., 2012). VQ proteins were aligned using the BioEdit program. A neighbor-joining (NJ) phylogenetic tree was constructed using these proteins through MEGA7.0 software (Tamura et al., 2011). Bootstrapping was performed with 1,000 replications. Genes were classified according to the distance homology with A. thaliana and O. sativa genes (Cheng et al., 2012; Li et al., 2014).

Sequence alignment, motif prediction and gene structure of GmVQ genes

Multiple alignments of the VQ full length proteins were conducted using Jalview software with default parameter settings. The online MEME analysis used to identify the unknown conserved motifs (http://meme.ebi.edu.au/meme/intro.html) using the following parameters: site distribution: zero or one occurrence (of a contributing motif site) per sequence, maximum number of motifs: 20, and optimum motif width ≥6 and ≤200 (Bailey et al., 2015). A gene structure displaying server program (http://gsds.cbi.pku.edu.cn/index.php) was used to show the structure of Glycine max VQ gene.

Gene duplication and collinearity analysis

The physical locations of the GmVQ genes on the soybean chromosomes were mapped by using MG2C website (http://mg2c.iask.in/mg2c_v2.0/). The analysis of synteny among the soybean genomes was conducted locally using a method similar to the one developed for the PGDD (http://chibba.agtec.uga.edu/duplication/) (Krzywinski et al., 2009). First, BLASTP, OrthoMCL software (http://orthomcl.org/orthomcl/about.do#release) and MCScanX software (Wang et al., 2012) were used to search for potential homologous gene pairs (E < 1 e⁻⁵, top five matches) across multiple genomes. Then, these homologous pairs were used as the input for the PGDD database (http://chibba.agtec.uga.edu/duplication/). Ideograms were created using Circos (Krzywinski et al., 2009).

Calculating Ka and Ks

The Ka and Ks were used to assess the selection history and divergence time of gene families (Li, Gojobori & Nei, 1981). The number of synonymous (Ks) and nonsynonymous (Ka) substitutions of duplicated VQ genes was computed by using the KaKs_Calculator 2.0 with the NG method (Xu et al., 2018). The divergence time (T) was calculated using the formula T = Ks∕(2 × 6.1 × 10⁻⁹) × 10⁻⁶ million years ago (MYA) (Kim et al., 2013a; Kim et al., 2013b).

VQ genes expression analysis of soybean

The expression data of VQ genes in different tissues, including seed, pod, SAM, stem, flower, leaf, root, root hair and nodule, is available in Phytozome V12.1 database (https://phytozome.jgi.doe.gov/pz/portal.html). The expression profile for VQ genes was utilized for generating the heatmap and k-means clustering using R 3.2.2 software (Gentleman et al., 2004).

Plant material and treatments

Glycine max (Williams 82) was used in this study. Seeds were planted in a 3:1 (w/w) mixture of soil and sand, germinated, and irrigated with half-strength Hoagland solution once every 2 days. The seedlings were grown in a night temperature of 20 °C and day temperature of 22 °C, relative humidity of 60 %, and a 16/8 h photoperiod (daytime: 05:00–21:00). After 4 weeks, the germinated seedlings were treated with 20% PEG6000 (drought), 250 mM NaCl solution (salt), 4 °C (cold), 100 µM abscisic acid (ABA), 100 µM salicylic acid (SA) solutions. Control and treated seedlings were harvested 1 h, 6 h, 12 h, and 24 h after treatment. All samples were frozen in liquid nitrogen and stored at −80 °C until use.

RNA extraction and Quantitative real-time PCR (qRT-PCR)

Total RNA was extracted from G. max using RNAiso Plus (TaKaRa, Tokyo, Japan) according to manufacturer’s instructions. The cDNA synthesis was carried out with approximately 2 µg RNA using PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa, Tokyo, Japan). Quantitative Real-time PCR (qRT-PCR) was performed using SYBR Premix Ex Taq II (TaKaRa, Tokyo, Japan) on an ABI Prism 7000 sequence detection system (Applied Biosystems, USA) with the primers listed in Table S1. PCR amplification was performed in accordance with SYBR Premix Ex Taq (TaKaRa, Tokyo, Japan) response system. For each sample, three technical replicates were conducted to calculate the averaged Ct values. Relative expression was calculated by the 2^−ΔΔCt method (Livak & Schmittgen, 2001). The actin and GAPDH genes were used as internal control.

Gene Ontology Enrichment

Once the sequences were obtained ran a BLASTX search against the UNIPROT database at a 1e-30 significance level. The matches were extracted and compared to the GO annotation generated against UNIPROT hits located at EBI. The GO annotation of the GmVQ genes by using WEGO 2.0 website (http://wego.genomics.org.cn/).

Analyzed the cis-elements of GmVQ promoters

The cis-elements of GmVQ promoters were analyzed to further understand the GmVQ gene family. We examined the sequences within 1,500 base pairs (bp) upstream of initiation codons (ATG) for promoter analysis and searched for these sequences in the soybean genome. The cis-elements in promoters were subsequently searched using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

Gene interaction network

Protein sequence of GmWRKY transcription factors were obtained from the genome database of soybean, also were mapped to the WRKY proteins of Arabidopsis by BLASTP tool in the TAIR database. Subsequently, the interaction between GmVQs and GmWRKYs were forecasted based on the PAIR website (https://rc.webmail.pair.com/), and their network was drawn in Cytoscape 3.6.1.

Identification of GmVQs

Hidden Markov Model (HMM) of the VQ motif (PF05678) was used to search for putative VQs in soybean proteins database. A total of 75 VQs were identifiedand were named from GmVQ1 to GmVQ75 based on their physical locations on the chromosomes. This is different from the previous study, which 74 GmVQs were identified before the database updated (Wang et al., 2014; Zhou et al., 2016). ExPasy predicted that these 75 VQ proteins have different physical and chemical properties whose amino acid lengths ranged from 89 aa (GmVQ37) to 486 aa (GmVQ18), with an average of 223 aa and most of them were less than 300 aa. The molecular weights of these 75 VQ proteins ranged from 10.03 kDa (GmVQ37) to 52.79 kDa (GmVQ18) and their isoelectric points ranged from 4.29 (GmVQ69) to 10.74 (GmVQ51) (Table 1).

Phylogenetic analysis and multiple alignment of the VQ genes

To explore the phylogenetic relationships among the VQ genes of soybean, A. thaliana and O. sativa, a NJ phylogenetic tree was constructed (Fig. 1). We found that soybean and A. thaliana have a closer relationship than rice. Based on their relationship with AtVQs and OsVQs and the characteristics of GmVQs’ core domain, they were divided into 7 groups, designated Group I-VII (Figs. 1 and 2). For the 75 GmVQ proteins, Group VI contains two VQ proteins; Group V has the biggest amount, with 17 VQ proteins. Groups I, II, III, IV, VII contain 7, 15, 8, 12, 14 members respectively. At the same time, we found 5 types of VQ specificity domain: FxxxVQxLTG (54/75), FxxxVQxFTG (16/75), FxxxVQxVTG (2/75), FxxxVQxLTR (1/75), FxxxVQxLTS (1/75), besides, there is also a GmVQ protein (GmVQ10) has partial domain deletion (Fig. 2). Different types of VQ domains indicate that they might have different biological functions.

Conserved motifs and gene structures of the VQ gene family

We predicted that the 75 GmVQs contained 20 conserved motifs, with the motif length ranged from 11 aa to 50 aa (Fig. S1). Every GmVQ member contains 1-7 conserved motifs (Fig. 3B). All of the proteins, excepted GmVQ22, show motif 1 which contains a specialty VQ domain. Additionally, an unrooted phylogenetic tree was constructed with VQ protein sequences, suggested that the motifs organization of VQ genes were consistent with the phylogenetic tree (Fig. 3A). Group V contains motif 4, Group IV contains motif 2. We found that most groups possess more than two motifs, suggested that every group might have special functions with a highly conserved amino acid residue. Through the VQ gene structures analysis, half of the group VI has introns; genes in group V have longer coding regions, while genes of group I have shorter coding regions than other groups (Fig. 3C). Interestingly, 78.67% (59/75) of GmVQ genes are intronless genes. It is speculated that a large number of introns might be lost in VQ genes during evolution. The phylogenetic tree shows that genes from same branches have similar gene structures, while those from different branches have different gene structures (Fig. 3A).

Chromosome location and gene duplication

We drew a chromosomal location map of GmVQ genes on each chromosome (Fig. 4). GmVQs are distributed on all soybean chromosomes, except chromosome 16, and were densely distributed on chromosome 8 and chromosome 13, containing 13 and 7 members, respectively (Fig. 4). Most of GmVQs are distributed on the two ends of chromosomes.

Segmental or tandem duplicate in many gene families are the main expanding way in plants. To better study the evolution of GmVQ genes, we further explored gene duplication events using the MCScanX software. We found that 52 pairs of genes originated from segmental duplication, and 4 pairs of genes involved in tandem duplication events (Table S2).

Evolution and divergence of the VQ gene family in soybean and Arabidopsis

With the OrthoMCL software, we found 56 paralogous pairs in soybean, 37 orthologous pairs between soybean and Arabidopsis. Some VQ genes have never had any homology genes. All the paralogous and orthologous pairs are listed in Table 2. At the same time, we found that two or more GmVQ genes match to one AtVQ gene, implying that they might promote the expansion of the VQ gene family during evolution. We calculated Ka/Ks ratios of 55 paralogous pairs in soybean (Table 3). Most Ka/Ks ratios are <1, however, the GmVQ54/GmVQ63 and GmVQ65/GmVQ36 pairs are <1. In addition, the genetic differentiation of the 55 gene pairs occurred between 5 and 30 MYA.

Expression analysis of GmVQ genes among various tissues

Sixty-seven GmVQ genes were investigated using available RNA-seq data from nine different tissues (seed, pod, SAM, stem, flower, leaf, root, nodule, and root hair) (Fig. 5). We found that the expression levels of the GmVQs varied significantly in different tissues.Most GmVQ genes were found expressed in more than one detected organ. As shown in Fig. 5, genes in group A are expressed in all analyzed tissues. The expression levels of group B in pod and stem tissues are higher. Genes in group C have specific expression in leaf and root.

Expression patterns of GmVQs under abiotic stress

We randomly selected 25 GmVQ genes from seven groups, and made sure their responses to the plant hormones-, cold-, salt-, and drought-stress (Figs. 6–10). Under ABA treatment, most genes were up-regulated whole treatment period and six genes (GmVQ6/8/31/33/59/71) were obviously down-regulated at some treatment time points (Fig. 6, Table S3). The expression levels of seven genes (GmVQ2/27/40/48/53/68/74) reached the peak at the 6 h treatment time point and four genes (GmVQ9/21/31/71) reached the lowest expression levels at the early treatment time points (0–1 h treatment). With SA treatment, the expression levels of most GmVQs were down-regulated throughout, while GmVQ7 was up-regulated at 1 h, 6 h and 12 h treatment time points (Fig. 7, Table S4). In addition, nine GmVQ genes (GmVQ5/6/8/23/31/68/70/71/74) were down-regulated under all abiotic stress.

With cold treatment, the expression levels of fourteen GmVQ genes (GmVQ2/7/9/28/29/31/33/40/46/48/53/59/68/74) were up-regulated throughout (Fig. 8, Table S5), while the expression levels of three genes (GmVQ27/64/65) were down-regulated and then up-regulated during treatment. Under salt stress, the results were similar to that with cold stress treatment, most genes were up-regulated, eight genes (GmVQ9/23/27/33/65/68/70/71) were down-regulated throughout (Fig. 9, Table S6). On the contrary, under drought (PEG) stress, most genes were down-regulated, only eight genes (GmVQ2/6/7/8/21/29/33/48) were up-regulated during the treatment (Fig. 10, Table S7).

Cis-elements in GmVQ promoters

We found many hormone- and stress- related promoter’s cis-elements in GmVQ genes. Enhancer regions (CAAT-box) and core promoter element are around −30 bp of transcription start (TATA-box). Cis-acting regulatory element (A-box) are the common cis-acting elements in the promoter. Others cis-elements that were found in the 75 GmVQ s can be classified into three groups (Fig. 11). Twelve cis-elements involve in the hormone responsiveness; five cis-elements are stress-related elements: ARE/GC/LTR/MBS/TC; some GmVQ genes containe plant growth and development elements, such as CAT-box/circadian/GCN4/HD-Zip 1/MSA-like/RY-element. In addition, some GmVQ genes containe W-box motif, which is binding site for WRKY transcription factor.

Gene Ontology Enrichment

To further understand the functions of the GmVQs, we performed GO annotation and GO enrichment analyses (Fig. S2 and Table S8). The GO terms included three categories, biological process (BP), molecular function (MF) and cellular component (CC). GO enrichment confirmed that these GmVQs were enriched in the biological process (GO:0008150), regulation of biological process (GO:0050789) and biological regulation (GO:0065007) terms of the BP category. Cellular component (GO:0005575), intracellular (GO:0005622) and cell (GO:0005623) were the most abundant functions in the CC category (Table S8). MF was enriched in molecular function (GO:0003674) and binding (GO:0005488). The GO ebrichment suggested that GmVQs were play curcial roles in regulated of biological process.

Gene interaction network analysis

Based on the PAIR tool, we found the functions and their interactions of the GmVQs and GmWRKYs. As shown in Fig. 12A, 3 GmWRKYs are supposed to interact with GmVQ proteins, included GmWRKY115, GmWRKY149 and GmWRKY156, all of them belong to WRKY’s IIc group. In the Fig. 12B, we found that GmWRKYs and AtWRKYs are quite similar in their core domains, indicated that they might have same functions, such as interacted with VQ proteins.