Skip to main content

Genome-wide identification and expression analysis of the GhIQD gene family in upland cotton (Gossypium hirsutum L.)



Calmodulin (CaM) is one of the most important Ca2+ signaling receptors because it regulates diverse physiological and biochemical reactions in plants. CaM functions by interacting with CaM-binding proteins (CaMBPs) to modulate Ca2+ signaling. IQ domain (IQD) proteins are plant-specific CaMBPs that bind to CaM by their specific CaM binding sites.


In this study, we identified 102 GhIQD genes in the Gossypium hirsutum L. genome. The GhIQD gene family was classified into four clusters (I, II, III, and IV), and we then mapped the GhIQD genes to the G. hirsutum L. chromosomes. Moreover, we found that 100 of the 102 GhIQD genes resulted from segmental duplication events, indicating that segmental duplication is the main force driving GhIQD gene expansion. Gene expression pattern analysis showed that a total of 89 GhIQD genes expressed in the elongation stage and second cell wall biosynthesis stage of the fiber cells, suggesting that GhIQD genes may contribute to fiber cell development in cotton. In addition, we found that 20 selected GhIQD genes were highly expressed in various tissues. Exogenous application of MeJA significantly enhanced the expression levels of GhIQD genes.


Our study shows that GhIQD genes are involved in fiber cell development in cotton and are also widely induced by MeJA. Thw results provide bases to systematically characterize the evolution and biological functions of GhIQD genes, as well as clues to breed better cotton varieties in the future.


Calcium signaling is one of the most important cytosolic second messages that mediates various developmental processes and the responses to biotic and abiotic stresses (Reddy et al. 2011). Cytoplasmic Ca2+ signals exert their functions through changes in the Ca2+ concentration with spatiotemporal specificity (Liu et al. 2019) and can be induced by extracellular stimuli such as drought, salt-alkali stress, and light, or intracellular stimuli such as plant hormones and pathogenic factors (Hernandez et al. 2004; Yuan et al. 2017; Salveson et al. 2019). Most Ca2+ receptor proteins contain helix-loop-helix fold (EF-hand) motifs that act as Ca2+-binding domains. In higher plants, the calcium receptor proteins can be divided into four categories (Wei et al. 2019): Calmodulin (CaM), CaM-like proteins (CML), calcineurin B-like proteins (CBL), and calcium-dependent protein kinases (CDPK) (Ma et al. 2014), all of which contain EF-hand motifs.

CaM is one of the most important Ca2+ signaling receptors (Zhou et al. 2018), and it regulates diverse physiological and biochemical reactions. However, CaM has no enzymatic activity and it transmits signals by interacting with CaM-binding proteins to modulate cellular physiology (DeFalco et al. 2009). CaM-binding proteins (CaMBPs) play an important role between Ca2+ and CaM and are the target proteins of the direct action of CaM. CaMBPs can be divided into two classes that are Ca2+-dependent or Ca2+-independent (Reddy et al. 2011). The IQ motif was the first identified Ca2+-independent CaM-binding motif. In plants, the proteins containing IQ motifs include the myosin protein family, the calmodulin-binding transcription activator (CAMTA) protein family, the cyclic nucleotide-gated channel (CNGC) protein family, the IQ-motif (IQM) containing protein family, and the IQ67-domain (IQD) containing protein family (Steffen et al. 2013). IQD gene family members, which encoding plant-specific CaM/CML-binding proteins (CaMBPs), were firstly reported in Arabidopsis and rice (Abel et al. 2005). They are characterized by domains consisting of 67 amino acid residues (aa), such as the IQ67 domain, that are defined by a unique repetitive arrangement of the IQ motif; the Ca2+-dependent CaM recruitment motifs exhibit 1–5-10 and 1–8-14 arrangements (Burstenbinder et al. 2013).

Plant-specific IQD gene families have been analyzed in Populus trichocarpa (Ma et al. 2014), Arabidopsis thaliana, Oryza sativa (Abel et al. 2005), Phyllostachys edulis (Wu et al. 2016), Glycine max (Feng et al. 2014), and Solanum lycopersicum (Huang et al. 2013), and the functions of a few IQD genes have been reported. In tomato, the SUN genes and other members of the IQD gene family exert their effects on organ shape by interacting with microtubules (Steffen et al. 2013). In A. thaliana, AtIQD5 regulates pavement cell morphogenesis via Ca2+ signals (Liang et al. 2018). The tomato IQD gene SUN24 regulates seed germination through the abscisic acid (ABA) signaling pathway (Bi et al. 2018). The AtIQD1 protein localizes to microtubules and interacts with kinesin in light chain-related protein 1 (KLCR1) to facilitate the cellular transport of specific cargo (Burstenbinder et al. 2013). PtIQD genes showed tissue-specific expression patterns and could also be regulated by drought and methyl jasmonate (MeJA) stresses (Ma et al. 2014). Additionally, some of the GmIQD genes expressed specifically and could be regulated by MeJA stress (Feng et al. 2014).

