Characterization and Stress Response of the JmjC Domain-Containing Histone Demethylase Gene Family in the Allotetraploid Cotton Species Gossypium hirsutum

Histone modification is an important epigenetic modification that controls gene transcriptional regulation in eukaryotes. Histone methylation is accomplished by histone methyltransferase and can occur on two amino acid residues, arginine and lysine. JumonjiC (JmjC) domain-containing histone demethylase regulates gene transcription and chromatin structure by changing the methylation state of the lysine residue site and plays an important role in plant growth and development. In this study, we carried out genome-wide identification and comprehensive analysis of JmjC genes in the allotetraploid cotton species Gossypium hirsutum. In total, 50 JmjC genes were identified and in G. hirsutum, and 25 JmjC genes were identified in its two diploid progenitors, G. arboreum and G. raimondii, respectively. Phylogenetic analysis divided these JmjC genes into five subfamilies. A collinearity analysis of the two subgenomes of G. hirsutum and the genomes of G. arboreum and G. raimondii uncovered a one-to-one relationship between homologous genes of the JmjC gene family. Most homologs in the JmjC gene family between A and D subgenomes of G. hirsutum have similar exon-intron structures, which indicated that JmjC family genes were conserved after the polyploidization. All G. hirsutum JmjC genes were found to have a typical JmjC domain, and some genes also possess other special domains important for their function. Analysis of promoter regions revealed that cis-acting elements, such as those related to hormone and abiotic stress response, were enriched in G. hirsutum JmjC genes. According to a reverse transcription-quantitative polymerase chain reaction (RT-qPCR) analysis, most G. hirsutum JmjC genes had high abundance expression at developmental stages of fibers, suggesting that they might participate in cotton fiber development. In addition, some G. hirsutum JmjC genes were found to have different degrees of response to cold or osmotic stress, thus indicating their potential role in these types of abiotic stress response. Our results provide useful information for understanding the evolutionary history and biological function of JmjC genes in cotton.


