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

Abstract

Background

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

Results

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.

Conclusions

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²⁺ signals exert their functions through changes in the Ca²⁺ 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 Ca²⁺ receptor proteins contain helix-loop-helix fold (EF-hand) motifs that act as Ca²⁺-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 Ca²⁺ 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²⁺ and CaM and are the target proteins of the direct action of CaM. CaMBPs can be divided into two classes that are Ca²⁺-dependent or Ca²⁺-independent (Reddy et al. 2011). The IQ motif was the first identified Ca²⁺-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 Ca²⁺-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 Ca²⁺ 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, 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 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 (http://pfam.xfam.org), and the SMART databases (http://smart.embl-heidelberg.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). 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, 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) (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 ddH₂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 (Sun et al. 2018).

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

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^-△△C_T method. The online software Omicshare (https://www.omicshare.com/tools/Home/Soft/heatmap) 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^-ΔΔC_T method

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²⁺ 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 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²⁺ 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.

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.