Gossypium hirsutum L. (AADD, 2n = 4x = 52) is one of the most important economic crops worldwide because it provides major raw materials for the textile industry as well as edible oil. The fiber length is determined during its elongation stage (0 to 26 days post-anthesis (DPA) to reach its final length (Zhao et al. 2019). The cellulose of the fiber cell wall is mainly synthesized from 15 to 40 DPA, which determines the fiber strength (Song et al. 2019). Therefore, G. hirsutum L. fiber production and quality are very important. However, the production is greatly affected by abiotic and biotic stresses. For example, in China, aphid infestation has been found to reduce G. hirsutum L. production by 30% (Li et al. 2018), and salt stress has been shown to decrease cotton fiber production by 20% (Shaban et al. 2018). IQD proteins are CaMBPs that play an important role in plant stress signal transduction. However, these proteins have not been reported in G. hirsutum L. In this study, we identified 102 GhIQD genes and determined their chromosomal locations, predicted protein physicochemical properties, duplications, phylogenetic relationships, and expression patterns during fiber development at 5, 10, 15, and 25 DPA. Twenty selected GhIQD genes were used for the analysis of tissue-specific expression patterns and their response to MeJA stress. These preliminary results for the GhIQD genes provide the foundation for further researches on the physiological and biochemical functions of IQD proteins in G. hirsutum L.

Materials and methods

Identification of GhIQD gene family members in G. hirsutum L.

To identify GhIQD gene family members in G. hirsutum L., the G. hirsutum L. genome (Hu et al. 2019) sequences were downloaded from the Cotton Functional Genomics Database (CottonFGD, The AtIQD protein sequences from A. thaliana were downloaded from The Arabidopsis Information Resource (TAIR,, and the IQD protein sequences from Glycine max and Solanum lycopersicum were downloaded from Phytozome ( The G. hirsutum L. genome sequences were searched using a Blastp search against reported IQD proteins with e-value less than 1e-5 (Chen et al. 2017; Sun et al. 2017), including 34 SlSUN, 27 OsIQD, 66 GmIQD, and 33 AtIQD proteins. We further used a Hidden Markov Model (HMM) with the default parameters to search the G. hirsutum L. genome for IQD proteins (PF00612) in the Pfam database (, and the SMART databases ( were used to confirm that all the candidate sequences were members of the IQD family (Lin et al. 2014). The identified GhIQD genes were named according to their physical locations on the chromosomes from G. hirsutum and visualized by MapChart 2.2 software (Voorrips 2002).

The isoelectric point (pI) and molecular weight (MW) were predicted for each protein with the online software ExPASy ( The subcellular localization of GhIQD proteins was predicted using WoLF PSORT (

GhIQD protein sequence alignment and phylogenetic analysis

Multiple alignments of all the predicted IQD protein sequences from maize, soybean, bamboo, Arabidopsis, tomato, and G. hirsutum L. was performed with MEGA version 7 (Kumar et al. 2016; Qin et al. 2018).

Gene duplication and synteny analysis of GhIQD genes

The duplication pattern of each GhIQD gene was analyzed using MCScanX software according to the manual (Wang et al. 2012). The whole-genome BLASTP analysis of G. hirsutum L. was performed using local Blast software considering e-values less than 1e-5, and output was produced (Knip 2012). The Blast search outputs and the positions of all protein-coding genes were imported into MCScanX software (, and the genes were classified into the various types of duplications, including segmental, tandem, proximal, and dispersed duplications (Jia et al. 2019), using the default parameters (You et al. 2015). Synteny relationships were visualized with CIRCOS software (Krzywinski et al. 2009). Non-synonymous (Ka) and synonymous (Ks) substitution rates and the Ka/Ks ratio were estimated using DnaSP v5 software (Librado and Rozas 2009).

Expression analysis of the GhIQD genes

In the present study, RNA-seq data were downloaded from the public Cotton Functional Genomics Database (CottonFGD,, and the data were then used to survey the expression of the GhIQD genes (Zhu et al. 2017). The accession numbers of the RNA-Seq data for 5, 10, 20, and 25 DPA are SRR1695191, SRR1695192, SRR1695193, and SRR1695194, respectively, and all the expression values were standardized to fragments per kilobase per million (FPKM) values (Dong et al. 2018). The heatmap was performed to visualize gene expression patterns using OmicShare tools ( R software was used to visualize the gene expression profiles, and the TCseq package (Feng et al. 2019) was used to cluster the GhIQD gene expression patterns.

Plant materials and treatments

The G. hirsutum L. cultivar TM-1 was grown in the field in Anyang, Henan province, China. Leaves, stems, roots, and hypocotyl tissues were collected at the seedling stage. Stigma, petal, pollen, and calyx samples were collected at the flowering stage. The G. hirsutum L. cultivar TM-1 was also grown in a greenhouse under a 14 h light / 10 h dark photoperiod at 30 °C (day) and 28 °C (night) (Wang et al. 2018). Seedlings at the five-leaf stage were sprayed with 0.5 mM MeJA and the fifth leaf was sampled at seven time points (0 h, 3 h, 6 h, 12 h, 24 h, 48 h, and 72 h) after treatments (Yang et al. 2017).

The cotton cultivars TM-1 was obtained from the State Key Laboratory of Cotton Biology, Institute of Cotton Research of Chinese Academy of Agricultural Sciences. All tissue samples were immediately frozen in liquid nitrogen and stored at − 80 °C, and three biological replicates were conducted for each sample.

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

Total RNA from the collected samples was extracted using the Tiangen RNAprep Pure Plant Plus Kit (Tiangen, China) as directed by the manufacturer. First-strand cDNA was synthesized via reverse transcription of 2 μg of total RNA using the PrimeScript™ RT reagent Kit with gDNA Eraser (TaKaRa, Japan). Oligo 7 software was used to design gene-specific primers for qRT-PCR (Additional file 1: Table S1). The GhHis3 gene (AF024716) was used as an internal reference control for gene expression (Wang et al. 2013). The qRT-PCR experiments were performed with the TB Green™ Premix Ex Taq™ II RNaseH Plus kit (TaKaRa, Japan) on an ABI7500 real-time PCR system (Applied Biosystems, USA) with three replicates per sample. qRT-PCR assays were performed in a volume of 20 μL, which contained 2 μL of each primer, 1 μL of cDNA and 7 μL of ddH2O. The amplification conditions were as follows: initial denaturation at 95 °C for 2 min (Step 1), followed by 40 cycles of 10 s at 95 °C, 15 s at 58 °C, and 15 s at 72 °C (Step 2). The relative expression levels of the GhIQD genes were calculated using the 2-CT method (Livak and Schmittgen 2001). Statistical analyses were conducted using the one-way analysis of variance as implemented in SPSS software (Sun et al. 2018).


Identification of GhIQD gene family members in G. hirsutum L.

In order to identify IQD gene family members in G. hirsutum L., 33 Arabidopsis IQD protein sequences were used as queries to search in the upland cotton genome database. After the further selection of conserved domains, a total of 102 IQD genes from the G. hirsutum L. genome were identified as members of the GhIQD gene family. The chromosomal locations of the GhIQD genes were then determined using the upland cotton genome information (Hu et al. 2019). As a result, all of the GhIQD genes were mapped to the G. hirsutum L. chromosomes and named GhIQDA01.1 to GhIQDD13.12 based on their relative positions on the chromosomes (Table 1, Fig. 1). The number of amino acids (aa) in the predicted GhIQD protein sequences ranged from 120 (GhIQDA02.3) to 900 aa (GhIQDA05.6) with an average length of 458 aa, and the open reading frames (ORFs) ranged from 363 base pairs (bp) to 2 703 bp with an average length of 1 377 bp. The molecular weights (MWs) of the proteins encoded by these proteins varied from 13 689.84 Da (Daltons) (GhIQDA02.3) to 99 360.68 Da (GhIQDD05.6), with an average MW of 51 151.02 Da. Based on isoelectric point (pI) analysis, the calculated pIs of the 96 GhIQD genes were > 7.0 (with an average of 10.27), whereas six GhIQD genes were predicted to encode proteins with pIs < 7.0 (average of 6.05), including GhIQDA01.2, GhIQDA05.6, GhIQDA06.1, GhIQDD01.2, GhIQDD05.6, and GhIQDD06.1. The predicted subcellular localizations showed that 82 GhIQD proteins localized to the nucleus, nine GhIQD proteins localized to the mitochondria, nine GhIQD proteins localized to the chloroplasts, and two GhIQD proteins were found localized to the endoplasmic reticulum (ER) (Table 1).

Table 1 Information of GhIQD gene family members in G.hirsutum L., their sequence characteristics and subcellular location
Fig. 1

Chromosomal locations of cotton GhIQD genes on the G. hirsutum L. chromosomes. A01-A13 and D01-D13 indicate chromosomes from the A-subgenome and D-subgenome, respectively. The chromosome number is shown above each chromosome, and the relative locations of the GhIQD genes are indicated on the chromosomes

Phylogenetic analysis of GhIQD proteins

To examine the molecular evolutionary relationships among plant IQD proteins, the amino acid sequences of the IQD proteins from Arabidopsis, tomato, soybean, and G. hirsutum L. were used in phylogenetic analysis. As shown in Fig. 2, a phylogenetic tree was constructed with the Neighbor-joining (NJ) method from an alignment of all complete IQD protein sequences. The NJ tree showed that the IQD proteins group could be divided into four clusters (I, II, III, and IV). Cluster III is the largest one with 18 GhIQDs and cluster Ib is the smallest subcluster with 2 GhIQDs.

Fig. 2

Phylogenetic and evolutionary analysis of IQD proteins from different plant species. Note, green dots indicate G. hirsutum L. IQD proteins; pink dots indicate soybean IQD proteins; red dots indicate IQD proteins from tomato; blue dots indicate Arabidopsis IQD proteins. The phylogenetic tree was generated from an alignment of the IQD protein sequences using the Neighbor-joining (NJ) method in the MEGA 7 software package

IQD proteins are reported to specifically bind to calcium via CaM-binding sites (Cunwu et al. 2018). To better explore the biological functions of GhIQD proteins, the CaM-binding sites of the GhIQD proteins were predicted using the online Calmodulin Target Database software. As a result, GhIQD proteins are predicted to contain CaM-binding sites. Multiple consecutive strings of amino acid residues with scores > 7 are given in Additional file 1: Table S2. This result suggests that all GhIQD proteins contain CaM-binding sites with 1–3 strings of high-scoring amino acid residues.

Evolutionary analysis of GhIQD genes

To analyze the evolution of IQD gene family in G. hirsutum L., gene duplication events were analyzed. Based on the whole-genome duplication analysis in G. hirsutum L., 5 926 (8.14%) and 55 707 (76.56%) genes originated from tandem and segmental duplication, respectively. Therefore, we investigated the role of duplication events in the evolution of GhIQD genes. As shown in Fig. 3, among the 102 GhIQD genes identified in G. hirsutum L., 100 (98.04%) derived from segmental duplication events, and only two genes (GhIQDA13.3 and GhIQDD13.3) derived from proximal duplication. In contrast, none of the GhIQD genes was found to have arisen from tandem duplication events (Additional file 1: Table S3 and 4). These results indicate that segmental duplication is the main driving force in the expansion of the GhIQD genes.

Fig. 3

Circos figure of GhIQD gene pairs that arose from segmental duplication. The G. hirsutum L. chromosomes from the A- and D-subgenomes are shown in different colors, and the gene pairs involved in segmental duplication are linked by blue lines. The CIRCOS genome visualization tool was used to construct this figure

A Ka/Ks ratio > 1 indicates that paralogous gene pairs were produced by positive selection, a ratio < 1 indicates that paralogous gene pairs were under purifying selection and a ratio equal to 1 indicates that paralogous gene pairs were not subjected to selection pressure (Verma et al. 2017). To explore the type of selection pressure experienced by the duplicated GhIQD genes, paralogous GhIQD gene pairs were used to calculate synonymously (Ks) and non-synonymously (Ka) substitution rates to assess the ratio of non-synonymous to synonymous substitutions. As shown in Fig. 2, 50 paralogous gene pairs were identified. The Ka/Ks ratios of 48 members were < 1.0, and the Ka/Ks ratios for the remaining two paralogous gene pairs were > 1 (Additional file 1: Table S5), suggesting that the GhIQD paralogous gene pairs were mainly produced by purifying selection. Furthermore, the Ks ratio is stable and is usually used to estimate the evolution divergence time. The Ks ratios of 50 GhIQD gene pairs ranged from 0.009 40 to 0.379, and duplication events occurred from approximately 1.81 million years ago (MYA) to 72.87 MYA.

Transcriptome analysis of GhIQD genes during fiber development

Gene expression profiling can provide us with clues about the possible biological functions of genes. Therefore, we analyzed the gene expression profiles of GhIQD genes using the transcriptome data downloaded from the publicly available CottonFGD database. As shown in Fig. 4, a total of 20 GhIQD genes expressed during the developmental process in fiber cells. Cotton fiber development is divided into four overlapping stages: initiation, elongation, secondary cell wall deposition, and maturation (Tuttle et al. 2015). Based on the heatmap, six clusters of GhIQD genes are predominately expressed in cotton fiber cells (Fig. 4a). In detail, the 13 GhIQD genes in cluster 2 were highly expressed in fiber cells at 5 DPA; the expression levels of 17 GhIQD family members in cluster 5 were up-regulated in the 10 DPA samples; 21 genes in cluster 3 were significantly expressed in fibers at 20 DPA; genes in cluster 1 with 21 members were highly expressed in 25 DPA fiber cells; in cluster 4, the transcripts of seven genes were abundant in fibers at 20 and 25 DPA; and the remaining 10 GhIQDs in cluster 6 were highly expressed in the 5 DPA and 25 DPA samples (Fig. 4b). These results imply that GhIQD genes may function in fiber cell development in cotton.

Fig. 4

Expression profiling of GhIQD genes during four developmental stages of the cotton fiber cell. a A hierarchical clustering heatmap of GhIQD gene expression in fiber cells at 5, 10, 20, and 25 DPA. These data were obtained from the transcriptome data available on the Cottongen website ( The color scale represents the expression values, which have been normalized by the online software Omicshare b Gene expression of 89 GhIQD genes in six different clusters during four developmental stages of cotton fiber cells

Tissue-specific expression analysis of GhIQD genes by qRT-PCR

The GhIQD gene family has 102 members; according to the expression abundance of GhIQDs in transcriptomes, 20 genes were selected to investigate their expression patterns in different tissues, including the calyx, leaf, stigma, stem, root, petal, pollen, and hypocotyl. As shown in Fig. 5, GhIQDA13.1, GhIQDD12.1 and GhIQDD13.1 were predominantly expressed in pollen (Fig. 5d), indicating that these genes may play pivotal roles in pollen development. GhIQDD01.3, GhIQDD01.2 and GhIQDD05.2 showed stem-specific expression (Fig. 5b), and GhIQDA01.1, GhIQDA05.2 and GhIQDA08.1 expressed preferentially in leaves (Fig. 5a). The GhIQDA06.1, GhIQDD06.1 and GhIQDD09.1 genes showed higher expression levels in leaves and stems (Fig. 5c). Most genes investigated were abundantly expressed in different tissues (Fig. 5e), as observed for GhIQDA02.1 and GhIQDD02.1 that were highly expressed in all tissues with similar expression patterns. The cluster II genes, GhIQDA12.3 and GhIQDD12.4, were highly expressed in the leaf, petal, and hypocotyl (Fig. 5e). The expression patterns in specific tissues are strong evidence for roles of particular GhIQD genes in these specific locations and developmental processes.

Fig. 5

Tissue-specific expression analysis of 20 selected GhIQD genes by qRT-PCR. The relative expression levels of 20 GhIQD genes were examined by qRT-PCR assays and normalized to the expression level of the house-keeping gene GhHis3. The calyx, leaf, stigma, stem, root, petal, pollen, and hypocotyl tissues are indicated on the X-axis. Relative gene expression levels compared with the calyx are indicated on the Y-axis, * represents a significant difference at the p < 0.05 level, ** represents an extremely significant difference at the p < 0.01 level

Expression profiling of GhIQD genes in response to MeJA treatment

According to previous studies, the expression of most IQD genes can be induced by MeJA stress in plants (Bi et al. 2018). In this study, the expression patterns of GhIQD genes in plants exposed to MeJA treatment were examined in a qRT-PCR experiment. The results showed that the expression levels of the 20 selected GhIQD genes were significantly increased by MeJA treatment in the leaf (Fig. 6). As the time of treatment increased, the transcript levels for most genes increased significantly. In detail, the expression levels of GhIQDA09.3, GhIQDA13.1, GhIQDD01.3, GhIQDD02.1, GhIQDD05.2, GhIQDD09.1, and GhIQDD13.1 were induced from 0 h, with the highest expression levels detected at 12 and 24 h after the MeJA treatment. The GhIQDA06.1, GhIQDD06.1, and GhIQDD09.4 genes also exhibited the highest expression levels at 24 h after the MeJA treatment. The expression levels of GhIQDA01.1, GhIQDA02.1, GhIQDA08.1, GhIQDA12.3, GhIQDA13.4, GhIQDD01.2, GhIQDD12.1, and GhIQD12.4 peaked at 6 h after the treatment. Compared with the other genes, the maximum expression of GhIQDA05.5 occurred at 72 h. These results showed that the GhIQD genes are widely induced by MeJA treatment.

Fig. 6

The expression profiles of 20 selected GhIQD genes in response to MeJA treatment at 0, 3, 6, 12, 24, 48, and 72 h. The gene expression values were determined by qRT-PCR and calculated using the 2-CT method. The online software Omicshare ( was used to normalize the expression values and draw the heatmap. The color scale on the top of the figure represents expression values computed by the 2-ΔΔCT method


Calcium is one of the most important cytosolic second messengers, and calcium levels can be induced by intracellular and extracellular stimuli. CaM is one of the most important Ca2+ signaling receptors regulating diverse physiological and biochemical reactions. CaM functions by interacting with CaM-binding proteins (CaMBPs) to modulate cellular physiology. IQD proteins are plant-specific CaM/CML CaMBPs that are characterized by 67-amino acid domains. In this study, we identified 102 GhIQD genes in G. hirsutum and analyzed their chromosomal locations, protein physicochemical properties, gene duplication events, phylogenetic relationships, and expression patterns during the development of fiber cells. Twenty selected GhIQD genes were used for the analysis of tissue-specific expression patterns and their response to MeJA treatment.

The GhIQD gene family expanded by segmental duplication

A number of IQD genes have been reported in different plants; there are 33 AtIQD genes in A. thaliana, 28 OsIQD genes in O. sativa (Abel et al. 2005), 38 PtIQD genes in P. trichocarpa (Ma et al. 2014), 29 PeIQD genes in P. edulis (moso bamboo) (Wu et al. 2016), 67 GmIQD genes in Glycine max (Feng et al. 2014), and 34 SlSUN/IQD genes in Solanum lycopersicum (Huang et al. 2013). In this study, a total of 102 GhIQD genes were identified in G. hirsutum L. The number of IQD genes in G. hirsutum L. is greater than that found in other plant species, possibly because G. hirsutum L. is an allotetraploid cotton species that originated from the hybridization of G. arboreum and G. raimondii and subsequent polyploidization 12 million years ago (Wang et al. 2019). A whole-genome duplication analysis showed that 100 of the GhIQD genes arose from segmental duplication, and other two genes, GhIQDA13.3 and GhIQDD13.3, originated from proximal duplications. Additionally, Ka/Ks analysis indicated that most of the GhIQD genes were under purifying selection, which indicates that the segmentally duplicated GhIQD genes were subjected to strong purifying constraints during evolution. The divergence time of GhIQDs was earlier than 1.81 MYA, which indicates that the duplication events of GhIQDs occurred before the polyploidization events in G. hirsutum.

GhIQD genes participate widely in the regulation of growth in G. hirsutum L.

In Arabidopsis, the AtIQD proteins were reported to widely link calcium signaling to microtubules, membrane subdomains, and the nucleus (Bürstenbinder et al. 2017a). From the results of the transcriptome analysis, the GhIQD gene expression patterns could be clustered into six groups. In the fiber development process, the 5 DPA ovule stage is the primary cell wall synthesis stage of fiber cells; 10 DPA corresponds to the elongation stage of fiber development; and 2025 DPA is the transition stage of fiber development from elongation to secondary wall synthesis, which is important for fiber strength (Hu et al. 2013). PdIQD10 gene was found to be involved in the secondary cell wall biosynthesis and biomass formation in Populus (Badmi et al. 2018). Therefore, the genes in cluster 6 might participate in primary cell wall synthesis, the GhIQD genes in cluster 5 may contribute to fiber elongation, and the GhIQD genes in clusters 1, 3, and 4 may be involved in fiber strength development.

In Arabidopsis, AtIQD genes function as hubs in Ca2+ signaling to regulate growth and development tissue-specifically (Bürstenbinder et al. 2017b). To further elucidate the possible functions of the GhIQD genes, their expression patterns were investigated in various tissues (Fig. 5). The results show that some GhIQD genes are predominantly expressed in specific tissues, and the paralogous gene pairs exhibit similar expression patterns, such as GhIQDA13.1 and GhIQDD13.1, which are predominantly expressed in pollen, indicating that GhIQD gene pairs are functionally redundant.

MeJA is an ester of jasmonic acid and is widely present in plants (Yan et al. 2018). MeJA triggers the biosynthesis of plant defensive compounds and initiates the expression of pathogenesis-related genes involved in systemic acquired resistances and local resistances (Pei 2017). In this study, the expression levels of all 20 selected GhIQD genes were increased (Fig. 6) in response to MeJA treatment, which is consistent with previous studies showing that most of the IQD genes in P. trichocarpa and moso bamboo are induced by MeJA treatment (Ma et al. 2014; Wu et al. 2016). These results imply that the IQD gene family plays an important role in tolerance to MeJA abiotic stress in plants.


In conclusion, we identified 102 GhIQD genes in G. hirsutum L. Segmental duplication was the main driving force behind the expansion of the GhIQD family. Ka/Ks analysis showed that GhIQD genes were under purifying selection. Based on the expression analysis, 89 genes could be detected during the stages of fiber development, tissue-specific expression analysis showed that some of the GhIQD genes were specifically expressed, and all 20 selected GhIQD genes could be induced by MeJA treatment. These preliminary results provide the foundation for further research on the physiological and biochemical functions of IQD proteins in G. hirsutum L.

Availability of data and materials

The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.


  1. Abel S, Savchenko T, Levy M. Genome-wide comparative analysis of the IQD gene families in Arabidopsis thaliana and Oryza sativa. BMC Evol Biol. 2005;5:72. .

  2. Badmi R, Payyavula RS, Bali G, et al. A new calmodulin-binding protein expresses in the context of secondary cell wall biosynthesis and impacts biomass properties in Populus. Front Plant Sci. 2018;9:1669. .

  3. Bi L, Lin W, Zhuyan J, et al. The tomato IQD gene SUN24 regulates seed germination through ABA signaling pathway. Planta. 2018;248:919–31. .

  4. Bürstenbinder K, Mitra D, Quegwer J. Functions of IQD proteins as hubs in cellular calcium and auxin signaling: a toolbox for shape formation and tissue-specification in plants? Plant Signal Behav. 2017b;12:1692–708. .

  5. Bürstenbinder K, Möller B, Plötner R, et al. The IQD family of calmodulin-binding proteins links calcium signaling to microtubules, membrane subdomains, and the nucleus. Plant Physiol. 2017a;173:01743.02016. .

    CAS  Article  Google Scholar 

  6. Burstenbinder K, Savchenko T, Muller J. Arabidopsis calmodulin-binding protein IQ67-domain 1 localizes to microtubules and interacts with kinesin light chain-related protein-1. J Biol Chem. 2013;288:1871–82. .

  7. Chen M, Li K, Li H, et al. The glutathione peroxidase gene family in Gossypium hirsutum: genome-wide identification, classification, gene expression and functional analysis. Sci Rep. 2017;7:44743. .

  8. DeFalco TA, Bender KW, Snedden WA. Breaking the code: Ca2+ sensors in plant signalling. Biochem J. 2009;425:27–40. .

    CAS  Article  PubMed  Google Scholar 

  9. Dong C, Wang J, Yu Y. Identifying functional genes influencing Gossypium hirsutum fiber quality. Front Plant Sci. 2018;9:1968. .

  10. Feng L, Chen Z, Ma H, et al. The IQD gene family in soybean: structure, phylogeny, evolution and expression. PLoS One. 2014;9:e110896. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  11. Feng S, Xu M, Liu F, et al. Reconstruction of the full-length transcriptome atlas using PacBio Iso-Seq provides insight into the alternative splicing in Gossypium australe. BMC Plant Biol. 2019;19:365. .

  12. Sebastià HC, Hardin SC, Clouse SD, et al. Identification of a new motif for CDPK phosphorylation in vitro that suggests ACC synthase may be a CDPK substrate. Arch Biochem Biophys. 2004;428:81–91. .

  13. Hu G, Koh J, Yoo MJ, et al. Proteomic profiling of developing cotton fibers from wild and domesticated Gossypium barbadense. New Phytol. 2013;200:570–82. .

  14. Hu Y, Chen JD, Fang L, et al. Gossypium barbadense and Gossypium hirsutum genomes provide insights into the origin and evolution of allotetraploid cotton. Nat Genet. 2019;51(4):739–48. .

  15. Huang Z, Van Houten J, Gonzalez G, et al. Genome-wide identification, phylogeny and expression analysis of SUN, OFP and YABBY gene family in tomato. Mol Gen Genomics. 2013;288:111–29. .

  16. Jia Y, Li B, Zhang Y, et al. Evolutionary dynamic analyses on monocot flavonoid 3′-hydroxylase gene family reveal evidence of plant-environment interaction. BMC Plant Biol. 2019;19:347. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  17. Knip M. The SLEEPER genes: a transposase-derived angiosperm-specific gene family. BMC Plant Biol. 2012;12:192. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  18. Krzywinski M, Schein J, Biro I, et al. Circos: an information aesthetic for comparative genomics. Genome Res. 2009;19:1639–45. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  19. Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33:1870. .

    CAS  Article  Google Scholar 

  20. Li JH, Wu YK, Zhang Q. Aphid parasitism and parasitoid diversity in cotton fields in Xinjiang, China. PLoS One. 2018;13:e0207034. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  21. Liang H, Yi Z, Pablo M, et al. The microtubule-associated protein IQ67 DOMAIN5 modulates microtubule dynamics and pavement cell shape. Plant Physiol. 2018:00558–2018. .

  22. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2. .

    CAS  Article  PubMed  Google Scholar 

  23. Liu T, Zhuang L, Huang B. Metabolic adjustment and gene expression for root sodium transport and calcium signaling contribute to salt tolerance in Agrostis grass species. Plant Soil. 2019;1:14. .

    CAS  Article  Google Scholar 

  24. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods. 2001;25:402–8. .

  25. Ma H, Feng L, Chen Z, et al. Genome-wide identification and expression analysis of the IQD gene family in Populus trichocarpa. Plant Sci. 2014;229:96–110. .

  26. Pei Y. Role of hydrogen sulfide in the methyl jasmonate response to cadmium stress in foxtail millet. Front Biosci (Landmark Ed). 2017;22:530–8. .

    Article  Google Scholar 

  27. Qin L, Mo N, Muhammad T, et al. Genome-wide analysis of DCL, AGO, and RDR gene families in pepper (Capsicum Annuum L.). Int J Mol Sci. 2018;19:1038. .

  28. Reddy AS, Ali GS, Celesnik H, et al. Coping with stresses: roles of calcium- and calcium/calmodulin-regulated gene expression. Plant Cell. 2011;23:2010–32. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  29. Salveson I, Anderson DE, Hell JW, et al. Chemical shift assignments of a calmodulin intermediate with two Ca2+ bound in complex with the IQ-motif of voltage-gated Ca2+ channels (CaV1.2). Biomol NMR Assign. 2019;13:233. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  30. Shaban M, Ahmed MM, Sun H, et al. Genome-wide identification of lipoxygenase gene family in cotton and functional characterization in response to abiotic stresses. BMC Genomics. 2018;19:599. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  31. Song ZQ, Chen Y, Zhang CY. RNA-seq reveals hormone-regulated synthesis of non-cellulose polysaccharides associated with fiber strength in a single-chromosomal-fragment-substituted upland cotton line. Crop J. 2019;8(2):273–86. .

    CAS  Article  Google Scholar 

  32. Steffen A, Katharina B, Jens M. The emerging function of IQD proteins as scaffolds in cellular signaling and trafficking. Plant Signal Behav. 2013;8(6):e24369. .

    CAS  Article  Google Scholar 

  33. Sun HR, Hao PB, Ma Q, et al. Genome-wide identification and expression analyses of the pectate lyase (PEL) gene family in cotton (Gossypium hirsutum L.). BMC Genomics. 2018;19:661. .

  34. Sun Q, Wang G, Zhang X. Genome-wide identification of the TIFY gene family in three cultivated Gossypium species and the expression of JAZ genes. Sci Rep. 2017;7:42418. .

  35. Tuttle JR, Nah GJ, Duke MV. Metabolomic and transcriptomic insights into how cotton fiber transitions to secondary wall synthesis, represses lignification, and prolongs elongation. BMC Genomics. 2015;16:477. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  36. Verma G, Dhar YV, Srivastava D. Genome-wide analysis of rice dehydrin gene family: its evolutionary conservedness and expression pattern in response to PEG induced dehydration stress. PLoS One. 2017;12:e0176399. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  37. Voorrips RE. MapChart: software for the graphical presentation of linkage maps and QTLs. J Hered. 2002;93:77–8. .

    CAS  Article  PubMed  Google Scholar 

  38. Wang LL, Yang ZE, Zhang B, et al. Genome-wide characterization and phylogenetic analysis of GSK gene family in three species of cotton: evidence for a role of some GSKs in fiber development and responses to stress. BMC Plant Biol. 2018;18:330. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  39. Wang M, Tu L, Yuan D, et al. Reference genome sequences of two cultivated allotertraploid cottons, Gossypium hirsutum and Gossypium barbadense. Nat Genet. 2019;51(2):224-9.

  40. Wang M, Wang Q, Zhang B. Evaluation and selection of reliable reference genes for gene expression under abiotic stress in cotton (Gossypium hirsutum L.). Gene. 2013;530:44–50. .

  41. Wang Y, Tang H, DeBarry JD, et al. MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity. Nucleic Acids Res. 2012;40:e49. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  42. Wei CC, Fabry E, Hay E, et al. Metal binding and conformational studies of the calcium binding domain of NADPH oxidase 5 reveal its similarity and difference to calmodulin. J Biomol Struct Dyn. 2019;5:1–29. .

    CAS  Article  Google Scholar 

  43. Wu M, Li Y, Chen D, et al. Genome-wide identification and expression analysis of the IQD gene family in moso bamboo (Phyllostachys edulis). Sci Rep. 2016;6:24520. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  44. Yan C, Fan M, Yang M, et al. Injury activates Ca2+/calmodulin-dependent phosphorylation of JAV1-JAZ8-WRKY51 complex for jasmonate biosynthesis. Mol Cell. 2018;70:136–149.e137. .

    CAS  Article  PubMed  Google Scholar 

  45. Yang Y, Chen T, Ling X, et al. Gbvdr6, a gene encoding a receptor-like protein of cotton (Gossypium barbadense), confers resistance to Verticillium wilt in Arabidopsis and upland cotton. Front Plant Sci. 2017;8:2272. .

  46. You J, Zhang L, Song B, et al. Systematic analysis and identification of stress-responsive genes of the NAC gene family in Brachypodium distachyon. PLoS One. 2015;10:e0122027. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  47. Yuan P, Jauregui E, Du L, et al. Calcium signatures and signaling events orchestrate plant–microbe interactions. Curr Opin Plant Biol. 2017;38:173–83. .

    CAS  Article  PubMed  Google Scholar 

  48. Zuo CW, Liu YL, Guo ZG, et al. Genome-wide annotation and expression responses to biotic stresses of the WALL-ASSOCIATED KINASE - RECEPTOR-LIKE KINASE (WAK-RLK) gene family in apple (Malus domestica). Eur J Plant Pathol. 2018;248:919–31. .

  49. Zhao W, Dong H, Zahoor R, et al. Ameliorative effects of potassium on drought-induced decreases in fiber length of cotton (Gossypium hirsutum L.) are associated with osmolyte dynamics during fiber development. Crop J. 2019:619–34. .

  50. Zhou YP, Wu JH, Xiao WH, et al. Arabidopsis IQM4, a novel calmodulin-binding protein, is involved with seed dormancy and germination in Arabidopsis. Front Plant Sci. 2018;9:721. .

  51. Zhu T, Liang C, Meng Z, et al. CottonFGD: an integrated functional genomics database for cotton. BMC Plant Biol. 2017;17:101. .

    CAS  Article  PubMed  PubMed Central  Google Scholar 

Download references


We would like to thank the anonymous reviewers for their valuable comments and helpful suggestions which help to improve the manuscript.


This research was sponsored by the State Key Laboratory of Cotton Biology Open Fund (grant numbers CB2019A03 and CB2018A07), comprehensive Scientific research fund project of Xianyang Normal University (XSYK20002), the Innovation and Entrepreneurship Training Program for College Students in Shaanxi Province (S202010722071), the National Natural Science Foundation of China (grant number 31872175), and Key Research and Development Program of Shaanxi Province (grant number 2019NY-103). The funding bodies provided financial support to the research projects but didn’t involve in study design, data collection, analysis, or preparation of the manuscript.

Author information




Formal analysis, Kang YY and Zou CS; Investigation, Tian RJ. Resources, Li HZ; Software, Pang CY and Shang HH; Validation, Li SY, Liu FP, Cao LY, Jin YH, Li JY, Huang DQ and Liu Y; Visualization, Lv LM and Wang WB; Original draft, Dou LL; Reviewing and editing, Song GL and Xiao GH. The author(s) read and approved the final manuscript.

Corresponding authors

Correspondence to SONG Guoli or XIAO Guanghui.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no conflict of interest.

Supplementary Information

Additional file 1

: Table S1. Gene-specific primer pairs used in the qRT-PCR experiments. Table S2. Predicted calmodulin-binding sites of IQD proteins in G. hirsutum L. The calmodulin-binding sites were predicted with the online Calmodulin Target Database software, and the table shows strings of amino acid residues with a score of at least 7. Residues with a score of 9 are highlighted in bold. The numbers before the strings and after the strings indicate the locations of the first and the last amino acid residues of the strings in the GhIQD protein, respectively. Table S3. List of segmentally duplicated IQD gene pairs in the G. hirsutum L. genome along with their e-values identified from MCScanX. Table S4. List of tandem and segmentally duplicated IQD genes in the G. hirsutum L. genome identified with MCScanX software. Table S5. Divergence between paralogous IQD gene pairs in G. hirsutum L.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

DOU, L., LV, L., KANG, Y. et al. Genome-wide identification and expression analysis of the GhIQD gene family in upland cotton (Gossypium hirsutum L.). J Cotton Res 4, 4 (2021).

Download citation


  • Gossypium hirsutum L.
  • GhIQD genes
  • Segmental duplication
  • Expression analysis