Introduction
Epigenetics is the study of changes in gene expression that can be inherited without changes in the DNA sequence [1]. Gene expression is controlled through epigenetic regulation, which mainly includes DNA methylation, histone modification, chromatin remodeling and non-coding RNA regulation [2]. Eukaryotic genomic DNA is tightly packed into chromatin. Nucleosomes, which are the basic structural units of chromatin, consist of 146 bp of DNA wrapped around a histone octamer core composed of two molecules each of H2A, H2B, H3 and H4 [3]. Post-translational covalent modification of histones is an important mode of epigenetic regulation that regulates gene expression and other processes by affecting the state of chromatin [4]. The N-terminal amino acid tails of histones are often covalently modified, mainly through acetylation, methylation, ubiquitination, phosphorylation, glycosylation, ADP-ribosylation and SUMOylation [5]. Different covalent modifications occur at different sites of histones. For example, lysine can be acetylated, methylated and ubiquitinated, arginine can be methylated, and threonine and serine can be phosphorylated [6]. Histone methylation has important effects on gene transcription and chromatin structure that vary according to the location and degree of methylation. Histone methylation occurs mainly at the N-terminal arginine (R) and lysine (K) sites of histones H3 and H4, including K4, K9, K27, K36 of H3 and K20 for H4. Lysine residues can occur in three methylation forms, namely, monomethylated (Kme1), dimethylated (Kme2), or trimethylated (Kme3), whereas arginine can only undergo monomethylation (Rme1), symmetric dimethylation (Rme2s), or asymmetric dimethylation (Rme2a) [7]. Histone lysine methylation has several roles, not only affecting the static charge of modified residues, but also improving hydrophobicity and changing intramolecular or intermolecular interactions [5]. In general, methylation of histones H3K9 (H3K9me2/3) and H3K27 (H3K27me3) is associated with transcriptional inhibition, whereas methylation of H3K4 (H3K4me2/3) and H3K36 (H3K36me3) is associated with transcriptional activation [8]. Two histone lysine demethylases are currently known. The first, histone lysine-specific demethylase 1 (LSD1), reverses the monomethylation or dimethylation of lysine residues by flavin adenine dinucleotide (FAD), but does not act on trimethylated lysine [9]. Four proteins homologous to LSD1, namely, FLOWERING LOCUS D (FLD), LSD1-LIKE 1 (LDL1), LSD1-LIKE 2 (LDL2) and LSD1-LIKE 3 (LDL3), have been found in Arabidopsis and rice [10,11]. The second histone lysine demethylase, JumonjiC (JmjC) domain-containing histone demethylase, is functionally dependent on ferrous ion (Fe [II]) and α-ketoglutarate (α-KG) as prosthetic groups and is active on all three methylated lysines [12]. The first JmjC protein, which is essential for heart and brain development [13], was originally identified in mice and subsequently found to have demethylase activity [14]. More and more evidence suggested that the JmjC proteins constitute an important family of histone lysine demethylases and play an important role in maintaining stable histone methylation in vivo [9,15].
Moreover, the JmjC genes also play crucial roles in plant response to stresses. The overexpressed plants of AtJMJ15 showed greater salt tolerance compared with the wild-type Arabidopsis, while the functionally deficient mutant was more sensitive to salt stress [31]. AtJMJ17 participates in Arabidopsis response to dehydration stress and abscisic acid (ABA), and mutants with an inactivated form of this gene show dehydration stress tolerance and ABA hypersensitivity during stomatal closure [32]. In Medicago truncatula, MtJMJC5 undergoes cold specifically induced alternative splicing, and the cold-dependent alternative splicing could be reversed after ambient temperature returning to the normal [33]. In Brassica rapa, the transcripts of JmjC genes were induced under hot and cold treatments [34]. In maize, all 19 JmjC genes were responsive to heat stress treatment [35]. In rice, OsJMJ704 epigenetically suppressed defense negative regulator genes [36], while OsJMJ705 epigenetically activated defense positive regulator genes [37], to enhance resistance to pathogen attack.
Gossypium constitutes approximately 50 species, including 45 diploid and 5 tetraploid species, with sporophytic chromosome counts of 26 (2n = 2x = 26) and 52 (2n = 4x = 52), respectively. On the basis of cytological studies, eight different genomes of diploid Gossypium species have been distinguished (A-G and K), while the five tetraploid species are allotetraploids possessing both A and D diploid genomes [38]. Two of the diploids, the wild species Gossypium raimondii and the cultivated cotton G. arboreum, harbor D and A genomes, respectively [39,40]. Approximately 1.5 million years ago, diploid cotton (G. arboreum) originating from Africa and diploid cotton (G. raimondii) from South America gave rise via hybridization and genome doubling to five allotetraploid species, namely, G. hirsutum, G. barbadense, G. tomentosum, G. darwinii and G. mustelinum [38,39]. Among them, G. hirsutum was domesticated and is now the main cultivated cotton species [39,41]. More than 95% of commercially obtained cotton fiber is from G. hirsutum, which is characterized by its high productivity and strong environmental adaptability [42,43]. Moreover, because cotton includes a wide variety of species with rich genetic diversity, it has important applications in the exploration of the origin, genomic evolution and polyploidy of plants [44,45].
Cotton is an economically important crop as well as a model system to studying polyploidization [45,46]. As one of the most important fiber crops, the allotetraploid cotton species, G. hirsutum, is the source of much of the natural fiber used in the textile industry worldwide [47,48]. However, the role and function of JmjC domain-containing histone demethylase genes in cotton is poorly understood, so requires further investigation. In this study, consequently, we identified the JmjC genes at the genome-wide level in G. hirsutum, and examined their phylogenetic relationships, chromosomal locations, collinearity, gene and domain structures, cis-elements in the promoters and spatiotemporal expressions. Our results provide insights into the evolution and function of JmjC genes in G. hirsutum.

