Genome-wide identification and expression analysis of the VQ gene family in Cucurbita pepo L.

Identification of VQ genes in C. pepo

In present study, two methods were used to identify complete VQ members in C. pepo. Firstly, The Hidden Markov Model (HMM) profiles of the VQ motif (PF05678) was downloaded from PFAM protein family database (http://pfam.xfam.org/) (Punta et al., 2004). Then the model as a probe to perform a BLASTp against the C. pepo genome database using HMMER 3.0. The cut-off E-value was set as 0.1. Secondly, The known VQ motif-containing members of A. thaliana were obtained from The Arabidopsis Information Resource (TAIR) and database (http://arabidopsis.org) were used as queries to conduct a local BLASTp against the protein database of C. pepo. The targeting genes with similarity of E-value less than 1e^–20 were retained for the further analysis. We set e as two values of 0.1 and 1e^–20 for filtering. We set E-value to 0.1 for the initial filtering, and then set E-value to 1e^–20 for in-depth filtering. After merging the two results, the candidate VQ genes were evaluated using the online program PFAM (http://pfam.xfam.org/search), Conserved Domains DataBase Search (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and SMART tool (http://smart.embl-heidelberg.de/) to confirm the presence of VQ motif. The biophysical properties of VQ members, such as, peptide length, molecular weight (MW), isoelectric point (pI) were predicted using the online ExPasy program (http://www.expasy.org/tools/) (Wilkins et al., 1999).

Multiple sequence alignment and phylogenetic analysis

In order to investigate the evolutionary relationships and classification of VQ members in C. pepo, we firstly conducted a multiple sequence alignment based on the full-length amino acid sequence of all CpVQs using the ClustalW tool with auto strategy parameters. According to the alignment results, we constructed a phylogenetic tree using the neighbor-joining (NJ) method with MEGAX software (Tamura et al., 2011). Branch support for the tree topology was estimated by using a bootstrap analysis with 1,000 replicates. The phylogenetic tree was illustrated using the online ITOL (http://itol.embl.de/help.cgi) online tool.

Gene structure analysis and conserved motif detection

The online program Gene Structure Display Server (GSDS) (http://gsds.cbi.pku.edu.cn/) was used to draw a diagrammatic sketch of the intron-exon structure by comparing their coding sequence (CDSs) with their corresponding genomic DNA sequences. The conserved motifs of CpVQs were detected by the online tool MEME (http://meme.ebi.edu.au/). The parameters were set as follows: zero or one occurrence per sequence; maximum number of motifs: 20; and other optional parameters was set default (Bailey et al., 2015).

Chromosome locations, gene duplication and collinearity analysis

For the further investigation of evolution in CpVQs, we firstly constructed the distribution map of VQ members in C. pepo. The physical positions of all VQ genes were retrieved from the GFF3 annotation file using a local script and the schematic diagram of chromosomal location was visualized by TBTOOLS (Chen et al., 2020a). For the identification of tandem duplication events and segmental duplication events in CpVQ family, all CpVQs amino acid sequences were aligned using BLASTp, with an e-value of 1e^–10. As previous research described, two or more genes located on the same chromosome were arranged in a 200 kb distance and shared more than 70% identity as analyzed with BLASTp can be defined as tandem duplication events. The segmental duplication events were identified by using the Multiple collinear scanning toolkits (MCScanX) with default parameters (Wang et al., 2012a). The Circos program was used to draw collinearity maps to exhibit duplicated gene pairs between CpVQ genes as well as the synteny blocks of the orthologous VQ genes between C. pepo and Arabidopsis thaliana, Cucumis melo, Oryza sativa, Zea mays.

Calculating Ka and Ks

The Ka and Ks were calculated to assess the selection history and divergence time of gene families. The values of synonymous (Ks) and nonsynonymous (Ka) substitutions of duplicated VQ genes were estimated by using the KaKs_Calculator 2.0 with the NG method (Wang et al., 2010b). The divergence time (T) was calculated using the formula T = Ks/(2 × 6.1 × 10⁻⁹) × 10⁻⁶ million years ago (MYA) (Kim et al., 2013a, 2013b).

VQ gene expression analysis of C. pepo

The transcriptome expression data of VQ genes under Podosphaera xanthii stress was available from our laboratory. The transcript abundance was represented by fragments per kilobase of exon per million mapped reads (FPKM) values which were calculated based on RNA-Seq reads. The heatmaps showing expression profiles were generated using log10-transformed FPKM values. The results were presented as heatmaps using TBTOOLS software (Chen et al., 2020a).

Plant materials and treatments

Cucurbita pepo (Self-delivered F2) was used in present study. Seeds were planted in a 1:1 (w/w) mixture of soil and sand, cultured in an artificial climatic chamber kept at 30/22 °C with a 18/6 h photoperiod (day/night). Seedlings that germinated after 8 weeks were subjected to different stress conditions: 200 mM NaCl solution, 20% PEG6000 (drought), 4 °C (cold), immersing (waterlogging) and cultured for a total of 24 h, leaves were collected after 0, 12 and 24 h. In order to analyze the expression of VQ genes in different tissues, we collected plant roots, stems, and leaves for RNA preparation. All of the samples were immediately frozen in liquid nitrogen and stored at −80 °C for subsequent total RNA extraction. All samples were tested with three technical replicates and three independent biological replicates.

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

Total RNA was extracted from C. peop using RNAsimple Total RNA Kit (TIANGEN BIOTECH, Beijing, China) according to the manufacturer’s instructions. First-strand cDNA synthesis was accomplished using TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix (Transgen Biotech, Beijing, China). Quantitative Real-time PCR (qRT-PCR) was performed using TB Green Premix Ex Taq II (TliRNaseH Plus) (RR420Q TaKaRa Biotechnology, Beijing, China) on an FTC-3000P system (Funglyn Biotech, Toronto, Canada) with the primers listed in Table S1. The reaction procedure was completed under the following program: 30 s of pre-denaturation at 95 °C, 40 cycles of 5 s at 95 °C, and 30 s at 60 °C, and 4 °C to finish. All samples were tested with three technical replicates and three independent biological replicates. The relative expression level was calculated while using the 2^−∆∆CT method (Livak & Schmittgen, 2001). The actin and CAC genes were used as internal control (Obrero et al., 2011).

Identification and sequence analysis of VQ genes in C. pepo

After manually removing redundant entries through screening (manually delete the genes without VQ motif) and validation of the search results, a total of 44 VQ genes were identified within the whole genome of C. pepo, which were named CpVQ1–CpVQ44 based on their physical locations on the chromosomes. The characteristics concerning each gene, including gene number, length of coding sequence and amino acid sequence, MW and PI of the proteins, are summarized in Table 1. Subsequent sequence analysis of these 44 CpVQs showed that the encoded CpVQ proteins ranged from 73 amino acids (aa) (CpVQ34) to 628 amino acids (CpVQ14) in length (average length: 219 amino acids). Similar to previous studies in Arabidopsis and O. sativa, most VQ proteins contained less than 300 amino acids. The calculated molecular weight and isoelectric points of these proteins varied from 8.34 kDa (CpVQ34) to 67.83 kDa (CpVQ14), and 4.68 (CpVQ40) to 10.67 (CpVQ19), respectively. The length of the coding sequences of this gene family was between 222 and 1,887 bp (average length: 660 bp).

Phylogenetic analysis of CpVQs

To detect the evolutionary relationships and classification of the VQ family in C. pepo, unrooted phylogenetic Neighbor-Joining trees were constructed with the 44 CpVQ proteins and the known VQ protein from Arabidopsis. Through the analysis of the phylogenetic and structural features of the VQ domains, these proteins were divided into eight clades (I–VIII) based on the nomenclature of the Arabidopsis VQs, with two proteins in I and VIII groups, three each in VI and VII, six members in III, eight in II, and group IV and V have the biggest amount of proteins with 10 VQs (Fig. 1). When compared with the other groups, the size of group IV and V were significantly larger. These results were not completely consistent with previous studies in A. thaliana, O. sativa, and Zea mays. In addition, a phylogenic tree was built using the 44 CpVQ protein sequences, which indicated that the groups IV and V in C. pepo had more VQ members than those in Arabidopsis. In addition, it is also possible that the expansion of the two subgroups may have resulted from gene duplications.

Gene structure and conserved motifs analysis

Gene structure can provide more information about the evolutionary relationship in a gene family; thus, the organization of exons and introns in CpVQs was investigated using the GSDS2.0 online tool. Structural analysis of VQs showed that half of the members of group V had introns, such as CpVQ14 that had four introns, which may be related to their long sequences; genes in group VII have longer coding regions; whereas genes of group I have shorter coding regions than other groups. Moreover, 75% of CpVQs were found to be intronless genes, suggesting that many introns were lost during the long evolution. To further analyze the characteristics of CpVQs, the MEME online tool was used to predict the potential conserved motifs. A total of 20 motifs were identified, with lengths ranging between nine and 48 amino acids. Noteworthy, every CpVQ contained 1–6 conserved motifs (Fig. 2). Motif 1 contained the domain present in every VQ member; thus, representing an essential element in the VQ family. In addition, every group was found to have a clearly identifiable motif structure that distinguished it from the other groups. For example, Motif 2 and Motif 5 were prominent in Group IV, Motif 17 was only observed in Group VII, and Motif 6 was only present in Group III. It is worth mentioning that the gene structure and motif organization of CpVQs were similar within the same group, but diverged between groups, which also indicated that the phylogenetic classification was reliable.

Chromosome mapping and duplication events analysis

Two VQ genes (CpVQ1 and CpVQ2) could not be mapped on any chromosome, but the remaining CpVQs were randomly distributed on the chromosomes, except chomosomes 15 and 18 (Fig. 3). Chromosome 19 had the highest number of CpVQs (n = 7), followed by chromosome 10 (n = 6), chromosomes 1, 2, 6, 8, and 14 (n = 3), chromosome 7 (n = 2), and the remaining chromosomes contained one VQ. In addition, most CpVQs were found to be mainly located on the two ends of the chromosomes.

To better understand the evolution mechanism of VQs in plants, the tandem duplication and segmental duplication events of the CpVQ family were evaluated. Based on BLSTAp results, no tandem duplication event was detected. A total of 21 pairs with 28 CpVQs were identified in the whole genome (Fig. 4). Moreover, groups V and IV contained six and five segmental duplication events, respectively. Hence, these results suggested that segmental duplications were the primary driving force responsible for the expansion of the VQ family in C. pepo, which is consistent with previous reports.

Synteny analysis of VQ genes

To further understand the expansion mechanisms of the CpVQ family, we constructed comparative syntenic maps of C. pepo associated with four representative species, including two dicots (Arabidopsis and Cucumis melo) and two monocots (O. sativa and Z. mays). Overall, the CpVQ genes had the most homologous gene pairs with C. melo (n = 27), followed by Arabidopsis (n = 16). However, no homologous gene pairs were observed between C. pepo and O. sativa or Z. mays (Fig. 5). The syntenic analysis results further suggested that the VQ genes are highly conserved in dicotyledonous plants.

We calculated the Ka, Ks and Ka/Ks ratios of all para-homologous gene pairs (Table 2) to explore the evolutionary constraints of repeated CpVQs. The Ka/Ks values of most gene pairs were less than 1.0, which indicated that these gene pairs had undergone purification selection pressure. In addition, 21 pairs of genes varied from 31 to 216 MYA.

Expression pattern of the CpVQ genes in different tissues

In order to explore the possible role of CpVQ gene in the growth and development of C. pepo, we performed qRT-PCR expression analysis in three tissues of roots, stems and leaves. Expression patterns varied among the randomly selected 15 CpVQ genes (Fig. 6). Five CpVQ genes (CaVQ1, 3, 12, 21 and 34) were highly expressed in the root; two CpVQ genes (CpVQ9 and 22) were highly expressed in leaf and root; one CpVQ genes CpVQ40 were highly expressed in stem and root; three CpVQ genes (CpVQ16, 26 and 39) were highly expressed in stem.

Expression pattern of the CpVQs under powdery mildew stress

Several studies have shown that VQs play an important role in plant response to abiotic and biotic stress. Therefore, the expression profiles of CpVQs under powdery mildew stress were evaluated next based on RNA-seq data. Most of the CpVQs were found to be downregulated at early treatment timepoints, such as the members in group I and IV (Fig. 7). In addition, some members were downregulated and then upregulated. A total of 15 CpVQs from different groups were randomly selected for further validation by qRT-PCR (Fig. 8). Accordingly, most of the tested CpVQs were confirmed to be downregulated after powdery mildew infection, with the expression of these genes being correlated with the RNA-seq data at different treatment timepoints. Among these genes, the expression of two genes (CpVQ1, CpVQ39) significantly increased (more than 2-fold) at 12 h, whereas the expression of eight genes (CpVQ3, CpVQ9, CpVQ12, CpVQ16, CpVQ21, CpVQ22, CpVQ36, and CpVQ40) was significantly downregulated. In addition, three genes (CpVQ26, CpVQ33, and CpVQ34) had a high expression value at 24 h after powdery mildew stress.

CpVQs gene expression following abiotic stress by qRT-PCR

To further investigate the role of CpVQs in abiotic stress responses, We randomly selected 15 CpVQ genes from eight groups, and made sure their responses to the drought-, cold-, salt-, and waterlogging-stress.

Under drought treatment (Fig. 9), all CpVQs were upregulated at different treatment timepoints. The expression of six genes (CpVQ1, CpVQ14, CpVQ26, CpVQ33, CpVQ34, and CpVQ40) significantly increased (more than 5-fold) at 12 h, whereas CpVQ21 was significantly downregulated at 12 h but it then increased at 24 h. In addition, five genes (CpVQ3, CpVQ8, CpVQ9, CpVQ16, and CpVQ39) had a high expression value at 24 h after drought stress.

During cold stress (Fig. 10), the expression of eight CpVQs (CpVQ1, CpVQ16, CpVQ22, CpVQ26, CpVQ33, CpVQ34, CpVQ39, and CpVQ40) were significantly upregulated (more than 2-fold) at 24 h. In contrast, four genes (CpVQ3, CpVQ9, CpVQ21 and CpVQ36) were downregulated at 12 h and their levels were further decreased at 24 h. Three members (CpVQ8, CpVQ14, and CpVQ34) were downregulated at early time points but their expression were significantly increased and peaked at 24 h.

Upon salt stress (Fig. 11), six genes (CpVQ16, CpVQ22, CpVQ26, CpVQ33, CpVQ34, and CpVQ40) shared a similar expression trend, rapidly rising at early time points but with subsequently decreased expression. In contrast, CpVQ8 and CpVQ12 showed a trend of downregulation from 0 to 12 h, but then their levels were rapidly increased. In addition, the expression of CpVQ21 showed a trend of downregulation from 0 to 24 h.

For the case of waterlogging treatment (Fig. 12). The four genes (CpVQ3, CpVQ9, CpVQ21 and CpVQ26) upregulated in the first 12 h and then downregulated in the next 24 h. In contrast, CpVQ8, CpVQ12 and CpVQ14 downregulated in the first 12 h and then upregulated in the next 24 h. Interesting, the expression of seven CpVQs (CpVQ1, CpVQ16, CpVQ22, CpVQ33, CpVQ34, CpVQ39, and CpVQ40) were upregulated from 0 to 24 h.

Conservation of the VQ gene family of C. pepo

In higher eukaryotes, genes without introns are very common (Louhichi, Fourati & Rebaï, 2011; Liu et al., 2021). Herein, according to the gene structure, most VQ genes in C. pepo were also found to have no introns. Only 11 VQs had introns, among which CpVQ14 contained four introns. This result is consistent with the lack of introns in 88.2% of Arabidopsis, 90% of Chinese Cabbage, and 92.3% of tomato (Cheng et al., 2012; Zhang et al., 2015; Ding et al., 2019). Subsequently, a phylogenetic tree was constructed based on the VQ proteins of C. pepo and A. thaliana. Through the analysis of the phylogenetic and structural features of the VQ domains, and these proteins were divided into eight clades (I–VIII) based on the nomenclature of the A. thaliana VQs. This analysis revealed that VQs with introns are located in different subfamilies, suggesting that these introns appear relatively independently. Comparative these plants (Cicer arietinum and Medicago truncatula, A. thaliana and lower plants, moss) indicate that most VQ genes have lost introns during the long evolutionary period. Taken together, VQs identified in C. pepo and in other species provide a certain reference value for the evolution of introns in plants. The average length of the CpVQ is 219.3 amino acids, a majority of the C. pepo VQ genes are intronless so they encode relatively small proteins with fewer than 300 amino acid residues, which is highly similar to those reported for Nicotiana tabacum, Arabidopsis, tomato, O. sativa, and other plants (Liu et al., 2020; Cheng et al., 2012; Ding et al., 2019; Kim et al., 2013a). Noteworthy, CpVQ4, CpVQ14, CpVQ23, and CpVQ32 were all found to be close to 400 amino acids.

The genetic structure of the VQ gene and the conserved Motif not only appear in higher plants such as Arabidopsis, O. sativa, Z. mays, among others, but also in lower plants such as moss, which also implies their ancient origins in the evolutionary history and their important role in plant development (Jing & Lin, 2015; Song et al., 2016). Moreover, all 44 VQ proteins of C. pepo were found to harbor conserved VQ domains and contained the same type of Motifs, implying that in the same branch of the evolutionary tree, the closer the genes of VQ protein are, the more similar the gene structure is. Motif 1 was identified as the core Motif that composes the VQ domain, which was included in all pumpkin VQ proteins, thereby endowing the pumpkin VQ protein with specific biological functions. These results were in agreement with those reported for tomato and soybean species (Ding et al., 2019; Wang et al., 2019). In addition, VQ proteins with similar Motif composition were also found to be located in the same sub-branch of the evolutionary tree, as the Motifs composition between different branches was different.

Expansion mechanism of the C. pepo VQ family

Genome replication events play an important role in expanding the size of the genome and diversifying gene functions, as replication events can result in genes with new functions (Rensing, 2014). In higher plants, tandem repeat events and chromosome fragment replication events are the main processes contributing to gene family expansion (Wang, Wang & Paterson, 2012b). Previous studies have shown that diverse WGD events lead to the different sizes of plant genomes (Adams & Wendel, 2005). In this study, 44 VQ genes were identified in C. pepo, which is a slightly high number of genes as compared with other species, such as Arabidopsis with 34 VQ members, or O. sativa and Vitis vinifera with 39 and 18 VQs, respectively (Cheng et al., 2012; Kim et al., 2013a; Wang et al., 2015a). Nevertheless, the number of VQ members in C. pepo was lower than initially predicted while considering its genome size of 261 Mb, as compared with the 125, 389, 486 Mb of Arabidopsis, O. sativa, and V. vinifera, respectively. Therefore, it can be concluded that there is no necessary connection between the size of the genome and the number of family members. Previous studies have shown that chromosome fragment replication is considered to be the main mechanism of VQ gene expansion (Tang et al., 2008); thus, the evolutionary process can explain the number of specific VQ genes in a species, not the size of the genome. The main process for the expansion of the C. pepo VQ family was found to be chromosome fragment duplication events, which is consistent with the results of previous studies. Among the 44 CpVQs, a total of 28 members participated in 21 chromosome fragment duplication events, accounting for 63.6%, and no tandem duplication event was identified. This is consistent with the expansion of the VQ family in Brassica napus, which is driven by chromosome fragment duplication, with tandem repeat events as the second driving force (Zou et al., 2021). Gene duplication can produce gene function redundancy, and these repeated genes can develop different gene expression patterns. We calculated the Ka, Ks and Ka/Ks ratios of all para-homologous gene pairs to explore the evolutionary constraints of repeated CpVQs. The Ka/Ks values of most gene pairs were less than 1.0, which indicated that these gene pairs had undergone purification selection pressure. In the present study, some synlinear genes showed different expression patterns, such as CpVQ3 and CpVQ33, CpVQ9 and CpVQ21, including under drought treatment, suggesting that these homologous genes in C. pepo may have different functions in regulation the normal growth and development of plants.

Expression patterns of VQ members of C. pepo

Previous studies have shown that members of the VQ family play an important role in the entire plant development process and respond to various biotic and abiotic stresses. Molecular genetic evidence suggests that the plant VQ protein may be an important regulator of disease resistance and tolerance (Wang et al., 2015a). In this study, we analyzed the expression level of CpVQ gene in different tissues of C. pepo (Fig. 6). The results showed that most genes were differentially expressed in the tissues we analyzed. It indicates that the CpVQ gene may play an important role in the growth and development of these organs or tissues. In addition, based on transcriptome data, the expression VQs in C. pepo leaves under powdery mildew treatments at different timepoints was evaluated (as observed in Fig. 7). Most VQ members showed a significant decline in expression after 12 h of powdery mildew treatment, however, the expression patterns of the third group is different from the other groups. In the third group, except CpVQ20, all were significantly up-regulated upon powdery mildew infection, and three genes (CpVQ18, 27 and 38) had higher expression at 24 h, which may be because these genes responded late to powdery mildew stress. Studies have shown that SIB1(AtVQ23) and SIB2(AtVQ16) function as activators in plant defense against necrotrophic pathogens (Lai et al., 2011), so the CpVQ genes classified in the third group may have the same function. In contrast, after 24 h of powdery mildew treatment, the expression of most members in other groups decreased significantly. Similar results have been reported in grapevine that VvVQ genes are quickly responsive to powdery mildew stresses (Wang et al., 2015a). Therefore, the members of the VQ family in C. pepo may play a key role in the powdery mildew signaling pathway. Moreover, analysis of VQs expression under different abiotic stresses further demonstrated that the members of the VQ family are also significantly induced under different stresses. In this study, four CpVQ genes (CpVQ3, 9, 16, and 22) significantly upregulated during drought stress (as observed in Fig. 9), which results are similar to the OsVQ genes that 22 OsVQ genes are upregulated under drought stress (Kim et al., 2013a). Under cold stress (as observed in Fig. 10), five CaVQ genes (CaVQ1, 22, 33, 39 and 40) were significantly induced. Similar results have been reported in Chinese cabbage that BrVQ genes are quickly responsive to cold stresses (Zhang et al., 2015). In addition, the VQ genes are also sensitive to salt changes. Two CpVQ genes (CaVQ3 and 9) were upregulated during salt treatment (as observed in Fig. 11). Similar changes occurred in Arabidopsis. AtVQ9 and AtVQ15 were significantly expressed under salt stress. Seven CaVQ genes (CaVQ1, 16, 22, 33, 34, 39 and 40) were upregulated at waterlogging treatment (as observed in Fig. 12). In summary, CpVQ members may be involved in regulating the response of plants to various abiotic stresses and powdery mildew stress, and their response mechanisms may be complex and diverse.

In summary, this study has systematically analyzed the evolutionary relationship, conserved structure, and expression patterns of the members of the VQ family of C. pepo at the whole genome level. The selection of candidate genes can provide reference for future investigations.