To identify members of the MIKC gene family in G. hirsutum, Arabidopsis MIKC sequences were obtained from the TAIR database (http://www.arabidopsis.org) and used as queries for a BLASTP algorithm-based against the G. hirsutum genome database (https://www.cottongen.org/species/Gossypium_hirsutum/nbi-AD1_genome_v1.1) (Zhang et al., 2015). The MIKC protein domain was analyzed using the Hidden Markov Model (HMM) from the Pfam database (http://pfam.xfam.org/). The SRF-TF and K-box domains were confirmed by Pfam accessions (PF00319 and PF01486, respectively). All of the candidate proteins were manually checked using the above described methods to remove the redundant sequences.

To construct a MIKC-protein phylogenetic tree using MEGA 6.06, MIKC proteins from four plant species, Arabidopsis, O. sativa, V. vinifera, and G. hirsutum, were employed. The neighbor-joining method with amino acid p-distance was applied to construct the tree (Tamura et al., 2011), and the reliability was obtained by bootstrapping with 1,000 replicates.

The exon/intron structures of MIKC genes were retrieved by the alignment of predicted coding sequences with corresponding genomic sequences using the gene structure display server (GSDS) program (http://gsds.cbi.pku.edu.cn/).

The online program MEME (http://meme-suite.org/) was employed to determine the conserved motifs in GhMIKCs with the following optimum parameters: a motif width of 8–200 amino acids and a maximum of 13 motifs. The identified motifs were annotated using the program InterProScan (Quevillon et al., 2005).

The chromosomal distributions of MIKC genes were obtained based on genome annotation data. The MapInspect software was applied to draw images of their physical locations in G. hirsutum.

The expression of MIKC family genes were measured using RNA-sequencing method. The raw RNA-sequencing data of G. hirsutum TM-1 seven different tissues (root, stem, leaf, flower, ovule, seed, and fiber) was downloaded from the NCBI Gene Expression repository under the accession number PRJNA248163 (Table S4) (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA248163/). The relative data were normalized to calculate the expression levels. Hierarchical clustering was performed using Genesis 1.7.7 (Sturn et al., 2002).

Gossypium hirsutum L. (cv CCRI24) was cultivated in the field in Zhengzhou, China. Five different tissue parts of flower: sepal, petal, stamen, carpel, and ovule were sampled, respectively at full bloom stage. Arabidopsis (Columbia-0) was used as wild type; the leaves of wild type and transgenic lines grown for 25 days were harvested. All samples were frozen immediately in liquid nitrogen and kept at −80°C for total RNA extraction. Total RNA was extracted from each sample using the TRIzol reagent (TIANGEN, Beijing, China) and treated with RNase-free DNase I. Gel electrophoresis and a Nanodrop2000 nucleic acid analyzer were employed to detect the quality of RNA. The first cDNA strand was synthesized from 1 μg total RNA using the Transcriptor First Strand cDNA Synthesis Kit version DRR047A (TaKaRa, Dalian, China). The cDNA was diluted five times for the next experiments.

The gene-specific primers used for qRT-PCR were listed in Supplementary Table S2 and S3. The G. hirsutum His3 gene and Arabidopsis Actin2 gene were used as an internal control respectively. The qRT-PCR was performed using SYBR Green (Roche) on a LightCycler480 system (Roche). Each reaction was conducted in a 96-well plate with a volume of 20 μl. The PCR cycling parameters were as follows: 95°C for 5 min, 40 cycles of 95°C for 10 s, 60°C for 10 s, and 72°C for 10 s, followed by an increase from 60 to 95°C. The relative expression levels were analyzed using the LightCycler® 480 gene scanning software. Three biological replicates were measured and each biological replicate was run three times.

We amplified GhAGL17.9 using cDNA templates from the mix of CCRI24 root, stem, leaf and flower. The amplified product was cloned into vector pCambia2301 (CAMBIA) containing the CAULIFLOWER MOSAIC VIRUS (CaMV) 35S constitutive promoter, and then, the constructed vector was introduced into Agrobacterium tumefaciens GV3101 (Clough and Bent, 1998). Floral dip method was used for Agrobacterium-mediated transformation of Arabidopsis. Positive transgenic lines were selected on MS medium containing kanamycin. To grow the transgenic lines, seedlings were sown in plastic pots filled with a nutrient soil and vermiculite mix. Then, they were grown in a culturing room at 22°C under a 16-h light/8-h dark cycle for 1 month.

HMMER and BLASTP algorithm-based searches were used to identify MIKC protein HMM profiles based on the highly conserved MADS and K-box domains. To identify the maximum number of MIKC genes in G. hirsutum, HMMs of SRF-TF, and K-box domains (PF00319 and PF01486, respectively) were extracted from Pfam database to use as queries against protein sequences from the G. hirsutum genome (https://www.cottongen.org/species/Gossypium_hirsutum/nbi-AD1_genome_v1.1). A total of 145 putative MIKC proteins were identified. To verify the results, we conducted a multiple sequence alignment and removed 35 redundant sequences. Finally, 110 MIKC protein sequences were identified by confirming their conserved domains using the Pfam web server (Figure S3). From the sequences, 92 MIKC^C genes and 18 MIKC^* genes were identified. Thus, 84% of the MIKC genes were MIKC^C in G. hirsutum (Figure 1A). The identified MIKC genes were listed with their corresponding locus tag (Table 1). We named the MIKC^C genes on the basis of their assignment to the 13 previously classified Arabidopsis, O. sativa and P. tremula subfamilies (Parenicová et al., 2003; Leseberg et al., 2006; Arora et al., 2007). Subgroup AGL17 had the greatest (13%) number of GhMIKC^C genes; however, subgroups TM8 and AGL12 had the lowest (2%) number of GhMIKC^C genes (Figure 1B). The GhMIKCs' encoding amino acids were relatively conserved, and MIKC^C proteins were highly conserved, ranging from 200 to 300 amino acids in most cases. MIKC^* proteins generally possessed more than 300 amino acids. The chromosomal locations of the 108 GhMIKCs were distributed in different subgroups of the A and D genomes, while GhAGL17.12 and GhAP3.10 were located on scaffolds.

To examine the phylogenetic relationships among G. hirsutum MIKC proteins and to categorize them within the established subfamilies from other plants, we performed a multiple alignment analysis using the neighbor-joining method of 110 full-length MIKC proteins from G. hirsutum, 44 MIKC proteins from V. vinifera, 46 MIKC proteins from Arabidopsis, and 41 MIKC proteins from O. sativa (Table S1). The MIKC^C proteins were divided into 13 subfamilies (SVP, BS, AGL17, AGL15, AP3-PI, AGL12, SOC1, AG/SHP/STK, AP1/FUL, AGL6, SEP, TM8, and FLC; Figure 1C). The AGL17 subgroup was the largest, and the FLC subgroup was absent in the G. hirsutum genome. Additionally, no TM8 family members were found in Arabidopsis. TM8 constituted the smallest clade, having only four members, including two GhMIKCs, GhTM8.1, and GhTM8.2. The MIKC^* proteins were divided into two subfamilies.

A phylogenetic analysis revealed that our tree corresponded to those reported recently in V. vinifera and C. sativus (Díaz-Riquelme et al., 2009; Hu and Liu, 2012). The structures of the MIKC genes also helped to determine phylogenetic relationships (Figure 2). Most members had significant sequence identities in the same subfamily and similar exon-intron structures, indicating close evolutionary relationships. The most important differences were in the exon-intron lengths (Figure 2B). In general, most members contained eight exons in the SEP, AGL6, and AP1 gene families (except GhAGL6.1, GhAGL6.4, GhAP1.2, GhAP1.6, and GhAPI.8). The SVP (other than SVP1) and AGL12 subgroups had seven exons, whereas GhAGL15.1 and GhAGL15.4 of the AGL15 subgroup had four exons, which was consistent with GhSVP1 of the SVP subgroup. The AGL17 genes displayed relatively longer lengths compared with other subgroup genes. Additionally, GhBS4 had 11 introns and the first exon was meaningfully shorter, while in GhSOC1.5, the second of seven introns was longer than the others. The MIKC^* had much shorter gene lengths and more introns than the MIKC^C. GhMIKC^*12 had the fourth longest intron, which distinguished it from other members of the MIKC^* family.

We used MEME to analyze MIKC proteins, and 13 conserved motifs were identified (Figure 3B). Most of the closely related MIKC proteins had similar motif type distributions in the same subfamily (Figure 3A). The most striking divergence among the subgroups was in the composition of the C-terminal domains. Motif 1 contained the MADS domain in all of the MIKC families, except GhSEP8. The highly conserved sequence logs were showed in Figure S2. The differences between I regions and K-box domains were distinctly shown in the MIKC^C and MIKC^* proteins (Figure 3B). The K-box domain contained three motifs, 2, 4, and 9, in GhMIKC^C. However, motif 4, 5, 8, and 9 were present in the GhMIKC^* K-box domain, depending on the lengths. The I region in the MIKC^C subfamily contained Motifs 3 and 6, while members of the MIKC^* contained motifs 6 and 11, which resulted in a longer I region.

Among the 25 G. hirsutum chromosomes, MIKC genes were physically located on all of the 13 A chromosomes and on 12 of the 13 D chromosomes (Figure 4). Among the 110 MIKC genes, two genes, GhAP3.10 and GhAGL17.12, could not be distributed on the G. hirsutum chromosomes, but were located on unmapped scaffolds (7,246 and 4,768, respectively). The greatest numbers of genes were located on Dt-chr12 (eight genes), followed by Dt-chr2, At-chr12, At-chr13, Dt-chr11, and Dt-chr13 (seven genes on each). In contrast, two genes were located on chromosomes At-chr1, At-chr8, At-chr10, Dt-chr9, and Dt-chr10. Only one gene was mapped on At-chr9 and Dt-chr8, and no genes were located on Dt-chr1.

To explore the expression patterns of the MIKC family genes in G. hirsutum- specific developmental processes, the 110 genes' expression profiles were detected in seven different tissues (root, stem, leaf, flower, ovule, seed, and fiber) by transcriptome sequencing (Figure 5). A heat map showed that different genes shared similar expression patterns within subfamilies. For example, the SEP, AG, AP1, AP3/PI, TM8, and AGL6 subgroups were preferentially expressed in flowers. Similarly, the SVP subfamily was expressed especially in flowers. Additionally, some SEP members (GhSEP1, GhSEP3, GhSEP4, GhSEP7, and GhSEP8), and the AG and BS subgroups, were highly expressed in reproductive organs (ovules and fibers). Simultaneously, SEP, AP1 and five of the AGL6 genes (GhAGL6.1, GhAGL6.2, GhAGL6.5, GhAGL6.6, and GhAGL6.7) were also detected in roots. Interestingly, the SOC family displayed diverse expression profiles. GhSOC1.3 and GhSOC1.7 had high expression levels in roots. In addition, GhSOC1.8 was mainly expressed in flower, while GhSOC1.1, GhSOC1.4 and GhSOC1.9 were highly expressed in leaves. GhSOC1.6 and GhSOC1.10 were exclusively and highly expressed in stems. The GhAGL15 (GhAGL15.2, GhAGL15.4, and GhAGL15.5) and GhAGL17 (GhAGL17.7, GhAGL17.9, and GhAGL17.11) subfamilies were relatively highly expressed in roots, and GhAGL17.2, GhAGL17.8, and GhAGL17.12 were expressed in flowers. Four of the GhMIKC^* genes (GhMIKC^*6, GhMIKC^*7, GhMIKC^*17, and GhMIKC^*18) had high expression levels in flowers, and GhMIKC^*7 and GhMIKC^*17 were also highly expressed in seeds.

ABCDE model genes regulate the formation of five floral organs in Arabidopsis (Sánchez-Fernández et al., 2001; Dietrich et al., 2009; Kuromori, 2010). To validate the participation of MIKC genes in regulating flowering, we selected 16 of ABCDE model orthologous genes to test their expression in five parts of floral organs (sepal, petal, stamen, carpel, and ovule) by qRT-PCR in G. hirsutum (Figure 6). GhAP1.4 and GhAP1.11 (A class) showed high expression levels in sepals, petals, and carpel. Differently, GhAP1.8 was preferentially expressed in sepal. GhAP3.5, GhAP3.6, and GhAP3.8 of the AP3 subfamily, belonging to B class, were expressed in petals and stamens. GhAG4 of the C class displayed the highest expression level in stamen. GhAG7 and GhAG8 of the D class had higher expression levels in carpel and ovules. GhSEP1, GhSEP4, and GhSEP6 (E class) were expressed in four different floral organs. GhBS2 and GhBS3 (B sister class) were mainly expressed in carpel and ovules. SOC1 accelerates the flowering time, and thus, it is involved in the promotion of floral organ formation. Therefore, high expression levels of GhSOC1.2 and GhSOC1.8 were detected in sepals, stamens and carpel. These results were consistent with the ABCDE model.

AGL17 is the biggest subgroup (Figure 1B). To further investigate the role of the GhAGL17 subfamily in plant growth and development, we transformed GhAGL17.9 into Arabidopsis (Columbia-0) driven by the CAULIFLOWER MOSAIC VIRUS (CaMV) 35S promoter. We identified 12 T₃ generation transgentic lines that showed an early flowering phenotype. QRT-PCR results confirmed that GhAGL17.9 was overexpressed in transgenic lines L1 and L3 (Figure 7). Meanwhile, the numbers of rosette leaves were significantly decreased compared with WT (Table 2). To explore the molecular mechanisms that impact the flowering time in transgentic lines, qRT-PCR was used to detect the expression of flowering-related genes in transgentic lines. LFY is a flowering integration promoting factor, AGL17 can positively regulate the expression of LFY gene (Han et al., 2008), and CO is a photoperiod pathway regulator, AGL17 acts downstream of CO (Han et al., 2008). As shown in Figure 7, the expression levels of LFY gene in lines 35S-L1 and 35S-L3 were three times higher than in the wild type. CO gene expression was not significantly increased. SOC1 is a flowering promoter that regulates different signals of the flowering pathways (Lee and Lee, 2010; Ding et al., 2013). The up-regulation of SOC1 activates downstream targets, including LFY and promotes flowering in Arabidopsis (Schönrock et al., 2006; Lee et al., 2008). Approximate four-fold increases in SOC1 expression levels were observed in two transgenic lines.

In total, 110 MIKC genes were identified based on G. hirsutum genome sequences. Based on phylogenetic relationships with Arabidopsis and O. sativa orthologs (Figure S1), the G. hirsutum type II MADS family (MIKC^C) was divided into 13 subfamilies (Figure 1C). Interestingly, an FLC subfamily was not identified in the G. hirsutum genome. Similar results were found in O. sativa, C. sativus, Z. mays, and S. bicolor genomes as well (Arora et al., 2007; Zhao et al., 2010; Hu and Liu, 2012). The FLC genes are involved in controlling flowering time through the vernalization and autonomous pathways (Helliwell et al., 2006, 2011; Greb et al., 2007). Vernalization is not required for flowering in O. sativa, C. sativus, Z. mays, and S. bicolor (Arora et al., 2007; Zhao et al., 2010; Hu and Liu, 2012). Thus, vernalization might not be essential for cotton flowering as well. In addition, we found that in most subgroups, the numbers of proteins in G. hirsutum were not doubled, compared with in the diploids Arabidopsis and O. sativa. This implied that gene duplication could give rise to the amplification of MIKC subfamily genes in a variety of forms (Flagel et al., 2008; Hargreaves et al., 2014). As previously reported, multiple duplications and diversifications in the different clades of different species cause different evolutionary constraints (Lynch and Conery, 2000; Flagel and Wendel, 2009; Airoldi and Davies, 2012).

Chromosomal assignments indicated that the gene locations were equally divided among four pairs of chromosomes (At-chr6 and Dt-chr6, At-chr7 and Dt-chr7, At-chr10 and Dt-chr10, and At-chr13 and Dt-chr13) in A as well as in D genome (Figure 4). However, five D-genome chromosomes (Dt-chr2, Dt-chr4, Dt-chr9, Dt-chr11, and Dt-chr12) contained more genes compared with the corresponding A-genome chromosomes (At-chr2, At-chr4, At-chr9, At-chr11, and At-chr12). Additionally, large numbers of MIKC genes were located on the last three chromosomes (chr11, chr12, and chr13) of both genomes. This could indicate that the current phenomena were derived from differential rates of genomic evolution and inter-genomic hereditary information transfer (Paterson et al., 2000; Wendel and Cronn, 2003).

Global expression patterns analyses in seven different tissues showed that the API, AP3, AG, SEP, and BS subfamilies were almost all expressed in the flower development stage (Figure 6). Floral organ identities and flower meristem are regulated by five kinds of genetic functional genes (A-B-C-D-E) during flower development, from sepals to ovules (Díaz-Riquelme et al., 2009; Na et al., 2014). A qRT-PCR analysis showed the expression patterns of the orthologous genes of the ABCDE model in flower organogenesis (Figure 6), which were consistent with previous findings in Arabidopsis (Ó'Maoiléidigh et al., 2014; Xie et al., 2015). Further, the API subgroup of A class genes were not only expressed in sepals and petals, but also exhibited carpel expression profiles. Before and after pollination, the API-like gene may aid in the carpel development in Orchidaceae, which triggered ovary development (Mondragon-Palomino and Theissen, 2011; Acri-Nunes-Miranda and Mondragón-Palomino, 2014). Thus, AP1 subgroup genes may have similar expression patterns in Orchidaceae and allotetraploid cotton. A few GhMIKC^* genes were highly expressed in flowers and seeds, which was in accordance with previous results in Arabidopsis (Verelst et al., 2007) and O. sativa (Liu et al., 2013). These results indicated that the expression profiles of MIKC^* genes were involved in functional redundancy and conservation in the process of G. hirsutum evolution.

In Arabidopsis, AGL17 acts as a novel flowering promoter, which is involved in the photoperiod pathway. Under long-day conditions, the overexpression of AtAGL17 causes early flowering (Han et al., 2008). As the largest subgroup of the GhMIKC^C family, one member of the AGL17s, GhAGL17.9, was overexpressed in Arabidopsis to explore its biological functions. The transgenic lines displayed earlier flowering than wild type (Figure 7). The expression levels of the related positive marker genes, especially LFY and SOC1, which are involved in regulating the flowering process, were higher in p35S::GhAGL17.9 lines than in wild type. LFY overexpression can prematurely cause plant development and accelerate blossoming processes (Nilsson et al., 1998; Dornelas and Amaral, 2004). AGL17 targets LFY to promote flowering (Han et al., 2008). SOC1 encodes a MIKC protein, a floral pathway integrator, which is regulated by a variety of flower signaling pathways (Lee et al., 2000; Wang et al., 2009; Ding et al., 2013). However, the relationship between AGL17 and SOC1 in flowering is not clear, which remains to be functionally explored further in the future.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.