Identification of G. hirsutum JmjC Genes
Using a combination of methods, we identified 50 JmjC genes in the allotetraploid cotton species G. hirsutum, more than twice the number found in Arabidopsis (21) or rice (20). To better understand the evolutionary history of G. hirsutum JmjC genes, we searched for JmjC genes in the other two diploid cotton species, G. arboreum and G. raimondii, using the same method. We identified 25 G. arboreum and 25 G. raimondii JmjC genes, exactly half the number found in G. hirsutum. The newly identified JmjC genes in G. raimondii were first named GrJMJ1-25 according to their order on chromosomes. The JmjC genes in G. arboreum were named GaJMJs based on the homologs in G. raimondii with the same number as in G. raimondii. The JmjC genes in G. hirsutum were named GhJMJs, corresponding to the homologs in G. raimondii and G. arboreum, with suffixes D and A added to represent D and A subgenomes, respectively. Basic features of G. hirsutum JmjC genes, including chromosomal location, protein length, molecular weight (Mw), isoelectric point (pI), and predicted subcellular localization, are listed in Table S1. An analysis of physicochemical characteristics revealed that amino acid number, Mw, and theoretical pI vary widely among JmjC gene family members in G. hirsutum. According to the results of subcellular prediction, most JmjC gene family members are located in the nucleus, which is consistent with their histone demethylation function.

Phylogenetic Analysis of G. hirsutum JmjC Genes
To assess phylogenetic relationships of JmjC family genes, we performed multiple alignments of full-length protein sequences and JmjC domain sequences of JmjC genes from G. hirsutum, G. arboreum, G. raimondii, Arabidopsis and rice (Table S2). The evolutionary trees were constructed from the aligned sequences using MEGA 7.0 software, respectively. As shown in Figure 1 and Figure
Gossypium hirsutum is an allotetraploid cotton species with AA and DD genomes whose progenitors are widely believed to be the extant diploid cotton species G. arboreum and G. raimondii [40,43]. To explore the evolutionary history of the JmjC gene family in Gossypium, we performed synteny analysis to identify directly homologous JmjC genes in G. hirsutum, G. arboreum and G. raimondii. The locations of each JmjC gene and homologous gene pair uncovered in the analysis are shown in Figure 2. We identified 25 directly homologous gene pairs between the D subgenome of G. hirsutum and G. raimondii, 25 directly homologous gene pairs between the A subgenome of G. hirsutum and G. arboreum, and 25 directly homologous gene pairs between the D and A subgenomes of G. hirsutum. According to these results, a one-to-one relationship exists between homologs of the JmjC gene family in the two diploid cotton species and the two subgenomes of tetraploid cotton, which indicates that JmjC family genes have been highly conserved during cotton evolution.

Structural Features and Conserved Domains of G. hirsutum JmjC Genes
Exon-intron structures play a crucial role in the evolution of organisms [49]. To understand the structural diversity of the JmjC gene family in G. hirsutum, we further analyzed exon-intron structures of these genes using the online Gene Structure Display Server v2.0 ( Figure 3A,B). The number of exons in G. hirsutum JmjC genes varied from 2 to 33. Specifically, all genes possessed 7 to 16 exons except for GhJMJ20A/GhJMJ20D with 2 and GhJMJ6A/GhJMJ6D with 33. Most homologs in the JmjC gene family between A and D subgenomes of G. hirsutum have similar gene structures. For example, such gene pairs have highly consistent gene structures, exons/introns that are consistent in arrangement and number, and exons whose lengths are very similar. Exceptions exist, however, such as four gene pairs (GhJMJ11A/GhJMJ11D, GhJMJ15A/GhJMJ15D, GhJMJ18A/GhJMJ18D and GhJMJ25A/GhJMJ25D) have different intron numbers. Conserved domains refer to domains that are invariant or identical over biological evolution or within a protein family and that generally have important functions that cannot be changed [50]. In histone demethylases containing the highly conserved JmjC domain, the activity of these proteins can be predicted by focusing on the conserved Fe (II) and α-KG binding sites within the domain [16].
To gain a better understanding of domain diversity, we determined the domain architectures of 50 G. hirsutum JmjC proteins by comparison with the full-length protein sequences of JmjC genes in the Pfam and SMART databases. As shown in Figure 3A,C all members of the G. hirsutum JmjC gene family were found to have a conserved JmjC domain, which is associated in JmjC proteins with demethylation of the histone lysine site [15]. The second most widely distributed domain was the JmjN domain, which appeared in all members of groups KDM4/JHDM3 and KDM5/JARID. All members of the KDM5/JARID subfamily also contained a zf-C5HC2 domain, which in combination with genes downstream of REF6 (AtJMJ12) allows the latter to assume the role of a histone demethylase [51]. The zf-C5HC2 domain was not observed in some members of subfamilies KDM4/JHDM3, possibly because the zf-C5HC2 domain in these proteins was less conserved and difficult to detect under the analysis parameters of this study. In addition to GhJMJ2A/GhJMJ2D and GhJMJ6A/GhJMJ6D, other members of the KDM5/JARID subfamily possess the FYRC domain that may have chromatin binding activity or contribute to the function of the JmjC domain by interacting with other proteins [16]. It is also found that GhJMJ6A and GhJMJ6D have three unique structural domains, the ARID domain, PHD domain and PLU-1 domain. The ARID domain is related to sequence-specific DNA binding [52], the PHD domain specifically recognize methylated (modified) histone codes and may be important readers of histone codes [53], and the PLU-1 domain plays a critical role in gene regulation and chromatin stability [54]. In the JMJD6 subfamily, four JmjC genes (GhJMJ3A/GhJMJ3D and GhJMJ18A/GhJMJ18D) containing the F-box domain are consistent with previous study [55].

