Genome-wide identification and functional characterization of CDPK gene family reveal their involvement in response to drought stress in Gossypium barbadense

Identification of the CDPK family genes in Gossypium barbadense

The G. barbadense genome data comes from https://cottonfgd.org/ (Zhu et al., 2017), The Arabidopsis CDPK protein sequences of family members were downloaded from the Arabidopsis genome database TAIR website (https://www.arabidopsis.org/). Using Pfam online database (http://pfam.xfam.org/) to download the CDPK gene seed file EF-hand (Pfam ID: PF00036) and serine/threonine protein kinase (Pfam ID: PF07714). The hmmsearch program of HMMER3.0 was used to identify protein sequences containing CDPKs conserved domains. All CDPK protein sequences were extracted NCBI CDD and Pfam online database were used to verify the conservative domain. Finally, all members of CDPK gene family were obtained and named according to their different positions on chromosomes in G. barbadense.

Molecular structure and physicochemical properties

The Protparam tool (https://web.expasy.org/protparam/) was used to predict the molecular weight, isoelectric point and other relevant information of the CDPK protein. Using the online tool myristoylation (https://web.expasy.org/myristoylator/) predict N-myristoylation, and palmitoylation prediction tools is CSS-plam. Subcellular localization of the CDPK proteins were predicted using an online tool wolf POSRT (https://www.genscript.com/psort.html).

Phylogenetic tree, gene structure, conserved motif, and promoter region analyses

The ClustalX software was used to perform multiple sequence alignments on the CDPK proteins of G. barbadense and Arabidopsis, and the phylogenetic tree was constructed by MEGAX software: the Neighbor-Joining algorithm was used, and the Bootstrap repeat value was 1,000 times (Kumar et al., 2018). The online GSDS tool (http://gsds.cbi.pku.edu.cn) (Hu et al., 2015) was used to visualize the exon-intron structures of the CDPK genes. Using MEME (http://meme-suite.org/tools/meme) (Bailey et al., 2009) online software to analyze the conserved motifs of G. barbadense CDPK, the parameters were set as the maximum number of motifs found was 10, and the longest motif length was 100 nt (Nucleotide). The 1,500 bp DNA sequence upstream was extracted, and the PlantCARE (Lescot et al., 2002) database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) was used to predict the possible cis-acting elements.

Chromosome location and gene replication

The information of CDPK gene family chromosome physical location was extracted from G. barbadense genome annotation file, and the chromosome physical distribution map of the genes was drawn with MapInspect. The definition of gene duplication is as follows: (1) the sequence length is greater than 80% of the gene; (2) the identity of alignment region is greater than 80%; (3) only one repeat event is calculated for closely linked genes. According to the chromosomal location of CDPK gene, tandem duplication and segment duplication can be identified. The KaKs_Calculator2.0 of Tbtools (Chen et al., 2020) was used to calculate the non-synonymous mutation rate (Ka) and the synonymous mutation rate (Ks) of the replication gene pair in G. barbadense.

Expression profile analysis, RNA extraction and real-time PCR

The transcriptome data of GbCDPK genes in the tissues and the abiotic stresses were from Sequence read archive, SRA, and the accession number was PRJNA274882. The RPKM value (the reads per kilobase of transcript per million mapped reads) represents the abundance of gene expression at the transcription level. Using Tbtools (Chen et al., 2020) map the gene expression profiles of different tissues and different abiotic stresses. Total RNA is extracted with the plant total RNA extraction kit (TIANGEN, Beijing, China). RNA reverse transcription adopts HiScript II QRT SuperMix for qPCR (R222-01), a product of Nanjing Novazin. The cDNA was used as a template to perform RT-PCR and qPCR experiments. The internal reference gene is the ubiquitinated protein Ubiquitin7 (UBQ7), and the RT-PCR reaction uses 1.0% agarose gel to detect the amplification results. The quantitative instrument uses Roche products, and the calculated Ct value is used to calculate the relative expression of the gene. The calculation method is 2^−ΔΔCt (Pfaffl, 2001). Each gene amplification result contains at least 3 technical replicates and 3 biological replicates. The amplification system and procedure are shown in Appendix A and Appendix B.

Experiment treatment

The G. barbadense was grown for 4 weeks in a high humidity and high nutrient environment (nutritional soil: vermiculite: perlite = 3:1:1) at 28 °C, and then treated with abiotic stress to well growing cotton seedlings. Cold stress and high temperature stress are moved plants to incubators at 4 °C and 37 °C, respectively. For salt stress, seedlings are irrigated with 200 mM NaCl. Seedlings are irrigated with 20% PEG6000 for drought stress. One hour after the start of the treatment, the second true leaf of the treated cotton seedling was harvested and stored in a refrigerator at −80 °C. Each treatment included 3 seedlings, and each treatment had three experimental replicates.

VIGS experiment in G. barbadense

The silenced fragments of each CDPK gene were amplified from cotton and inserted into VIGS vector (TRV-00) to generate CDPK silent constructs. TRV-00 without the silent segment was used as control. Consequently, the vector was transformed into the GV3101. GV3101 contained the same amount of TRV-00 and TRV-CDPK, and infiltrated into the cotyledons of 14-days-old cotton seedlings by syringe infiltration method. Two weeks after inoculation, the silenced cotton seedlings were used for gene expression analysis and identification PEG treatment.

Methods of physiological index

RWC was measured by using the formula: $\mathrm{RWC}\left(\text{\%}\right)=\frac{\text{FW}-\text{DW}}{\text{TW}-\text{DW}}\times 100$. Where FW, stands for fresh weight TW for turgid weight and DW for dry weight (Loutfy et al., 2012). G. barbadense leaves were soaked in 30 mL deionized water at 25 °C and then shaken for 30 min. Heated the solution to 100 °C for 30 min and then at room temperature. Measure REL with conductivity meter (Yang et al., 2012). The content of MDA, Pro, soluble sugar measured using the fresh sample. The amount of MDA was measured by the thiobarbituric acid method (Heath & Packer, 1968). The content of Pro was determined by the Sulfosalicylic acid indandione method (Bates, Waldren & Teare, 1973). Soluble sugar content was measured as described by the anthrone method and determined at 620 nm by a spectrophotometer (Fales, 1951). Water loss rate (%) = $\frac{W1-W2}{W1}\times 100$. W1 stands for weight of detached blade before water loss. W2 weight of detached blade after 1, 2, 3, 4 and 5 h of water loss (Czyczyło-Mysza et al., 2018). All physiological results were repeated three times and a one-way analysis of variance (ANOVA) was performed with Prism to detect significant differences between both treatments (P < 0.05, P < 0.01).

Identification and proteins characteristic of CDPKs in Gossypium barbadense

CDPK protein sequences of Arabidopsis were provided as direct queries against the G. barbadense genome to perform BLASTP and hidden Markov Model (HMM) to identify the GbCDPK gene family members. After NCBI-CDD and Pfam software database verification, 84 CDPK genes were identified in G. barbadense genome (Table S1).

The biochemical characteristics of the CDPK gene, such as amino acid sequence, molecular weight and isoelectric point were analyzed by the ProtParamtool online software. The results show that the CDPK genes are conservative. The longest CDPK protein encodes 648 amino acids, and the shortest only encodes 155 amino acids. The molecular weight of GbCDPKs varies 17.99 KDa to 71.854 KDa, and the theoretical isoelectric point ranges from 4.313 to 9.481. Multiple sequence alignment analysis indicated that all of GbCDPKs have four core domains. In addition, all the GbCDPKs have EF-hands motifs, which allows binding of Ca²⁺. It reported that 71 GbCDPK proteins contain four EF-hands, nine proteins (GbCDPK25, GbCDPK26, GbCDPK39, GbCDPK49, GbCDPK53, GbCDPK67, GbCDPK68, GbCDPK77, GbCDPK81) have three EF-hands, and remaining four proteins (GbCDPK22, GbCDPK38, GbCDPK61, GbCDPK64) only contains 2 EF-hands. Our bioinformatic analysis revealed that 40 members harbored potential N-myristoylation sites, and 74 GbCDPK proteins contained S-palmitoylation sites (Table S1). Using TMHMMServerv2.0 prediction, it was found that none of the 84 CDPK proteins had a transmembrane structure. Using Wolf POSRT tool to predict subcellular localization, it was found that 22 CDPK proteins are in the chloroplast, 33 CDPK proteins are in the cytoplasm, eight are in the nucleus, and the others are in the mitochondria, extracellular, vacuole, endoplasmic reticulum, and peroxisome.

Phylogenetic analysis of GbCDPKs

To assess the phylogenetic relationship of CDPK protein in G. barbadense. A total of 118 CDPK proteins from G. barbadense and Arabidopsis were used to construct the phylogenetic tree by adjacency method (Fig. 1). According to the results of phylogenetic tree clustering, 84 CDPK genes can be divided into four subfamilies. Referring to the research results of Arabidopsis, four subfamilies were named as class I, class II, class III, and class IV. Further analysis of the clustering results of CDPK members revealed that the gene structures in the homologous gene pair were similar. The number of subgroup II proteins was the largest, including 30 G. barbadense and 10 Arabidopsis CDPK proteins, accounting for 35.71% and 29.41% of the total, respectively. The first subgroup includes 27 CDPK proteins from G. barbadense and 10 CDPK proteins from Arabidopsis. The third subgroup includes 19 CDPK proteins from G. barbadense and eight CDPK proteins from Arabidopsis. The fourth subgroup has the least number, including only eight CDPK proteins from G. barbadense and 3 CDPK proteins from Arabidopsis, indicating that it’s far from each other.

Chromosomal distribution and gene duplication

To investigate the chromosomal locations, chromosomal maps were drawn for the 84 CDPK genes from G. barbadense by MapInspect software. All 84 GbCDPK genes are mapped on 13 pairs of chromosomes (Fig. 2). Furthermore, cotton genome includes A-subgroup and D-subgroup that all contain 42 GbCDPKs, respectively, and the duplication events may illuminate the mechanism about the expansion of GbCDPK gene family. There are six genes on chromosome 4 and 13 of A-subgroup, five genes on chromosome 2 of A-subgroup and 5 and 13 of D-subgroup, respectively. In contrast, there is only one gene on chromosome 3 and 8 of A-subgroup and chromosome 8 of D-subgroup, and there are 2–4 genes on other chromosomes. These results indicated that CDPK genes were widely distributed in G. barbadense. In order to speculate the possible relationship between GbCDPK genes and potential gene replication in the G. barbadense genomes, we analyzed the occurrence of tandem duplication and large-scale segmental replication during the evolution of the gene family (Table S2). Among 84 genes, 41 pairs of homologous genes were found in phylogenetic tree (Fig. 1). There were 40 pairs of segmental replications and one pair of tandem replications (GbCDPK10/GbCDPK11).

To study the evolutionary replication relationship of genes in G. barbadense on chromosome segments, MCscanx software was used to detect the replication genes. 12, 16 and 60 pairs of whole genome replicates were detected among A, D and AD subgroups in G. barbadense, respectively (Fig. 2). KaKs_Calculator was used to perform Ka/Ks analysis replication gene pairs. It was found that the Ka/Ks of all gene pairs were less than 1, indicating that the CDPK gene of G. barbadense might have undergone strict purification and selection during the evolution process, implying that the replication gene is evolutionarily conserved, structurally stable, and might have consistent function.

Analysis of gene structure and conserved motifs

The CDPK genes were visualized for the gene structure and conserved motifs. The results showed that they have similar gene structure distribution patterns in the same subfamily. Most CDPKs of the first and the third subfamily have seven or eight exons, the CDPKs of the second subfamily have eight or nine exons, and the fourth subfamily has 12 exons. Most CDPKs of the family have 12 exons (Fig. 3). Generally, the structural conservation and similarity of the same type of genes are closely related to their evolutionary relationship, but there are some exceptions. GbCDPK38 has only five exons and GbCDPK61 has only three exons. The longest exon of most GbCDPKs is the first exon. However, there are still a small number of genes in the family members that have evolved significantly different structures. GbCDPK39 and GbCDPK81 are homologous genes from different sub-genomes, GbCDPK39 has a typical gene structure of this group, while the first exon of GbCDPK81 is relatively short. GbCDPK10 and GbCDPK11 are homologous genes of the same sub-genomes, but their gene structure is also different. GbCDPK10 has typical gene structure characteristics, but the first exon of GbCDPK11 is also shorter than that of GbCDPK10. This structural difference is likely to have formed randomly during evolution. The changed gene structure may lead to changes in gene splicing sites, thereby generating new transcripts or non-functional gene products, and ultimately increasing the number of members of the cotton GbCDPKs gene family.

The MEME online software was used to conduct conservative analysis on the amino acid sequence of CDPK, and 10 conservative motifs were found (Fig. 4). These 10 conserved motifs exist in subgroups I, II, III and IV, but there are four genes (GbCDPK11, GbCDPK14, GbCDPK38, GbCDPK61) whose conservative motifs are distributed differently in the first subgroup. GbCDPK11, GbCDPK14 contains the same conserved motifs, without Motif 2, Motif 3, Motif 5. While GbCDPK38 has only three conserved motifs (Motif 8, Motif 9, Motif 10), GbCDPK61 has two conserved motifs (Motif 8, Motif 10). However, in the same subfamily, CDPK has a similar number, type, and spatial distribution of motifs, demonstrating that the functions of CDPK genes in the same subfamily are similar.

Analysis of promoter cis-elements

To explore the regulatory mechanisms of GbCDPK genes, the 1.5kb sequence upstream of the start condon of the GbCDPKs was selected to analyze the constitutions of cis-regulatory elements. It was found that all CDPK genes contained CAAT-box conserved elements and TATA-box. In addition to these conserved elements, there are four types of cis-regulatory elements in the promoter of CDPK gene: (1) Light regulatory elements, including Box 4, G-Box, TCCC-motif and AE-box, etc. (2) Plant growth and development regulatory elements, including CAT-box, GCN4_motif, O2-site, MBSI and circadian, etc. (3) Phytohormone response elements, including methyljasmonate response elements TGACG-motif and CGTCA-motif, (ABRE), gibberellin responsive element P-box, etc. (4) Stress response elements include LTR, MBS, ARE and mechanical damage response element W-box, etc.

Expression profile of different tissues and abiotic stress responses

To explore the possible biological functions of GbCDPKs, the transcriptome data of G. barbadense were used to analyze the expression characteristics of gene family in ovule and fiber at 0, 1, 3, 5, 10 and 20 days after anthesis, as well as in root, stem, leaf, receptacle, petal, stamen, calycle and sepals (Fig. 5A). The results showed that the CDPK gene has tissue-specific expression and spatio-temporal expression characteristics. According to the expression characteristics, the CDPK family can be divided into three types: the first group of genes are generally expressed less, but some genes were expressed in some or specific tissue. The second group was highly expressed in most tissues and developmental stages. The third group of genes are mainly expressed in vegetative organs or floral organs, but less expressed in ovule and fiber. Generally, homologous genes have similar expression patterns. Some homologous genes, such as GbCDPK13, GbCDPK56 and GbCDPK30, GbCDPK72 are highly expressed in the stamens. GbCDPK28, GbCDPK70; GbCDPK33, GbCDPK75 and GbCDPK38, GbCDPK80 are predominantly expressed in ovules. GbCDPK1, GbCDPK43 and GbCDPK25, GbCDPK67 are predominantly expressed in receptacles and stamens. The above results indicate that the CDPK gene family is widely involved in the growth and development of G. barbadense.

CDPKs can regulate the response of plants to abiotic stress. The transcriptome data was used to analyze the expression characteristics of the CDPK gene of G. barbadense under low temperature, high temperature, drought (PEG6000), and salt treatment (Fig. 5B). According to the expression characteristics shown in response to abiotic stress, the CDPK gene family can be divided into four groups: the first group most genes are low expressed under adversity stress. The second group all genes is low expressed. The three groups are all highly expressed. The fourth group most genes are highly expressed. In addition, homologous genes have similar expression profiles. For example, some homologous genes GbCDPK31 and GbCDPK73 are highly expressed under all stresses.

Expression under drought stress

The transcript levels of CDPK genes in cluster 1 and 2 were decreased under drought. In order to further understand the expression of CDPK genes at transcriptional level in G. barbadense. qRT-PCR was performed using RNA isolated from drought stress (Fig. 6). After drought, the expression of GbCDPK4 showed a downward-increasing-decreasing trend, and the expression was the lowest at 6 h. GbCDPK19 and GbCDPK78 showed an overall upward-decreasing-increasing-decreasing trend. The former had the highest expression at 12 h, and the latter at 12 h. The expression level of GbCDPK23 was the highest at 3 h. The expression level of GbCDPK23 showed a trend of decrease-increasing-decreasing, and the expression level of 12 h was the highest. The expression levels of GbCDPK44 and GbCDPK66 showed a trend of increase-decrease-increasing, and the expression level of the former was the lowest at 6 h. The expression was highest at 24 h.

Silencing GbCDPKs compromises cotton resistance to drought stress

In order to further study the function of CDPKs drought tolerance in G. barbadense, we chose to design specific primers in the 3′-UTR region of CDPKs to construct VIGS vectors. The transcript levels of CDPK genes in cluster 3 and 4 were significantly increased under drought, with the highest peak observed all stage of the drought treatment. The empty vector (TRV-00) and the vector with target fragments (TRV-CDPKs) were activated and propagated. Then the cotton cotyledons were injected into the 10 days seedlings for gene silencing. After two weeks of injection, when the positive seedlings were albino (Fig. 7A), indicating success of the VIGS experiment. The observed gene expression patterns may reflect their functions. The interference efficiency of the target gene was detected by qPCR technology. The results show that the GbCDPKs genes (CDPK32, CDPK68, CDPK74, CDPK80 and CDPK83) in cotton plants injected with bacterial liquid containing the target gene fragments can indeed interfere effectively and specifically (Fig. 7B).

To determine the tolerance of GbCDPKs silenced plants under drought stress, leaf discs were made from the same parts of TRV-00 and TRV-GbCDPKs (TRV-CDPK32, TRV-CDPK68, TRV-CDPK74, TRV-CDPK80 and TRV-CDPK83). The leaves were incubated in 0 and 20% PEG6000 solution for 4 days. The results showed that the leaf chlorosis and wilting of gene silenced TRV-CDPK32, TRV-CDPK68, TRV-CDPK74, TRV-CDPK80 and TRV-CDPK83 plants were significantly more chlorosis and wilting than that of the control (Fig. 7C).

In addition, we measured several physiological indicators including RWC, REL, water loss rate and MDA, proline and soluble sugar content under drought treatment and normal growth conditions. As shown in Figs. 7D–7H, in normal conditions, relative water content, electrical conductivity, proline and soluble sugar content of TRV-00 were higher than TRV-GbCDPKs. However, the TRV-00 of MDA was lower than silencing. Under the drought stress, the relative water content, proline and soluble sugar content were decreased in TRV-CDPKs plants, compared with the TRV-00. The relative electrical conductivity, MDA of TRV-00 and TRV-CDPKs both increased. From Fig. 7I, TRV-CDPKs was significantly higher than TRV-00, indicating that the drought resistance of plants decreased remarkably after silenced. In conclusion, the drought tolerance of the plants is weaker after silenced.