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. hirsutumL., 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. hirsutumL., 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).

Gene ID	Gene name	length /aa	ORF length /bp	Mass /Da	pI	Predicted location	Exon number
GH_A01G1101	GhIQDA01.2	838	2 517	92 762.8	6.45	nucl	6
GH_A01G2332	GhIQDA01.3	466	1 401	51 680.07	10.04	nucl	5
GH_A02G0100	GhIQDA02.1	540	1 623	60 671.48	11.13	nucl	5
GH_A02G1393	GhIQDA02.2	385	1 158	43 210.99	10.19	nucl	3
GH_A04G0950	GhIQDA04.1	467	1 404	51 973.41	10.46	nucl	5
GH_A05G0263	GhIQDA05.1	411	1 236	45 861.25	10.34	nucl	3
GH_A05G0790	GhIQDA05.2	434	1 305	47 787.02	9.98	nucl	6
GH_A05G1339	GhIQDA05.3	403	1 212	44 740.05	10.14	nucl	3
GH_A05G1378	GhIQDA05.4	447	1 344	50 742.09	10.1	nucl	4
GH_A05G1577	GhIQDA05.5	306	921	34 962.1	10.42	nucl	5
GH_A05G1877	GhIQDA05.6	900	2 703	99 312.5	5.27	nucl	6
GH_A05G3022	GhIQDA05.7	339	1 020	37 749.85	10.8	nucl	4
GH_A05G3607	GhIQDA05.8	574	1725	52 858.66	10.51	nucl	5
GH_A05G4277	GhIQDA05.9	464	1 395	52 514.35	10.14	nucl	5
GH_A06G0050	GhIQDA06.1	659	1980	72 053.71	6.26	nucl	5
GH_A06G1837	GhIQDA06.2	442	1 329	50 167.53	9.86	nucl	5
GH_A06G1940	GhIQDA06.3	382	1 149	43 725.98	10.81	nucl	5
GH_A07G0094	GhIQDA07.1	414	1 245	45 320.56	10.21	nucl	3
GH_A07G0227	GhIQDA07.2	437	1 314	49 258.14	10.51	nucl	5
GH_A07G0720	GhIQDA07.3	651	1956	72 594.59	9.72	nucl	6
GH_A08G1336	GhIQDA08.2	539	1 620	60 952.52	9.78	nucl	4
GH_A08G2765	GhIQDA08.4	360	1 083	39 811.23	10.4	nucl	3
GH_A08G2816	GhIQDA08.5	418	1 257	47 014.72	10.52	nucl	5
GH_A09G1366	GhIQDA09.2	409	1 230	46 558.51	10.2	nucl	5
GH_A09G1685	GhIQDA09.3	452	1 359	50 179.15	10.12	nucl	5
GH_A09G2303	GhIQDA09.5	511	1 536	56 997.69	11.39	nucl	5
GH_A09G2394	GhIQDA09.6	404	1 215	45 045.51	9.94	nucl	3
GH_A10G0319	GhIQDA10.1	311	936	35 184.76	10.9	nucl	5
GH_A10G1560	GhIQDA10.2	523	1 572	58 681.88	10.66	nucl	5
GH_A11G0174	GhIQDA11.1	293	882	33 829.1	11.24	nucl	3
GH_A11G0321	GhIQDA11.2	767	2 304	84 830.18	9.77	nucl	6
GH_A11G1049	GhIQDA11.3	531	1 596	59 599.5	10.51	nucl	5
GH_A11G1587	GhIQDA11.4	381	1 146	43 192.29	9.28	nucl	6
GH_A12G1325	GhIQDA12.1	500	1 503	55 878.89	10.73	nucl	5
GH_A12G2604	GhIQDA12.2	519	1 560	59 126.62	9.88	nucl	5
GH_A12G2690	GhIQDA12.3	451	1 356	49 874.94	9.74	nucl	4
GH_A13G0768	GhIQDA13.2	549	1 650	60 060.1	9.99	nucl	6
GH_A13G0771	GhIQDA13.3	550	1 653	60 669.01	10	nucl	6
GH_A13G1240	GhIQDA13.5	306	921	34 548.98	10.4	nucl	3
GH_A13G1987	GhIQDA13.6	362	1 089	41 511.44	10.45	nucl	3
GH_A13G2167	GhIQDA13.7	326	981	38 128.29	10.38	nucl	5
GH_D01G1150	GhIQDD01.2	838	2 517	92 819.9	6.13	nucl	6
GH_D01G2410	GhIQDD01.3	466	1 401	51 732.91	10.07	nucl	5
GH_D02G0104	GhIQDD02.1	543	1 632	61 081.99	11.09	nucl	5
GH_D03G0607	GhIQDD03.2	385	1 158	43 307.07	10.25	nucl	3
GH_D04G0099	GhIQDD04.1	464	1 395	52 501.4	10.16	nucl	5
GH_D04G0861	GhIQDD04.2	474	1 425	52 953.74	10.51	nucl	5
GH_D04G1282	GhIQDD04.3	470	1 413	52 267.63	10.25	nucl	5
GH_D05G0268	GhIQDD05.1	411	1 236	45 920.33	10.28	nucl	3
GH_D05G0787	GhIQDD05.2	435	1 308	48 163.42	10.09	nucl	6
GH_D05G1342	GhIQDD05.3	403	1 212	44 653.93	10.36	nucl	3
GH_D05G1385	GhIQDD05.4	447	1 344	50 640.2	10.09	nucl	4
GH_D05G1607	GhIQDD05.5	306	921	34 947.05	10.39	nucl	5
GH_D05G1915	GhIQDD05.6	900	2 703	99 360.68	5.3	nucl	6
GH_D05G3028	GhIQDD05.7	427	1 284	47 736.21	10.29	nucl	5
GH_D06G0037	GhIQDD06.1	664	1995	72 524.53	6.89	nucl	5
GH_D06G1866	GhIQDD06.2	442	1 329	50 104.22	9.81	nucl	5
GH_D06G1971	GhIQDD06.3	386	1 161	44 138.29	10.65	nucl	5
GH_D07G0100	GhIQDD07.1	415	1 248	45 440.73	10.2	nucl	3
GH_D07G0234	GhIQDD07.2	437	1 314	49 274.22	10.48	nucl	5
GH_D07G0707	GhIQDD07.3	650	1953	72 097.97	9.75	nucl	6
GH_D08G1206	GhIQDD08.3	540	1 623	60 919.52	9.78	nucl	4
GH_D08G2761	GhIQDD08.5	359	1 080	39 873.37	10.43	nucl	3
GH_D08G2809	GhIQDD08.6	395	1 188	44 794.15	10.56	nucl	5
GH_D09G1141	GhIQDD09.1	436	1 311	47 553.9	10.08	nucl	6
GH_D09G1315	GhIQDD09.2	414	1 245	47 099.05	10.12	nucl	5
GH_D09G1630	GhIQDD09.3	446	1 341	49 661.4	10.1	nucl	5
GH_D09G2241	GhIQDD09.4	511	1 536	57 044.68	11.27	nucl	5
GH_D09G2334	GhIQDD09.5	406	1 221	45 218.73	9.98	nucl	3
GH_D10G0332	GhIQDD10.1	311	936	35 020.59	10.74	nucl	5
GH_D10G1333	GhIQDD10.2	523	1 572	58 604.9	10.63	nucl	5
GH_D11G0175	GhIQDD11.1	293	882	33 561.81	11.06	nucl	3
GH_D11G0334	GhIQDD11.2	767	2 304	84 784.04	9.8	nucl	6
GH_D11G1078	GhIQDD11.3	531	1 596	59 590.38	10.55	nucl	5
GH_D11G1617	GhIQDD11.4	411	1 236	47 034.9	9.46	nucl	7
GH_D12G1345	GhIQDD12.1	501	1 506	56 136.15	10.76	nucl	5
GH_D12G2626	GhIQDD12.3	518	1 557	58 961.41	9.89	nucl	5
GH_D12G2716	GhIQDD12.4	451	1 356	49 960.99	9.73	nucl	4
GH_D13G0741	GhIQDD13.2	522	1 569	60 358.36	9.98	nucl	6
GH_D13G0744	GhIQDD13.3	552	1 659	60 897.05	9.98	nucl	6
GH_D13G1947	GhIQDD13.5	376	1 131	43 085.2	10.43	nucl	3
GH_D13G2149	GhIQDD13.6	377	1 134	42 426.34	10.51	nucl	5
GH_A01G0547	GhIQDA01.1	285	858	32 302.82	10.15	mito	5
GH_A02G2016	GhIQDA02.3	120	363	13 689.84	11.79	mito	3
GH_A03G0236	GhIQDA03.1	488	1 467	53 976.96	9.96	mito	4
GH_A13G0956	GhIQDA13.4	480	1 443	53 068.34	9.98	mito	5
GH_D01G0541	GhIQDD01.1	285	858	32 287.85	10.11	mito	5
GH_D03G0046	GhIQDD03.1	133	402	15 440.04	11.65	mito	2
GH_D03G1733	GhIQDD03.3	489	1 470	54 209.22	10.02	mito	4
GH_D12G2109	GhIQDD12.2	315	948	36 162.76	10.74	mito	4
GH_D13G1022	GhIQDD13.4	480	1 443	52 890.99	9.94	mito	5
GH_A08G1119	GhIQDA08.1	445	1 338	49 676.63	9.38	E.R.	4
GH_D08G1102	GhIQDD08.1	445	1 338	49 457.39	9.38	E.R.	4
GH_A09G1181	GhIQDA09.1	436	1 311	47 752.06	10	chlo	6
GH_A08G2051	GhIQDA08.3	317	954	35 369.75	10.06	chlo	4
GH_A09G1949	GhIQDA09.4	461	1 386	51 179.31	10.33	chlo	4
GH_A10G2676	GhIQDA10.3	384	1 155	43 186.39	9.93	chlo	4
GH_A13G0121	GhIQDA13.1	385	1 158	44 219.64	10.37	chlo	5
GH_D08G1176	GhIQDD08.2	456	1 371	50 866.22	10.21	chlo	4
GH_D08G2064	GhIQDD08.4	508	1 527	56 482.43	9.83	chlo	6
GH_D10G2782	GhIQDD10.3	398	1 197	45 120.65	9.88	chlo	4
GH_D13G0120	GhIQDD13.1	385	1 158	44 188.66	10.27	chlo	5