Cis-Elements in the Promoter Regions of G. hirsutum JmjC Genes
A promoter region is a DNA sequence located upstream of the 5 -end of a gene that activates RNA polymerase to facilitate accurate binding to template DNA and the specificity of transcription initiation [56]. To further elucidate the possible regulatory mechanisms of G. hirsutum JmjC genes in abiotic stress response, we searched the PlantCARE database for cis-regulatory elements in JmjC promoter regions (Figure 4).
We detected 12 types of cis-regulatory elements related to hormones and stresses in the promoters of G. hirsutum JmjC genes. Hormone-related elements included ABREs (abscisic acid responsiveness), EREs (ethylene responsive element), CGTCA motifs (MeJA-responsiveness), AUXRR-core motifs (auxin-responsive element) and GAREs (gibberellin-responsive element). Plant hormones can affect plant growth and induce physiological responses to environmental stimuli [57,58]. ABREs, which are cis-acting elements involved in abscisic acid reactivity, promote stomatal closure to control water loss in dry conditions [48,59]. The GARE motif is a responsive element for gibberellin, which plays an important role during seed, reproductive and maturation stages [60]. The CGTCA motif, related to jasmonic acid, regulates plant growth and is involved in the development of germination, rooting, flowering, fruiting and senescence. This element also actively participates in mechanisms underlying defense against biotic and abiotic stress [61]. Other cis-acting promoter elements including LTR (low-temperature responsiveness), MBS (myeloblastosis binding sites) and Myb (v-myb avian myeloblastosis viral oncogene homolog) are pressure-dependent. For example, the Myb transcription factor family can improve plant resistance drought stress [62]. We also detected cis-acting elements known to respond to biological or abiotic stresses including G-box and GT1-motif (termed the glycosyltransferase 1 motif), circadian (cis-regulated elements involved in circadian rhythm) and AREs (anaerobic-induced essential regulated elements) in promoter regions of G. hirsutum JmjC genes. ABRE, cis-acting element involved in abscisic acid responsiveness; ARE, anaerobic-induced essential regulatory element; AUXRR-core, an auxin-related element; CGTCA-motif, cis-acting regulatory element involved in MeJA responsiveness; circadian, involved in circadian cis-regulating elements; ERE, ethylene-responsive element; GARE-motif, gibberellin-responsive element; G-box, cis-acting element with light effect; GT1-motif, involved in salt-stress response elements; LTR, cis-acting element involved in hypothermic stress response; MBS, related to drought induction; Myb, cis-acting element involved in response to drought, high temperature and low temperature.

