The cotton varieties CCRI50, Jiumian2, TM-1, and Yumian8 used in this study were preserved in Hebei Agricultural University and grown in a greenhouse (16 h light/8 h dark, 28°C day/25°C night) for experiments. The TM-1 variety was grown for tissue-specific analysis. The roots of the seedlings at the cotyledon stage and the stem, leaf, and SAM of seedlings at two true leaf stages (TLSs) were sampled, respectively. The calycle, sepal, petal, stamen, pistil, and ovule were collected 0 day postanthesis (DPA). Fibers were separated from ovules at 5 DPA. Two early-maturity cotton varieties CCRI50 [whole growth period (WGP) in Yellow River Basin Region (YRBR) is 110 days] and Jiumian2 (WGP in YRBR is 114 days) and two late-maturity cotton varieties TM-1 (WGP in YRBR is 135 days) and Yumian8 (WGP in YRBR is 136 days) were used for temporal expression analysis. A total of 10 shoot apexes were collected, respectively, at two, three, four, and five TLSs. The Arabidopsis and tobacco plants were grown in a plant growth chamber under 16 h light/8 h dark, 22°C. The whole seedlings of Arabidopsis homozygous plants were collected for expression analysis. All samples were immediately frozen in liquid nitrogen for further analysis.

The four Arabidopsis SEP1/2/3/4 protein sequences were downloaded from the TAIR (https://www.arabidopsis.org/) and then blast against the published genomes of G. hirsutum on CottonFGD (https://cottonfgd.org/). The protein sequences and functional annotations were filtered for the protein family database (Pfam) identifiers of the MADS and K domains (PF00319 and PF01486), respectively. The candidate GhSEPs were confirmed with HMMER 3.0 and the Batch CD-Search service. The homologous genes of SEP in Gossypium (G. raimondii, G. arboretum, and G. barbadense) were determined in the same way.

The phylogenetic trees were constructed by MEGA 7.0 using the neighbor-joining (NJ) method with 1,000 bootstrap replications and default parameters and then displayed with the online iTOL tool (https://itol.embl.de/) (Kumar et al., 2016; Letunic and Bork, 2019). Segmental and tandem duplications were detected by MCScanX with default parameters. Homologous genes between the At and Dt subgenomes were determined using the bidirectional best hit method in BLAST. The duplication events were fetched and then displayed with TBtools (Chen et al., 2020).

Plant total RNA was extracted using RNA Easy Fast Plant Kit (TIANGEN Beijing, China). DNase treatment was performed with RNase-free DNase (TIANGEN, Beijing, China) before purification. The cDNA was synthesized using cDNA Synthesis SuperMix (TRANS, Beijing, China). Gene transcription was detected by quantitative real-time PCR (qPCR) using AugeGreen™ Master Mix for HRM (US EVERBRIGHT, Suzhou, China) on the ABI 7500 PCR Detection System (USA). Gene-specific primers for qPCR were verified for its specificity according to the single peak in the melting curve and are listed in Table S1 . The relative expression level was calculated using the 2^-△Ct formula. The expression was normalized to AtTUB2 (AT5G62690) in Arabidopsis and GhHis3 (GhM_D03G0424.1) in cotton and shown as relative values to the maximal gene expression level set at 100%. Three biological repeats were applied on each sample, and three technical repeats were performed on each reaction. The standard deviation (SD) of three biological repeats was calculated. R.E.L. stands for the relative expression level.

The coding regions of GhSEPs were cloned (primers listed in Table S1 ) and purified by the E.Z.N.A.^® Gel Extraction Kit (Omega Bio-tek, Norcross, Georgia, USA). The resulting fragments were ligated into vector pGreen0229 with a 35S promoter. The constructs were introduced into Agrobacterium tumefaciens strain GV3101 and transformed to Arabidopsis wild-type (WT) plants (Columbia) using the floral dip method. For each GhSEP gene, at least 25 T1 individual lines were obtained and confirmed by genotyping for the exogenous fragment. Then, five homozygous T3 lines were randomly selected and saved for a detailed observation of phenotypes. Expression analysis was performed on the seedlings of homozygous plants.

The 2kb promoter sequence of each GhSEP was cloned and, finally, the genome regions of GhSEP1A, GhSEP2D, GhSEP3.1A, GhSEP4.1D, and GhSEP4.2D were obtained and inserted into the pGreen-GUS vector to generate pro:GhSEP-GUS constructs. Then, the construct was transformed into Arabidopsis using the methods described above. The homozygous lines were used for the histochemical assays of GUS activity. Different tissues were collected for GUS staining, including 13-day-old seedlings and cauline leaves, inflorescence, and floral organs from 6-week-old plants. Samples were immersed in the staining solution (Coolaber, Beijing, China) and incubated at 37°C overnight. Then, the samples were decolored in 70% (v/v) ethanol twice until the negative-control material (WT) turned white. Stained and cleared specimens were visualized and photographed using a stereoscope (AxioCam ICc 5; Carl Zeiss, Jena, Germany).

The SAM of VIGS plants were embedded in Optical Cutting Temperature (OCT) compound (Leica Wetzlar, Germany) and fast-frozen in liquid nitrogen. Then, the samples were cross-sectioned by the freezing microtome (7500; Leica, Wetzlar, Germany) and the serial sections were expanded on a plus on a slide (CITOGLAS, Beijing, China) followed by observation under a microscope (DM2500; Leica Wetzlar, Germany).

The SAM of the early- and late-maturity cotton mentioned above was sampled and dehydrated in ethanal serial solutions (from 20% (v/v) to 100%). The tissues were visualized by adding eosin in 70% ethanal. Then, the tissues were incubated in a serial solution of ethanal/histonclear and histonclear/paraffin chips. The samples were finally embedded in the paraffin for section using a microtome (Leica, HistoCore AUTOCUT, Germany). The paraffin ribbon was expanded on plus on slides, and the paraffin was removed by a serial solution of ethanal/histonclear. The morphological characteristics of the SAM at each period were recorded and photographed under a microscope (DM2500; Leica, Germany).

The coding regions of candidate genes were cloned into pGADT7 (AD) and pGBKT7 (BD) vectors (Clontech, San Francisco, California, USA). The two-hybrid assay was performed according to the Matchmaker^® Gold Yeast Two-Hybrid System (Clontech San Francisco, California, USA).

The open reading frames (ORFs) of full-length genes were inserted into separate pSAT1A-nEYFP-N1 and pSAT1-cEYFP-C1-B vectors and then transformed into Agrobacterium (GV3103). These Agrobacteria were coinfiltrated into Nicotiana benthamiana leaves. The florescence signals were detected under confocal microscopy (FV10i; Olympus, Tokyo, Japan).

The coding region of GhSEP2 and GhSEP4.2 were cloned into a 35S:-6HA vector (pGreen). The mesophyll protoplasts transient expression system (Li et al., 2021b) was applied for the ChIP assay. The fifth leaves and shoot apexes from TM-1 at five TLSs were collected and sliced into fine filaments followed by gentle digestion for 9 h (1.5% cellulose, 0.4% macerozyme, 0.5 mol·L^-1 mannitol, 20 mmol·L^-1 KCl, 20 mmol·L^-1 MES, 10 mmol·L^-1 CaCl₂, and 1.0 g·L^-1 BSA) to extract protoplasts. The plasmids 35S:GhSEP2-6HA and 35S:GhSEP4.2-6HA were transiently expressed in the protoplast by 40% PEG under isotonic pressure maintained by 0.5 mol·L^-1 mannitol. Then, the protoplasts were cultured with a WI buffer (4 mmol·L^-1 MES, 0.5 mol·L^-1 mannitol, 20 mmol·L^-1 KCl) in the dark for 16 h. More than 10^-15 million cells were collected and fixed in 1% formaldehyde on ice for 20 min. The nuclei were isolated and sonicated to produce DNA fragments approximately 500 base pairs. Nucleoprotein was then immunoprecipitated by an anti-HA antibody conjugated with agarose (Sigma, St. Louis, Missouri, USA). The DNA fragments were purified by the E.Z.N.A.^® Gel Extraction Kit (Omega Bio-tek, USA). Western blot was applied to detect the fusion protein using an HA antibody (Invitrogen, USA). The relative enrichment of each fragment was determined by qRT-PCR. Primers used for the ChIP assay are listed in Table S1 . The enrichment fold of each fragment was calculated first by normalizing the amount of a target DNA fragment against GhHis3 as an internal control and then by normalizing the value for transgenic protoplasts against non-transformed protoplasts. The SD of three biological repeats was calculated.

The MADS box family is highly conserved in the MADS-box domain. To distinguish SEP-like genes from other MADS family members in Gossypium, the four Arabidopsis SEP protein sequences were blasted against the published genomes. The results were filtered referring to the sequence similarity and the identified members of the cotton MADS-box family (Ren et al., 2017; Nardeli et al., 2018). Then, the selected sequences were confirmed by the phylogenetic tree constructed with SEP homologs ( Table S2 ) and named according to the members in the same subgroup ( Figure 1A and Table S3 ). Finally, 12 SEP-like genes were identified and distributed equally in the At and Dt subgenomes, which were named according to the sequence similarity ( Figure S2A and Table S4 ). Consistent with the protein sequence similarity of each homologous gene pair in the Gh genome, their gene structures were conserved ( Figures S2B, C ). The GhSEPs contained MADS-box and K-box domains, and their C-terminals possessed a characterized SEP I motif and SEP II motif ( Figures S2A, D ) (Zahn et al., 2005; Pan et al., 2014). We also found conserved motifs between the cotton and Arabidopsis SEP1/2 and SEP3 proteins, respectively, designated as the SEP1 motif and SEP3 motif, respectively ( Figure S2A ).

The phylogenetic tree was clustered into four groups ( Figure 1A ). The SEP3 group was closer to the SEP1 group, and the two groups were then rooted with the SEP2 group. The SEP4 group was located in a single evolutionary branch. Only the SEP3 group contained genes from all the chosen species, and other SEP-like genes in different species were clustered irregularly to other groups, suggesting that SEP3-like genes were mostly conserved in the angiosperm. Interestingly, SEP4-like genes in Gossypium were closely clustered with the rest of monocot SEP-like genes except for SEP3, indicating a different evolution of SEP4-like genes. These results supported that SEP-like genes originated from a single MADS gene and then extended during the whole genome duplication events followed by the evolution of diverse functions (Smaczniak et al., 2012; Chen et al., 2017; Zhao et al., 2017). To further elucidate the evolution of SEPs in the Gossypium species, their synteny was analyzed. The results showed the segmental duplications of SEP3 genes in the diploid Gr and Ga genomes ( Figure 1B ), and the diploid SEP3 sequences were highly similar to the SEP3 genes in the tetraploid Gh genome ( Figure S2A ). Therefore, we speculate that SEP genes in the tetraploid Gh genome originated from their diploid progenitors.

The expression pattern of GhSEPs was studied first using the transcriptomic data obtained from CottonFGD (Hu et al., 2019). The transcripts of GhSEP homologies from At and Dt subgenomes displayed a consistent expression pattern in the heat map ( Figure S3A ). GhSEP1/2/3.1 were highly transcribed during fiber development, whereas GhSEP4s were undetectable in the ovule or fiber. GhSEP1 transcribed relatively higher during seed germination, and GhSEP4.1 could be induced by cold, drought, and salt stress. However, the transcriptomic data lacked information about the GhSEP3.2 gene, so we further compared the transcription of GhSEP genes in our unpublished transcriptomic data of the leaf and SAM during cotton floral transition ( Figure S3B ). The similar expression of GhSEP homologies from At and Dt subgenomes was also observed. Additionally, GhSEP genes displayed higher expression after floral transition (at five TLSs), especially in the SAM. The transcripts of GhSEP3.2A and D were not expressed in the leaf and SAM. The expression of GhSEPs were further examined in vegetative and reproductive tissues. Due to the extraordinarily high sequence similarity of homologies in At and Dt subgenomes, their transcripts were detected together and shown as GhSEP1, GhSEP2, GhSEP3.1, GhSEP4.1, and GhSEP4.2. The expression pattern of GhSEPs was different in flowers from the outside to the inner whorl ( Figure S3C ), which is similar as their homologies in Arabidopsis. GhSEP4.1 and GhSEP4.2 were expressed in the calycle. Except for GhSEP3.1, the other GhSEPs were expressed in the sepal. All the GhSEPs participated in the development of the petal and stamen. Only GhSEP3.1 showed considerable transcripts in the pistil. Noticeably, GhSEP2 expression raised up in the fiber, suggesting a novel function of cotton SEP genes in fiber development. In the SAM, GhSEPs consistently displayed higher expression, especially GhSEP4.1 and GhSEP4.2 ( Figure 2A ). In addition, GhSEP2 and GhSEP3.1 transcription in the SAM was relatively lower. Only GhSEP4.2 showed a much higher expression in the vegetative tissues of leaves. Since SEP-like genes have been shown to control flowering time, we monitored the expression of GhSEPs in developing seedlings. Before expression analysis, the cytological morphology of the SAM was observed in early- and late-maturity cottons to illustrate the occurrence of floral transition. The eminence of the SAM was obvious in the early-maturity cotton CCRI50 and Jiumian2 when the fourth true leaf flattened, while the eminence was just visible in the SAMs of TM-1 and Yumian8 at five TLSs, suggesting a proximate 5-day delay of the floral transition in late-maturity cotton varieties ( Figure 2B ). Consistent with the occurrence of floral transition, the expression of GhSEPs was upregulated at the four TLSs in the early-maturity cotton varieties, while their expression only increased slightly in late-maturity cotton varieties. The GhSEP expression level was gradually increased along with the floral transition ( Figure 2C ), suggesting that they promote flowering time. In line with this, a higher expression level of GhSEPs was detected in the early-maturity cotton varieties than the late ones during floral transition ( Figure S3D ). The transcription of GhSEPs in tissue and developing seedlings suggested redundant roles in floral organ development and flowering time control.

Furthermore, we cloned the 2 kb upstream promoters of GhSEP1A, GhSEP2D, GhSEP3.1A, GhSEP4.1D, and GhSEP4.2D whose coding regions were obtained for functional study. GUS reporter lines were constructed driven by these promoters to mimic a detailed expression of GhSEP genes in Arabidopsis ( Figure 2D ). GUS signals were only observed in the SAM in pro:GhSEP1/2/3.1-GUS Arabidopsis lines. GhSEP4.2 was expressed in the vascular bundle. GUS signals in pro:GhSEP4.1-GUS were observed in the whole leaf. After bolting, GUS signals were all detected in cauline leaves, flowers, and siliques. Interestingly, GhSEP4.1 and GhSEP4.2 displayed a consistent expression pattern in the cauline leaf and sepal as in the true leaf. The GUS signal of pro:GhSEP3.1-GUS was strong in the reproductive organ. GhSEP2 showed unique transcription in ovules consisting of the qPCR results. It worth noticing that the GhSEP4.1 transcriptional level in Arabidopsis leaves shown by GUS staining is not consistent with the low expression level in TM-1 cotton leaves detected by qPCR ( Figures 2A, D ). The dynamic transcription of GhSEP genes in a plant life cycle indicated at least partially independent functions.

Then, the coding region of each GhSEP gene was cloned and overexpressed with GFP in the tobacco epidermal cells. Consistent with their roles as a transcription factor, the fluorescent signals of the GhSEPs-GFP were detectable in the nuclei ( Figure 2E ).

In summary, the overlapping and distinguished expression patterns of GhSEPs suggested functional redundancy and diversity. GhSEP4s might function dominantly during floral transition, and each GhSEP gene was likely to participate differently in the formation of each whorl of flowers. The prevalent expression of GhSEP2 in ovule indicated a novel function.

To elucidate the function of GhSEPs, their coding sequences were amplified and the resulting fragments shared the same sequences with GhSEP1A, GhSEP2D, GhSEP3.1A, GhSEP4.1D, and GhSEP4.2D after blasting against the G. hirsutum genome. Considering the high sequence similarity of GhSEP homologies in At and Dt subgenomes, the above sequences of GhSEP genes were overexpressed in Arabidopsis for functional study. There were more than 25 individual T1 lines of each GhSEP gene that flowered uniformly earlier than WT plants ( Figures 3A, B ). Some individuals even flowered with two rosette leaves ( Figure S4A ). Five lines of each transgene were randomly selected for transcriptional and phenotypic analysis. The results showed that the flowering time of transgenic plants reduced up to half of the WT plants ( Figure 3B ) and the flowering phenotype was closely associated with the expression levels in 35S:GhSEPs ( Figure 3C ), suggesting that GhSEPs promoted flowering time in a dosage-dependent manner.

Next, due to the similarity of the GhSEPs’ coding sequences, specific fragments from the 3’ UTR of each GhSEP gene were amplified and constructed into the pTRV vector for the virus-induced gene silencing (VIGS) assay applied on an early-maturity cotton CCRI50 to silence the native GhSEP genes. The silencing effects were examined in the SAM ( Figure 4A ). VIGS caused a decreased expression of the relative gene in the silencing plants of TRV : GhSEP3.1, TRV : GhSEP4.1, and TRV : GhSEP4.2. However, the silencing of GhSEP1 and GhSEP2 was synchronal due to high sequence similarity. The silencing efficiency of each gene was above 70%. Compared with the control group, the silencing of GhSEPs delayed cotton squaring and flowering time ( Figures 4B and S5A ). The cross-section of the shoot apex at four TLSs demonstrated floral and inflorescence meristems of the control plants indicating the stage of reproductive growth, whereas the SAMs of TRV : GhSEPs were still flattened and only leaf primordia were observed, suggesting that the floral transition has not yet occurred ( Figure 4C ). A statistical analysis of the traits reflecting cotton maturity suggested that the silencing of GhSEPs resulted in a significant delay of maturity and the effects of GhSEP3.1/4.1/4.2 were comparable ( Figure 4D ). The late-maturity phenotypes of TRV : GhSEP1 and TRV : GhSEP2 were more severe, which might be resulted from the synchronal silencing of GhSEP1 and GhSEP2. These results described the participation of cotton SEP-like genes in floral transition in a dosage-dependent manner.

SEP genes are well known as the E function in floral development. Therefore, the organ morphological phenotypes were observed in the overexpression of Arabidopsis and the silencing of cotton plants. Some homozygous plants of 35S:GhSEPs produced curly leaves ( Figure S4B ). The expression of GhSEPs in those lines were much higher compared with the lines that possessed normal rosette leaves and flowered later ( Figure 3 ). However, the homozygous 35S:GhSEPs consistently showed no significant perturbation of a flower ( Figure S4C ). Similarly, the silencing of any GhSEPs in the cotton plants had no effects on floral development ( Figure S5B ). However, some heterozygotes of the Arabidopsis 35S:GhSEPs T1 plants produced abnormal flowers with absent sepals from a leaf-like structure that was similar to bracts ( Figure S4A ). Moreover, terminal flowers were observed in some individuals. These plants flowered extremely early and failed to produce any seeds. These severe phenotypes could be caused by the higher expression level of any GhSEPs. These results suggested that GhSEPs acted in a dosage-dependent manner.

To elucidate the regulatory mechanism of GhSEPs, the expression of flowering-time genes was first detected in 35S:GhSEPs Arabidopsis plants. The transcription of LFY and AP1 was dramatically upregulated in each of the GhSEP overexpression lines. In addition, the expression of the other flowering time regulators FT, SOC1, CONSTANS (CO), and AGL24 was increased significantly ( Figure S6 ). Although GhSEPs shed similar effects on LFY and AP1 transcription, they upregulated the other detected genes differently in the overexpression lines. GhSEP1 and GhSEP4.1 were likely to perform similarly in flowering time control that SOC1 expression change was weaker than FT, CO, and AGL24. GhSEP2 had an even effect on the examined flowering time genes. Only FT transcriptional increase was violent in the 35S:GhSEP3.1 lines. Additionally, GhSEP4.2 regulated FT and SOC1 at the similar level that was stronger than the regulation of CO and AGL24. The discrepancy of the effects of GhSEPs on flowering time regulators indicated the distinct roles of each GhSEP in flowering time control. Then, the expression trends of LFY, AP1, FT, and SOC1 were monitored in Arabidopsis developing seedlings. AP1 and LFY expression increased rapidly 5 days after germination of the 35S:GhSEPs plants, while they remained at a relatively low level in the WT plants, suggesting that the floral transition occurred earlier in the 35S:GhSEPs plants. However, the expression trends of FT and SOC1 were similar in the GhSEP overexpression and WT plants, except for the FT expression trend in the 35S:GhSEP4.1 ( Figure 5A ). These results suggested that the upregulation of SOC1 and FT in the 35S:GhSEP plants might be indirect.

Furthermore, the expression of flowering time genes was analyzed in the TRV : GhSEP silencing plants ( Figure 5B ). Consistent with the expression change in Arabidopsis, the silencing of GhSEPs all caused significant expression decreases of GhAP1 and GhLFY, among which the same expression variation was found in TRV : GhSEP1 and TRV : SEP2 that was greater than that in TRV : GhSEP3.1 and TRV : GhSEP4s, which was probably due to the cosilencing effects on GhSEP1/2. However, GhFT and GhSOC1 expression was undisturbed by the effective silence of GhSEPs ( Figure S5C ), which supported our speculation on the transcriptional increases of SOC1 and FT in the 35S:GhSEPs of Arabidopsis.

Therefore, LFY and AP1 might be the downstream target genes of GhSEPs in the regulation of flowering time. To prove this, ChIP analysis was performed to examine the direct binding of GhSEPs to the GhLFY and GhAP1 genomic DNA. In the upland cotton genome, there were four copies of GhAP1 and two copies of GhLFY whose transcriptional regulation regions were different from each other. Thus, the CArG-box elements were predicted on their whole genome regions. The results showed that GhAP1-A04 and -D04 contained denser CArG-boxes, especially in the first introns ( Figure 5C ). Considering the higher expression of GhSEP4.2 in cotton leaves ( Figure 2A ) and the specific expression of GhSEP2 in the seed and fiber ( Figure 2D and S3C ), the DNA fragments bound with GhSEP2-HA and GhSEP4.2-HA were extracted from cotton protoplasts overexpressing the fusion proteins ( Figure S7 ), respectively. The results suggested that GhSEP2 was strongly bound to the genomic region of GhAP1-A04 around the start codon and the other two regions distributed in balance of the up- and downstream of the start codon ( Figure 5D ). The binding distribution of GhSEP2 to the GhAP1-D04 was similar to that of GhAP1-A04, but the fold enrichment was quite lower than GhAP1-A04. Fold enrichments were also detected at a low level around the start codon of GhAP1-A13 and the middle region of the first intron of GhAP1-D13. GhSEP2 was also strongly bound to GhLFY-A07 at the end of second intron and weakly bound to the genomic region near the stat codon of GhLFY-D07. Different from GhSEP2, GhSEP4.2 did not bind to GhAP1-D13 and displayed no binding superiority to the other GhAP1 and GhLFYs. Two binding regions were detected located at the first and last introns of GhAP1-A04 and 1.8 kb up- and downstream of the start codon of GhLFY-A07. GhSEP4.2 was bound to the transcription start region of GhAP1-D04 and GhAP1-A13 and the 1 kb upstream of the start codon in the GhLFY-D07 genome. Additionally, it seemed that GhSEP2 binding affinity was stronger than GhSEP4.2, which is supported by in vitro EMSA results (Jetha et al., 2014), and GhSEPs were likely to loop DNA that was explained by the separated CArG-boxes in two binding sites (Melzer et al., 2009).

DNA looping is formed with the binding of MADS tetramers (Theissen and Saedler, 2001). Hence, the protein interaction of GhSEPs were analyzed. First, we failed to study the dimerization between GhSEPs in two directions in yeast cells due to the strong self-activation of GhSEP2. Subsequently, the homodimers and heterodimers of GhSEPs were visualized by the bimolecular fluorescence complementation (BiFC) assay ( Figure 6A ). All the GhSEP proteins were capable to form heterodimers between each other, and homodimerizations were observed, except for GhSEP2. Then, the interaction between GhSEPs and flowering time regulators were tested, including GhAP1 (GhM_A13G0958), GhSOC1.4 (GhM_A11G0098), and GhSVP.5 (GhM_A12G1181). Strong interactions between GhSEPs and GhAP1, GhSOC1.4, and GhSVP.5 were detected both in the yeast two-hybrid assay and BiFC analysis ( Figures 6B, C ). These results suggested that GhSEPs were able to form tetramers with flowering time regulators and themselves.

Plant MADS-box transcription factors play essential roles in almost every developmental process. Many MADS-box genes have conserved functions across the angiosperms, but some have evolved novel functions in specific species. The morphological diversity of flower organs is closely related to functional divergence within the MADS-box genes. An evolutionary analysis of the MADS genes in plants points out a close link of SEP-like genes’ origination with angiosperm evolution, especially in the formation of bisexual flowers (Zahn et al., 2005; Theissen and Melzer, 2007; Ruelens et al., 2017). The overexpression of cotton SEP-like genes affected inflorescence and floral development ( Figure S4A ), which is conserved in the angiosperm. Aside from floral development, the flowering time variation was observed in the overexpressing and silencing plants of a single GhSEP ( Figures 3 , 4 ), revealing the roles of GhSEPs in promoting flowering. However, the mutation or silencing of single SEP-like genes in Arabidopsis, petunia, orchid, rice, and tomato does not affect flowering time although an overexpression of them causes the early flowering of Arabidopsis (Ditta et al., 2004; Matsubara et al., 2008; Pan et al., 2014; Wu et al., 2018; Zhang et al., 2018). These phenotypes suggested that the function of GhSEPs in floral development is largely redundant, but their function in flowering time control is less redundant.

SEP-like genes share a most recent common ancestor with AGL6 and AP1 genes 296 MYA (Shen et al., 2019) and duplicated in the whole genome duplication event in flowering plants (Zhao et al., 2017). It could be observed in the phylogenetic tree of MADS-box genes that the lineages of SEP-, SQUA-, FLC-, and AGL6-like genes are nested within a strongly supported superfamily by evolutionary conserved tandem duplications between SEP1 and SQUA, as well as SEP3 and FLC (Ruelens et al., 2013; Zhao et al., 2017). Syntenic relationships are also supported by the conserved tandems of SEP-SOC1 and SEP-SVP in the angiosperm (Zhao et al., 2017). The evolutionary evidence all points to a close relation of SEP-like genes with MADS flowering time regulators. The members of the SOC1-, SVP-, FLC-, and SQUA-like genes are generally the regulators of floral transition that determines the flowering time and reproductive success (Bowman et al., 1993; He, 2012). Additionally, the functional characterization of gymnosperm AGL6- and SOC1-like genes CjMADS14 and CjMADS15 show their conserved function in flowering time control (Katahata et al., 2014). Consistently, the clusters of SVP, AG, AP1, and SOC1 genes diverge sequentially in this order to form the SEP genes in the lineage of cotton MADS genes (Nardeli et al., 2018), suggesting that cotton MADS transcription factors evolute conservedly. Therefore, we speculate that the function in the flowering time control of GhSEP genes was preserved from the ancestor during evolution.

It has always been considered that SEP-like genes in angiosperms are clustered into two groups: the SEP1/2/4 and SEP3 (Shen et al., 2019). However, cotton SEP4 genes diverged independently with AtSEP4 and monocot SEP-like genes other than SEP3 orthologs ( Figure 1A ). In the SEP4 clade, OsMADS34, OsMADS1, and OsMADS5 contribute to the origin of distinct grass inflorescences and spikelets (Gao et al., 2010; Meng et al., 2017), whereas the other cotton SEP genes were clustered with SEP1/2/3 orthologs that function as flowering time promoters and floral organ determinators (Chung et al., 1994; Ferrario, 2003; Zhang et al., 2017; Lin et al., 2020). The sequence diversity and different conserved motifs in the C terminal of GhSEPs also supported the functional divergency ( Figure S2 ) (Vandenbussche et al., 2003a). Consistently, the expression of each GhSEP gene differed spatially and temporally ( Figure 2 ), suggesting that they have diverse functions. Hence, we propose that the neofunctionalization and redundancy that occurred following gene duplication jointly contribute to the functions of GhSEPs (Zahn et al., 2006).

Dosage dependency is common in the regulation of quantitative traits. This principle coincided with the early flowering phenotype of Arabidopsis overexpression lines varied along with the expression level of GhSEPs ( Figure 3C ). The sterile 35S:GhSEPs plants flowered with only two or fewer leaves ( Figure S4A ), which was reasonably caused by the excessive transcription levels of GhSEPs. Under these circumstances, floral defects were visible. In Arabidopsis, the inactivation of one SEP1 allele in a sep2 sep3 double-mutant background converts the normal flowers to severe abnormalities in ovule development (Favaro et al., 2003). A further reduction of SEP1 activity in the sep1 sep2 sep3 triple mutant results in flowers consisting of sepals only (Pelaz et al., 2000). The flower formation is fully disturbed in sep1 sep2 sep3 sep4 quadruple mutants that only leaf-like organs exist (Ditta et al., 2004). The defects of rice spikelet progress in osmads1-z osmads5-3 osmads34-1 triple mutants compared with double mutants osmads1-z osmads5-3 and osmads1-z osmads34-1 (Wu et al., 2018). The phenotypic evidence all reveals that the function of SEP proteins in floral development depends on their concentration.

The in vitro binding assay of SEPs and CArG fragments proved that protein concentration affects the binding affinity of SEPs and their cooperative complexes (Jetha et al., 2014; Rumpler et al., 2018), which lay a molecular foundation of the concentration-dependent regulation of SEP proteins to their target genes. Furthermore, the regulation of zygotic genes by maternal-effect genes during the Drosophila development sets a paradigmatic example for concentration-dependent target gene regulation (Lander, 2007). Whether SEP protein concentration is responsible for different target gene regulation lacks evidence. The Arabidopsis GhSEP overexpression lines could be summarized into three types according to the phenotypes ( Figure 3 ). Type I only flowered earlier, and GhSEP expression levels were the lowest in the plants compared to the other transgenic lines. In Type II plants, the transcripts of GhSEPs were more abundant than Type I plants and curled leaves were additionally observed. Type III plants further demonstrated floral organ defects due to the excessive expression. The phenotypic variation is controlled by different genes downstream of GhSEPs. Therefore, a relatively minor accumulation of GhSEP proteins suffice to promote flowering, but the function in floral development requires more abundant GhSEPs.

Although the FQM describes the mechanism of SEP in floral development, the regulatory mechanism in flowering time control remains unclear. The overexpression of SEP3 causes similar phenotypes as 35S:AP1 to promote the flowering of Arabidopsis, which could be masked by an elevated expression of LFY. As a floral meristem identity gene, LFY is a downstream of SEP3 (Castillejo et al., 2005), and the activation of LFY will, in turn, accelerate floral transition (Moon et al., 2005). OsMADS7 (also named OsMADS45) is a SEP3 homolog. The overexpression of OsMADS7 can overcome photoperiod to upregulate Hd3a and RFT1 simultaneously resulting in early flowering (Wang et al., 2013). The expression changes of FT, SOC1, CO, LFY, and AP1 in 35S:GhSEPs were consistent with the overexpression lines of other SEP orthologs ( Figure S6 ) (Tzeng et al., 2003; Pu et al., 2020). However, the expression trend of FT and SOC1 in the 35S:GhSEP plants remained the same as in the WT, whereas the transcriptional elevation of LFY and AP1 advanced suggesting earlier floral transition ( Figure 5A ). Meanwhile, the late flowering phenotype induced by the decrease of GhSEP expression was only related to the expression changes of GhAP1 and GhLFY ( Figures 5B and S5C ). Thus, AP1 and LFY were the downstream genes of GhSEPs in flowering time control. A genome-wide identification of AtSEP3 and OsMADS1 (SEP4 homolog) targets both demonstrate the direct binding of them to the gibberellic acid (GA) biosynthesis genes to regulate the GA level that is involved in regulating flowering time (Kaufmann et al., 2009; Khanday et al., 2016). Thus, the expression changes of other flowering time genes in 35S:GhSEP plants might be due to the regulation of GA biosynthesis.

Previous studies revealed the transcription of GhAP1 in the leaf and SAM (Cheng et al., 2021), whereas GhLFY expression was mainly detected in the SAM (Li et al., 2013). Only GhSEP4.2 displayed commensurate expression in leaf with GhAP1 ( Figure 7A ). Thus, we hypothesize that GhAP1 is one of the targets of GhSEP4.2 in the leaf, where GhSEP4.2 interacts with GhSOC1, GhSVP, and itself to form tetramers ( Figure 7B ). In the SAM, GhSEPs were highly expressed at different time forming dynamic tetramers ( Figure 7B ). GhSEP4.1 and GhSEP4.2 were mainly expressed in the inflorescence meristem and maintained reproductive growth. They were also expressed earlier in the floral meristem to determine the formation of the calycle, sepal, and petals. GhSEP1 was then expressed to participate in the growth of the sepal, petal, and stamen. GhSEP2 and GhSEP3.1 were expressed latest to regulate petal and stamen formation. In addition, the pistil determination only involved GhSEP3.1. GhSEP2 expression was upregulated in the ovule when the fiber initiates. The fine-tuning of the transcription of GhSEP genes was essential for their function in reproductive growth, and the resulting protein tetramers, in turn, regulated different target genes in the SAM.

MADS proteins form tetramers to loop DNA, and this process requires only a limited amount of interaction energy between the DNA-bound dimer (Melzer et al., 2009; Jetha et al., 2014). According to this principle, GhSEP4.2 binds to GhAP1-A04 easier than other GhAP1 regions as two separated binding sites were identified with an equal enrichment ( Figures 5D and Figure 7 ). GhSEP4.1/4.2 accumulated dramatically due to their numerous transcripts, which is the same as Arabidopsis SEP4 that transcribes earliest in the SAM (Ditta et al., 2004). The transcriptional superiority usually contributes to a dominant role. However, GhSEP4.2 bindings in GhAP1 genome region were comparatively weaker than GhSEP2 binding to GhAP1-A04 ( Figure 5D ). Moreover, GhSEP2 was unable to self-dimerize ( Figure 6A ). GhSEP2 binding to GhAP1-A04 must be coordinated. GhSEP1, whose expression was higher in the SAM, was considered rather than GhSEP4.2 because the binding sites of GhSEP2 did not overlap with that of GhSEP4.2. Inferring from these, GhSEP1 and GhSEP2 regulate GhAP1 dominantly on A04 loci ( Figure 7B ). Furthermore, the binding sites of GhSEP4.2 coincide with that of GhSEP2 on GhAP1-D04 and -A13. They form a heterodimer to regulate these two loci. Separated binding sites suggested a loop of GhAP1-D04. However the weak binding of GhSEP2 on GhAP1-A13 should be strengthened via other GhSEPs. GhSEP4.1, GhSEP4.2, or GhSEP1, rather than GhSEP2, might cooperate with the regulation on GhAP-A13 due to their stronger interaction ability. Meanwhile, the binding on GhAP1-D13 was only detected by GhSEP2 with a relatively low enrichment indicating a weak regulation of GhSEPs on this locus. Next, the binding superiority for GhLFY was also in the At subgenome. Heterodimerization is responsible for the strong binding. The similar enrichment and nearby binding sites of GhSEP2 and GhSEP4.2 suggested a cobinding on GhLFY-D07 ( Figure 5D ). This regulatory model also needs the participation of flowering time regulators, such as GhSOC1 and GhSVP ( Figures 6B,C ). Members in a tetramer constrained each other to achieve a fine-tuning of the reproductive success of cotton to produce fibers.

Overall, the proposed model described a sequential developmental regulation from vegetative to reproductive growth linked by the temporal and spatial expression of GhSEP genes whose function is largely dependent on the cofactors to form different complexes dynamically to target different downstream genes. It provides new insights into the sequential regulation of cotton flowering and reproductive growth that contribute to fiber production.