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.


Background
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 Ca 2+ signals exert their functions through changes in the Ca 2+ concentration with spatiotemporal specificity ) and can be induced by extracellular stimuli such as drought, saltalkali 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 Ca 2+ receptor proteins contain helix-loop-helix fold (EFhand) motifs that act as Ca 2+ -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) , all of which contain EF-hand motifs.
CaM is one of the most important Ca 2+ 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 Ca 2+ and CaM and are the target proteins of the direct action of CaM. CaMBPs can be divided into two classes that are Ca 2+ -dependent or Ca 2+ -independent (Reddy et al. 2011). The IQ motif was the first identified Ca 2+ -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/CMLbinding 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 Ca 2+ -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 , Arabidopsis thaliana, Oryza sativa (Abel et al. 2005), Phyllostachys edulis (Wu et al. 2016), Glycine max , 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 Ca 2+ 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 . Additionally, some of the GmIQD genes expressed specifically and could be regulated by MeJA stress .
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, https:// cottonfgd.org/about/download.html). The AtIQD protein sequences from A. thaliana were downloaded from The Arabidopsis Information Resource (TAIR, https://www.arabidopsis.org/index.jsp), and the IQD protein sequences from Glycine max and Solanum lycopersicum were downloaded from Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html). The G. hirsutum L. genome sequences were searched using a Blastp search against reported IQD proteins with evalue less than 1e-5 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 (http://pfam.xfam. org), and the SMART databases (http://smart.emblheidelberg.de) 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 (https://www.expasy.org/tools/). The subcellular localization of GhIQD proteins was predicted using WoLF PSORT (https://wolfpsort.hgc.jp).

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 (http://chibba.pgml.uga.edu/mcscan2/), 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). Nonsynonymous (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, https://cottonfgd.org/about/download. html), 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 (https://www.omicshare.com/tools/ Home/Soft/heatmap). 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) . 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 ddH 2 O. 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 -△△C T method (Livak and Schmittgen 2001). Statistical analyses were conducted using the one-way analysis of variance as implemented in SPSS software .

Results
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 (   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).

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. 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 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 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.
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,  (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.

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 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 (https://www.cottongen.org/). 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 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.

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.

Discussion
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 Ca 2+ signaling receptors regulating diverse 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 physiological and biochemical reactions. CaM functions by interacting with CaM-binding proteins (CaMBPs) to modulate cellular physiology. IQD proteins are plantspecific 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 , 29 PeIQD genes in P. edulis (moso bamboo) (Wu et al. 2016), 67 GmIQD genes in Glycine max , 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 1∼2 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 20∼25 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 Ca 2+ signaling to regulate growth and development tissuespecifically (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 Wu et al. 2016). These results imply that the IQD gene family plays an important role in tolerance to MeJA abiotic stress in plants.

Conclusions
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.
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.