Expression Profiling of G. hirsutum JmjC Genes in Different Tissues and under Abiotic Stresses
To explore JmjC gene expression patterns during cotton development, we used reverse transcription-quantitative polymerase chain reaction (RT-qPCR) to measure JmjC gene expressions in roots, stems, leaves, ovules (0, 1, 5, 10 and 20 DPA) and fibers (10, 15, 20 and 30 DPA) ( Figure 5A,B and Figure S3). As shown in the heat map revealing the expression patterns of different G. hirsutum JmjC genes, the majority of the homologous genes in the JmjC family have similar expression patterns. Most JmjC family members had their highest expression in 20 and 30 DPA fibers. Some family members had higher expression levels in 10 and 15 DPA fibers (GhJMJ5D, GhJMJ8A, GhJMJ13A, GhJMJ17D, GhJMJ18D and GhJMJ24A). A few genes were also highly expressed in 0 and 1 DPA ovules (GhJMJ1D, GhJMJ10D, GhJMJ11D, GhJMJ13D, GhJMJ14D, GhJMJ15A, GhJMJ15D and GhJMJ19D). In contrast, expression levels in root, stem and leaf tissues were generally low. According to the RT-qPCR results, G. hirsutum JmjC genes may play a role in the development of G. hirsutum fibers. To further confirm the response of G. hirsutum JmjC genes to abiotic stress, we carried out a RT-qPCR analysis on all genes in the G. hirsutum JmjC family using cotton subjected to cold and polyethylene glycol (PEG) treatments ( Figure 6A,B; Figures S4 and S5). According to the expression profiling, all genes were differentially expressed in leaves in response to cold stress (4 • C) or osmotic stress (PEG) compared with normal conditions (0 h). Most members of the G. hirsutum JmjC gene family were up-regulated in response to cold stress or osmotic stress. GhJMJ1A, GhJMJ3A, GhJMJ4A, GhJMJ6A, GhJMJ14A, GhJMJ14D, GhJMJ18A, GhJMJ18D, GhJMJ19D, GhJMJ22A, GhJMJ22D, GhJMJ24A and GhJMJ24D were significantly up-regulated under cold stress compared with non-treated controls. Under osmotic stress conditions of different durations, GhJMJ3A, GhJMJ5D, GhJMJ9A, GhJMJ9D, GhJMJ18A, GhJMJ18D, GhJMJ19D, GhJMJ22D, GhJMJ24A and GhJMJ24D were obviously up-regulated. Seven G. hirsutum JmjC genes (GhJMJ3A, GhJMJ18A, GhJMJ18D, GhJMJ19D, GhJMJ22D, GhJMJ24A and GhJMJ24D) were thus obviously up-regulated both under cold and osmotic stress conditions. These results suggest that these G. hirsutum JmjC genes respond to cold and osmotic stresses and have potential roles during cold and drought conditions.