1. nucl indicates nucleus, mito indicates mitochondria, chlo indicates chloroplast, E.R. indicates endoplasmic reticulum

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. hirsutumL., 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. hirsutumL., 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 [less than] 0.05 level, ** represents an extremely significant difference at the p [less than] 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 (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-ΔΔCT 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.

Availability of data and materials

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

References

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. https://doi.org/10.1186/1471-2148-5-72 .

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. https://doi.org/10.3389/fpls.2018.01669 .

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. https://doi.org/10.1007/s00425-018-2950-6 .

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. https://doi.org/10.1080/15592324.2017.1331198 .

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. https://doi.org/10.1104/pp.16.01743 .

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. https://doi.org/10.1074/jbc.M112.396200 .

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. https://doi.org/10.1038/srep44743 .

DeFalco TA, Bender KW, Snedden WA. Breaking the code: Ca2+ sensors in plant signalling. Biochem J. 2009;425:27–40. https://doi.org/10.1042/BJ20091147 .

Dong C, Wang J, Yu Y. Identifying functional genes influencing Gossypium hirsutum fiber quality. Front Plant Sci. 2018;9:1968. https://doi.org/10.3389/fpls.2018.01968 .

Feng L, Chen Z, Ma H, et al. The IQD gene family in soybean: structure, phylogeny, evolution and expression. PLoS One. 2014;9:e110896. https://doi.org/10.1371/journal.pone.0110896 .

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. https://doi.org/10.1186/s12870-019-1968-7 .

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. https://doi.org/10.1016/j.abb.2004.04.025 .

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. https://doi.org/10.1111/nph.12381 .

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. https://doi.org/10.1038/s41588-019-0371-5 .

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. https://doi.org/10.1007/s00438-013-0733-0 .

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. https://doi.org/10.1186/s12870-019-1947-z .

Knip M. The SLEEPER genes: a transposase-derived angiosperm-specific gene family. BMC Plant Biol. 2012;12:192. https://doi.org/10.1186/1471-2229-12-192 .

Krzywinski M, Schein J, Biro I, et al. Circos: an information aesthetic for comparative genomics. Genome Res. 2009;19:1639–45. https://doi.org/10.1101/gr.092759.109 .

Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33:1870. https://doi.org/10.1093/molbev/msw054 .

Li JH, Wu YK, Zhang Q. Aphid parasitism and parasitoid diversity in cotton fields in Xinjiang, China. PLoS One. 2018;13:e0207034. https://doi.org/10.1371/journal.pone.0207034 .

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. https://doi.org/10.1104/pp.18.00558 .

Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2. https://doi.org/10.1093/bioinformatics/btp187 .

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. https://doi.org/10.1007/s11104-019-04140-8 .

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. https://doi.org/10.1006/meth.2001.1262 .

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. https://doi.org/10.1016/j.plantsci.2014.08.017 .

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. https://doi.org/10.2741/4500 .

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. https://doi.org/10.3390/ijms19041038 .

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. https://doi.org/10.1105/tpc.111.084988 .

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. https://doi.org/10.1007/s12104-019-09883-0 .

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. https://doi.org/10.1186/s12864-018-4985-2 .

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. https://doi.org/10.1016/j.cj.2019.11.003 .

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. https://doi.org/10.4161/psb.24369 .

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. https://doi.org/10.1186/s12864-018-5047-5 .

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. https://doi.org/10.1038/srep42418 .

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. https://doi.org/10.1186/s12864-015-1708-9 .

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. https://doi.org/10.1371/journal.pone.0176399 .

Voorrips RE. MapChart: software for the graphical presentation of linkage maps and QTLs. J Hered. 2002;93:77–8. https://doi.org/10.1093/jhered/93.1.77 .

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. https://doi.org/10.1186/s12870-018-1526-8 .

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. https://doi.org/10.1038/s41588-018-0282-x.

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. https://doi.org/10.1016/j.gene.2013.07.084 .

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. https://doi.org/10.1093/nar/gkr1293 .

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. https://doi.org/10.1080/07391102.2019.1633409 .

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. https://doi.org/10.1038/srep24520 .

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. https://doi.org/10.1016/j.molcel.2018.03.013 .

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. https://doi.org/10.3389/fpls.2017.02272 .

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. https://doi.org/10.1371/journal.pone.0122027 .

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. https://doi.org/10.1016/j.pbi.2017.06.003 .

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. https://doi.org/10.1007/s10658-018-1591-8 .

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. https://doi.org/10.1016/j.cj.2019.03.008 .

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. https://doi.org/10.3389/fpls.2018.00721 .

Zhu T, Liang C, Meng Z, et al. CottonFGD: an integrated functional genomics database for cotton. BMC Plant Biol. 2017;17:101. https://doi.org/10.1186/s12870-017-1039-x .