Discussion
In plants, histone methylation plays a crucial role in a variety of biological processes, ranging from transcriptional regulation to heterochromatin formation. Different levels of methylation are responsible for different biological functions and are, therefore, essential for plant growth and development [5]. Histone demethylases containing the JmjC domain, which is highly conserved in plants, exert a key role in maintaining homeostasis of histone methylation in vivo [5,63]. Proteins containing the JmjC domain constitute an important family of histone lysine demethylases in animals and plants and play an important role in histone modification, which is an essential aspect of epigenetics [64]. At the genome-wide level, the JmjC gene family has been analyzed to account for its evolutionary history and biological function in certain plant species. For example, the evolution of histone demethylase has been studied in some green plant lineages [65], and 21 members in Arabidopsis [66], 48 members in soybean [67], 22 members in woodland strawberry [68], 28 members in Rosa chinensis [69], 20 members in rice [70] and 19 members in maize [35] have been identified and analyzed for the JmjC gene family, respectively. In contrast, no systematic study has been carried out on the JmjC gene family in G. hirsutum.
In the present study, we conducted a comprehensive analysis of G. hirsutum JmjC genes, including their phylogenetic relationships, gene structures, domain structures, chromosomal locations, collinearity, promoter elements and expression profiles.
In this study, we identified 50 JmjC genes in the G. hirsutum genome for the first time and also detected 25 genes each in G. arboreum and G. raimondii, which is higher than the number found in Arabidopsis (21) and rice (20). Phylogenetic analysis of JmjC genes from five plant species divided cotton JmjC genes into five main subclasses. The KDM3/JHDM2 subfamily was found to contain the largest number of genes, accounting for 36% of the total. We also analyzed the structure of G. hirsutum JmjC proteins and thereby uncovered some preliminary details of their structural characteristics. Members of the KDM5/JARID subfamily, except for GhJMJ2A, GhJMJ2D, GhJMJ6A GhJMJ6D and GhJMJ21A, have the same domain, and they also share FYRC domains. AtJMJ14 has been shown to play a role in Arabidopsis flowering time control through the interaction of FYRC domains with transcription factors [27]. The G. hirsutum JmjC genes GhJMJ6A and GhJMJ6D in the KDM5/JARID subfamily also contain three domains (ARID, PHD and PLU-1) that are related to DNA binding and gene regulation [53,54,71]. In Arabidopsis the homologous gene of GhJMJ6A and GhJMJ6D is AtJMJ17, which binds directly to the chromatin of OPEN STOMATA 1 (OST1) and demethylated H3K4me3 to regulate OST1 mRNA abundance, thereby regulating response to dehydration stress [32]. GhJMJ6A and GhJMJ6D may, therefore, have similar functions in cotton.
Gossypium hirsutum, as an AADD allotetraploid species, evolved from hybridization and genome doubling of A genome diploid species (G. arboreum) and D genome diploid species (G. raimondii) [40,43]. Through analysis of the locations and homologous relationships of JmjC genes in G. raimondii, G. arboreum, and A and D subgenomes of G. hirsutum, the 50 JmjC genes in G. hirsutum A and D subgenomes have a one-to-one direct homology with the 25 JmjC genes located in G. raimondii and G. arboreum, which indicates that these 50 genes are derived from the common ancestor of G. raimondii and G. arboreum. The number of JmjC genes in G. hirsutum is twice that of G. raimondii and G. arboreum, which indicates that all JmjC genes were retained during the process of polyploidization. Studies have shown that the JmjC gene family has diverse functions, including involvement in growth and development pathways, such as those related to flowering, photoperiod, hormone, and cell division and differentiation, and can reverse methylation at multiple sites, such as H3K4, H3K9, and H3K27. Moreover, some JmjC demethylases can remove methylation from multiple different sites [12]. The expansion of the number of JmjC genes may thus be related to the special characteristics of cotton or its adaptation to the environment.
To further study the possible functions of JmjC genes during G. hirsutum growth and development, we used RT-qPCR to reveal the spatio-temporal transcription patterns of these genes in 12 tissues at different stages of G. hirsutum development. Cotton fibers are the longest single cells found in plants.
Many studies have shown that the morphological development of cotton fibers takes place in four stages: fiber initiation (on or near the day of anthesis), cell elongation (from the day of anthesis to approximately 21 to 26 DPA), secondary wall deposition and maturation (At approximately 45 to 60 DPA, the seed capsule dehisces and the thin fiber cells quickly dehydrate). Each of these stages is critical to fiber development [72]. According to the RT-qPCR results, most G. hirsutum JmjC genes clearly had similar expression patterns, that is, higher expression levels during fiber development occurring at 10-30 DPA, especially 20 and 30 DPA. In Arabidopsis, mutations of the REF6 (AtJMJ12) gene cause impaired cell elongation, and REF6 (AtJMJ12) and ELF6 (AtJMJ11) can also interact with the transcription factor BRI1 EMS SUPPRESSOR 1 (BES1) in the brassinosteroid (BR) signaling pathway to regulate the expression of BR reaction genes [24]. In addition, REF6 (AtJMJ12) promotes lateral root formation by eliminating the inhibition of PIN1/3/7 genes [73]. The homologs of REF6 (AtJMJ12) and ELF6 (AtJMJ11) in G. hirsutum, GhJMJ10A and GhJMJ19A, are highly expressed in 20 and 30 DPA fibers, which suggests that these G. hirsutum JmjC genes are involved in cotton fiber development.
Cis-acting elements in gene promoter regions play an important role in gene transcription initiation and participate in the regulation of gene expression [56]. We analyzed promoter sequences of G.
hirsutum JmjC genes and uncovered important cis-acting elements, such as those related to hormone and abiotic stress responses. We thus speculated that the G. hirsutum JmjC family genes have functions in plant response to abiotic stress. To confirm this hypothesis, we used RT-qPCR to analyze expression patterns of JmjC genes in G. hirsutum plants subjected to cold and osmotic treatments. In Arabidopsis, AtJMJ17 affects the expression of such genes by regulating the H3K4me3 level of the stress response genes encoding stomatal opening factor (OST1) and ABA INSENSITIVE 5 (ABI5) and plays a negative regulatory role in osmotic stress and ABA response [32]. According to the expression analysis in this study, the AtJMJ17 gene homologs GhJMJ6A/GhJMJ6D were highly expressed under both cold and osmotic stresses. At the same time, the cis-acting element analysis shows that GhJMJ6A/GhJMJ6D also contains ABRE cis-acting element. We thus speculate that these genes are closely related to the response of G. hirsutum to cold or osmotic stress. Our RT-qPCR analysis of these differentially expressed genes has thus provided novel insights into plant responses to cold and osmotic stresses.

Identification of JmjC Genes in G. hirsutum and other Cotton Species
Gossypium hirsutum, G. arboreum and G. raimondii genomic data were downloaded from the CottonFGD database [74] and used to create a local BLAST database. To identify JmjC genes in cotton, BLASTP and TBLASTN searches were carried out using the 21 Arabidopsis JmjC protein sequences [66] as queries. All identified candidate genes were analyzed online using Pfam [75] and SMART [76] databases to check for JmjC domains (PF02373).
Physicochemical parameters of G. hirsutum JmjC genes, including molecular weight (Mw) and isoelectric point (pI), were calculated using the compute pI/Mw tool in EXPASy [77], with the parameter set to average. Subcellular localization predictions of G. hirsutum JmjC proteins were carried out online using CELLO v2.5 [78].

Phylogenetic Analysis of G. hirsutum JmjC Genes
The amino acid sequences were aligned using the ClustalW program [79] with default parameters. According to previous published studies [35,67,69], the phylogenetic tree was then constructed using the neighbor-joining method as implemented in MEGA 7.0 [80] with the pairwise deletion option. Bootstrapping based on the p-distance model and 1000 replicates was conducted to assess internal branch reliability.

Chromosomal Mapping and Synteny Analysis of G. hirsutum JmjC Genes
The physical location data of JmjC genes were retrieved from the Gossypium genomic databases. The MapInspect software was used to generate chromosome distribution images [81]. The synteny analysis was performed as previously reported [82]. The chromosomal distribution and collinearity of G. arboreum, G. raimondii and G. hirsutum JmjC genes were then illustrated using the Advanced Circos circle diagram function in TBtools [83].

Structural Analysis of G. hirsutum JmjC Genes
To construct exon-intron structures that included the location and length of each G. hirsutum JmjC gene, coding sequences were compared with corresponding genomic sequences using the online Gene Structure Display Server v2.0 [84]. Conserved domains in JmjC domain-containing proteins encoded by G. hirsutum JmjC genes were identified through searches of Pfam [75] and SMART [76] databases with default parameters.

Prediction of Cis-Acting Elements in Promoter Regions of G. hirsutum JmjC Genes
To identify potential stress-response and hormone-related cis-acting elements in promoter sequences of G. hirsutum JmjC genes, the 2-kb sequence upstream of the initiation codon (ATG) of each G. hirsutum JmjC gene was compared against the PlantCARE database [85]. Adobe Illustrator software was used for the final figure construction.

Plant Materials and Stress Treatments
Gossypium hirsutum acc. TM-1 were used in the study. Roots, stems and leaves were collected from two-week-old seedlings, ovules were collected at 0, 1, 5, 10 and 20 DPA, and fibers were collected at 10, 15, 20 and 30 DPA. Cotton seedlings were planted in a liquid culture medium and grown until the third true leaf appeared in a plant growth chamber at 28 • C under a 16/8 h light/dark photoperiod. For the cold treatment, seedlings were incubated at 4 • C. For PEG treatment, seedlings were irrigated with 20% PEG 6000 solution. After cold and PEG treatments, the leaves were collected at 0, 1, 3, 6 and 12 h. All samples immediately were frozen in liquid nitrogen and stored at −80 • C.

Expression Analysis of G. hirsutum JmjC Genes
The RT-qPCR analysis was performed to reveal expression patterns of JmjC genes in representative tissues, namely, roots, stems, leaves, ovules (0, 1, 5, 10 and 20 DPA), and fibers (10,15,20 and 30 DPA). We also used the RT-qPCR analysis to examine expression patterns of JmjC genes in leaves of plants subjected to cold or PEG stress for 0, 1, 3, 6 and 12 h. Total RNA was extracted from different samples using an RN533-EASYspin Plus RNA rapid extraction kit (Adlai, Beijing, China). First-strand cDNA synthesis was carried out using 1 µg of total RNA and a PrimeScript first-strand cDNA synthesis kit (Takara, Dalian, China). Gene-specific primers (Table S3) were designed based on JmjC nucleotide sequences with Oligo software (v7.60, Molecular Biology Insights, Cascade, CO, USA). Reverse-transcription PCR was performed using the following protocol: 94 • C for 5 min; followed by 35 cycles of 94 • C for 30 s, 60 • C for 30 s, and 72 • C for 1 min; and the 72 • C for 10 min. The amplified fragment was purified with a FastPure Gel DNA Extraction mini kit (Vazyme, Nanjing, China) and cloned into a PMD18-T Vector (Takara, Dalian, China) for sequence verification. RT-qPCR was carried out on a RocheLightCycler480 instrument (Roche, Basel, Switzerland) the amplification parameters were as follows: stage (1): 95 • C, 5 min; stage (2): 40 cycles of 95 • C for 10 s, 60 • C for 10 s, 72 • C for 10 s; stage (3): extension at 72 • C for 10 min. The GhHIS3 gene was used as an internal control, with three biological replicates per sample. The 2-∆CT was used to calculate the relative mRNA levels of each gene in different tissues. The 2−∆∆CT was used to calculate the relative mRNA level of each gene in leaves under the cold and osmotic stress. A heatmap was generated using the HeatMap IIIustrator function in TBtools [83].

Conclusions
In this study, we conducted the first-ever comprehensive analysis of the G. hirsutum JmjC gene family. Using bioinformatics methods, we carried out a detailed investigation of phylogenetic relationships, chromosomal locations, collinear relationships, exon-intron structures, conserved domains and promoter elements of G. hirsutum JmjC family members. We also analyzed their expression patterns in different tissues and under cold and osmotic stresses.
The 50 G. hirsutum JmjC genes identified in this study were divided into five groups based on the result of phylogenetic analysis. The chromosomal localization showed that these JmjC genes are distributed at different densities on 23 of the 26 chromosomes of G. hirsutum. A collinearity analysis revealed that a one-to-one homologous relationship exists between JmjC genes in the two diploid cotton species G. arboreum and G. raimondii or in the two subgenomes of G. hirsutum. Most homologs in the JmjC gene family between A and D subgenomes of G. hirsutum have similar exon-intron structures, suggesting that the JmjC genes are strongly conserved in G. hirsutum. All G. hirsutum JmjC genes possess a typical JmjC domain, and other special domains were also found in some genes. According to an analysis of cis-acting promoter elements, G. hirsutum JmjC genes contain elements responsive to hormonal or abiotic stresses and might play an important role in plant growth and development.
A RT-qPCR analysis demonstrated that G. hirsutum JmjC genes might have an important role in the development of fibers. This analysis also revealed differing responses of G. hirsutum JmjC genes to cold and osmotic stresses. Most family members exhibited significantly higher expression levels under cold or osmotic stress, which indicates that these genes are closely related to G. hirsutum response to cold or osmotic stress.
These results provide valuable clues for future identification of the specific functions of this gene family and the genetic diversity of cotton and other dicotyledonous plants.
Supplementary Materials: The following are available online at http://www.mdpi.com/2223-7747/9/11/1617/s1, Table S1. Basic information on JmjC genes in G. hirsutum. Table S2. Protein sequences of JmjC genes from G. hirsutum, G. raimondii, G. arboreum, Arabidopsis and rice. Table S3. Primers used in this study. Figure S1. Phylogenetic tree of JmjC domain-containing histone demethylases from G. hirsutum (Gh), G. arboreum (Ga), G. raimondii (Gr), Arabidopsis (At) and rice (Os). Phylogenetic analysis was performed in MEGA 7.0 using the neighbor-joining method based on JmjC domain sequences of JmjC genes from different plant species. Figure S2. Physical locations of JmjC genes in the G. hirsutum genome. Figure S3. Expression profiles of G. hirsutum JmjC genes in different tissues. The error bars indicate the standard deviation estimated from three biological replicates. Figure S4. Expression profiles of G. hirsutum JmjC genes under cold treatments. The error bars indicate the standard deviation estimated from three biological replicates. Figure S5. Expression profiles of G. hirsutum JmjC genes under PEG treatments. The error bars indicate the standard deviation estimated from three biological replicates.

Conflicts of Interest:
The authors declare no conflict of interest.