Identification of Morus notabilis MADS-box genes and elucidation of the roles of MnMADS33 during endodormancy

The MADS-box genes encode transcriptional regulators with various functions especially during floral development. A total of 54 putative Morus notabilis MADS-box genes (MnMADSs) were identified and phylogenetically classified as either type I (17 genes) or type II (37 genes). The detected genes included three FLOWERING LOCUS C-like (MnFLC-like) genes, MnMADS33, MnMADS50, and MnMADS7. MnFLC-like proteins could directly or indirectly repress promoter activity of the mulberry FLOWERING LOCUS T-like (MnFT) gene. Transgenic Arabidopsis thaliana overexpressing MnFLC-like genes exhibited delayed flowering and down-regulated expression of FT and SUPPRESSOR OF OVEREXPRESSION OF CONSTANS1 (SOC1). The gene expression analyses in floral bud indicated that MnMADS33 expression increased, while MnFT expression decreased during the induction of dormancy in response to cold conditions. Dormancy release resulted in the down-regulation of MnMADS33 expression and the up-regulation of MnFT expression. Furthermore, abscisic acid promoted the transcription of MnMADS33 and MnFT, although the expression level of MnFT gradually decreased. These results are consistent with the hypothesis that MnMADS33 negatively regulated the expression of MnFT to repress dormancy release and flowering in mulberry. This study may be relevant for future investigations regarding the effects of MnMADS genes on mulberry flowering development.

The MADS-box gene family is not specific to plants as it is also commonly found in animals and fungi. The MADS-box genes, which encode transcription factors, are conserved among eukaryotes 1,2 . The DNA-binding domain at the N-terminus of these transcription factors consists of a conserved region (approximately 58 amino acids) called the MADS-box domain, which is also involved in dimerization reactions and the binding of accessory factors 3 [4][5][6] . In the model plant Arabidopsis thaliana, there are two types of MADS-box genes that differ regarding their conserved sequences and structures. Type I MADS-box genes encode only the MADS-box domain, and are further divided into the Mα, Mβ, and Mγ subclasses. In addition to the MADS-box domain, type II genes encode the intervening (I), keratin-like (K), and C-terminal (C) domains that ultimately form the so-called MIKC-type domain structure 7,8 . The I and K domains affect protein-protein interactions, resulting in the formation of homodimeric or heterodimeric complexes, while the highly variable C domain influences transcriptional regulation and contributes to the formation of protein complexes 2,9 . Type II genes are classified into the MIKC c and MIKC* subfamilies 7 . Based on their sequence characteristics and functions, the MIKC c genes have been further divided into 13 subcategories, which were named after their first identified members 1,10,11 .

. The CArG-box [C(A/T)8G, C(C/T)(A/T)G(A/T)4(A/G), or C(C/T)(A/T)6(A/G)G] is a consensus MADS-box transcription factor-binding site
There is a limited amount of information regarding the A. thaliana type I MADS-box genes. In addition, the type II genes include the floral homeotic genes, whose functions are related to the determination of floral organ identities, which can be explained by the well-studied ABCDE model [12][13][14] . All ABCDE model genes belong to the MIKC c subfamily, with the exception of AP2 15 . The APETALA1/FRUITFULL (AP1/FUL), APETALA3/PISTILLATA (AP3/PI), AGAMOUS (AG), SHATTERPROOF/SEEDSTICK (SHP/STK), and SEPALLATA (SEP) subfamily genes encode proteins that interact with one another to form a complex regulatory network, that controls the development of the floral organ (e.g., sepals, petals, stamens, carpels, and ovules) 2 .
Phylogenetic analyses of the Arabidopsis thaliana and mulberry MADS-box genes. To analyze the phylogenetic relationships among MnMADS genes, an unrooted neighbor-joining (NJ) tree was constructed based on the alignment of the full-length amino acid sequences of 54 putative MnMADS proteins and 103 A. thaliana MADS-box transcription factors. A phylogenetic analysis revealed that 17 type I and 37 type II MADS-box genes encoded the 54 MnMADS proteins (Fig. 1). According to the A. thaliana MADS-box subgroup classifications, the type II proteins were further divided into the following 14 groups: SEP, AGL6, AP1/FUL, AGL12, AG, SOC1, FLC, AGL15, SVP, AGL17, ARABIDOPSIS BSISTER (ABS), PI, AP3, and MIKC*. Most of these subgroups were named according to extensively investigated MADS-box genes. To further validate the evolutionary relationships among MnMADSs, an unrooted NJ tree was constructed using the MADS-box protein sequences of A. thaliana, grape, and mulberry ( Supplementary Fig. S1). Moreover, the evolutionary relationships among the 54 MnMADS proteins were investigated using a maximum likelihood (ML) method ( Fig. 2A). It is noteworthy that the FLC subfamily was adjacent to AP1, SEP, and AGL6 in the ML tree. Furthermore, the evolutionary relationships of the MnMADS proteins in Fig. 2A were consistent with the relationships in the NJ phylogenetic tree (Fig. 1).
Mulberry MADS-box protein domains, gene structures, and classes. All MnMADS proteins consisted of a MADS-box domain. As shown in Fig. 2B, except for MnMADS7, MnMADS29, MnMADS36, and MnMADS52, all the type II MnMADS proteins contained a K-domain. Additionally, the MnMADS gene structures were investigated using comparative alignments of the corresponding genomic and coding sequences. The number of introns varied from 0 to 10, with genes in the same family containing a similar number of introns. Most of the type I genes had no introns, while the majority of type II genes had six. The genes indicated by broad red lines in Fig. 2C are manually corrected cDNA sequences. For most MADS-box gene families and subfamilies, there were more A. thaliana genes than mulberry genes (or an equal number). For example, we identified 17 type I mulberry genes (31% of the total number of MnMADS genes), while there were 58 type I A. thaliana genes (56% of the total number of A. thaliana MADS-box genes) ( Supplementary Fig. S2). In contrast, the AP3 and PI subfamilies contained more mulberry genes than A. thaliana genes. We observed that A. thaliana carried only one AP3 and one PI gene, while mulberry contained four AP3 and two PI genes. Additionally, mulberry also had four SVP/AGL24 genes, while A. thaliana consisted of only two.
Expression profiles of mulberry MADS-box genes in various organs and flower-related tissues. Analyses of the MnMADS expression profiles in the roots, leaves, bark, winter buds, and male flowers resulted in the identification of 44 MnMADS genes that were expressed in at least one of these organs, including all 37 type II genes. The expression level of seven type I genes was generally low, while the expression of the remaining 10 type I genes was undetectable in all examined organs. Overall, the MnMADS expression levels were lowest in the roots, followed by the leaves. Most genes with more than 50 RPKM were associated with the ABCDE flowering model and were preferentially expressed in winter buds and male flowers. Several genes exhibited high expression levels in various tissues. The most highly expressed gene was MnMADS44, with an expression level of 535 RPKM in male flowers. It was also highly expressed in winter buds. MnMADS46 was the next most highly expressed gene, with an expression level of 304 RPKM in male flowers. The most highly expressed gene in bark was MnMADS15 (147 RPKM). Additionally, the expression level of MnMADS33 (i.e., FLC homolog) was biased in winter buds (>80 RPKM), while the expression level of the SVP subfamily MnMADS genes was relatively low (Fig. 3, Table S3).
To further associate MnMADS functions with specific developmental processes, a quantitative real-time polymerase chain reaction (qRT-PCR) was used to analyze expression profiles in various flower-related mulberry tissues. Of the 54 MnMADS genes, 39 were expressed in at least one tissue of the mulberry species 'Jinqiang 63' (JQ63) ( Table S4). Gene expression patterns for the buds of differentiating inflorescence primordia (buds I-III), during endodormancy and dormancy release (buds IV-VII), as well as in the different catkin development stages (catkins I-III) are shown in Fig. 4 III but also expressed in buds I-III and buds IV-V. Meanwhile, MnMADS1 (i.e., SVP homolog) was abundantly expressed in bud II and III, as well as in catkin I. The mulberry genome includes four MnMADS homologs of SVP/AGL24 (i.e., MnMADS1, 4, 42, and 43). We observed that MnMADS1 was expressed at a continuously high level during the endodormancy and dormancy release periods. The remaining genes in the SVP/GAL24 family were expressed at very low expression levels during these periods.

FLOWERING LOCUS C homologs in mulberry.
To provide additional evidence for the relationship between FLC and MnMADS33/MnMADS50/MnMADS7, we reanalyzed the previously reported phylogenetic relationships among the FLC-like proteins from mulberry and other species 53 . MnMADS33, MnMADS50, and MnMADS7 were clustered in the same clade and were closely related to the Mimulus guttatus FLC homologs (Fig. 5A). An investigation of the subcellular localization in N. benthamiana leaves indicated that MnMADS50 and MnMADS7 were mainly located in the nucleus. Meanwhile, MnMADS33 localized to the epicyte and organelles including guard cells, but primarily the nucleus (Fig. 5B). Analyses of protein structures revealed that MnMADS33 and MnMADS50 contained the classic MIKC motifs. Additionally, MnMADS7 comprised the MADS-box and I-domains, but the K-and C-domains were incomplete (Figs 2 and 5C). An amino acid sequence alignment confirmed that there are 31 amino acids before the MnMADS33 MADS-box domain. Moreover, we observed a 51.08% amino acid identity between MnMADS33 and MnMADS50. In contrast, A. thaliana FLC shared only 29.00% and 37.00% amino acid sequence identity with MnMADS33 and MnMADS50, respectively (Fig. 5C). Dual luciferase assays indicated that MnFLC-like proteins (i.e., MnMADS33, 50, and 7) could directly or indirectly suppress the promoter activity of MnFT and repressed the expression of the luciferase gene (Fig. 6A). According to the yeast two-hybrid assays, the interaction between MnMADS50 and MnMADS1 is strong. In addition, a very weak interaction between MnMADS33 and MnMADS1, and no interaction between MnMADS7 and MnMADS1 were detected (Fig. 6B). β-galactosidase activity assays were

Ectopic expression of mulberry FLOWERING LOCUS C-like genes in Arabidopsis thaliana resulted in a late-bolting phenotype.
To investigate the potential functions of MnMADS33, MnMADS50, and MnMADS7 genes regarding the flowering transition, transgenic A. thaliana plants expressing these genes under the control of the CaMV 35S promoter were examined. We selected three transgenic lines for each gene and analyzed their flowering phenotypes and gene expression levels ( Fig. 7A-C). Gene-specific primers were used to conduct qRT-PCR analyses of WT and transgenic A. thaliana to validate the phenotypic assessment of the transgenic lines and the expression of flowering-related genes (Fig. 7D,E). The FT and SOC1 expression levels were lower in all transgenic lines than in the WT (Fig. 7E). Meanwhile, the transgenic lines exhibited varying degrees of delayed flowering and produced more rosette leaves than those of the the WT plants (Fig. 7F).

Relationships between endodormancy and the FLOWERING LOCUS C-and FT-like genes in mulberry.
A qRT-PCR was conducted to clarify the relationship between the seasonal periodicity of floral bud development and the expression of MnFLC-like genes and MnFT-like gene (Fig. 8A,B, Supplementary  Fig. S4). An expansion gene (MnEXP2; Morus011237) and a cyclin D-type gene (MnCYCD3; Morus027452) were used to provide evidence of the endodormancy cycle activities 54 . The MnEXP2 and MnCYCD3 expression levels were rapidly down-regulated in November and recovered during dormancy release (Fig. 8C,D). Meanwhile, the MnMADS33 and MnFT expression levels fluctuated before endodormancy, after which MnMADS33 expression was sharply up-regulated, peaking in November, and then down-regulated until dormancy release in the   Bud I, II, and III represent initial, mid-term, and later stage of differentiating inflorescence primordial, respectively. Paraffin sections were prepared to ensure the correct stages ( Supplementary Fig. S7). Bud IV and V are dormant buds. Bud VI and VII were dormancy release buds. Catkin I, catkin II, and catkin III were before, mid-term, and later stage of pollination, respectively. The data presented here represent two biological replicate experiments. The qRT-PCR is repeated three times for each biological replicate. RPL15 is used as the internal control. This heat map is created by method mentioned in  Fig. S4). An analysis of floral buds of plants grown under artificially controlled cold conditions revealed that MnMADS33 expression was rapidly up-regulated after 10 days of chilling, while MnFT exhibited the opposite expression pattern (Fig. 8E,F). We also investigated the expression levels of genes related to ABA biosynthesis and catabolism [e.g., ABA 8-hydroxylase gene (MnCYP707A4; Morus027528), RING-H2-type zinc-finger gene (MnXERICO; Morus006867), and 9-cis-epoxycarotenoid dioxygenase gene (MnNCED1; Morus012507)] in mulberry floral buds exposed to low temperatures under field and artificially controlled conditions. The genes exhibited similar expression profiles under the two experimental conditions (Supplementary Fig. S5). The MnCYP707A4 and MnNCED1 expression levels were up-regulated after endodormancy, and then down-regulated before dormancy release. The expression of MnNCED1 coincided with that of MnMADS33. Furthermore, MnXERICO expression was down-regulated during endodormancy and up-regulated before dormancy release. Mulberry floral buds were treated with ABA to confirm the relationships between ABA and the expression of MnMADS33 and MnFT (Fig. 9). The ABA treatment up-regulated the expression of MnMADS33 and MnFT; however, while the increased MnMADS33 expression level was stable, the ABA-induced MnFT expression level gradually decreased (Fig. 9A,B). The expression of MnMADS50 and MnMADS7 was up-regulated at 24 h and then down-regulated at 48 h and 96 h after ABA treatment ( Supplementary Fig. S6). Furthermore, bud break was repressed by ABA treatment (Fig. 9C).  (Table S1). The MnMADS gene structures and encoded protein domains were highly similar to those of the corresponding genes from other species 2,3,56-61 . The ABCDE model for annual plants is the most well-known floral development model. We detected AP1, AP3, PI, AG, STK, and SEP homologs among the identified 54 MnMADS genes (Figs 1 and 2, Supplementary  Fig. S2). In most cases, homologs of the ABCDE model genes produced similar expression profiles (Figs 3 and 4). For example, four SEP-like genes in the E class were expressed in a floral-related-organ-biased manner. This result was consistent with the SEP function in A. thaliana floral organ development 10,13 . The A. thaliana MADS-box gene family is regulated by a complex network and is crucial for floral transitions. We detected the A-, B-, and E-class genes in the differentiating inflorescence primordia, suggesting a functional conservation among floral signal integration pathways.

Conserved evolution of mulberry type I and type II MADS-box genes.
Despite the lack of an assembled genome, 10 clusters containing 22 genes were detected based on their locations on the scaffolds (Table S1). A previous study involving synteny network analyses indicated that AGL6-SOC1 and SEP-AP1 tandems are conserved across all angiosperms 62 . In the present study, we observed one AGL6-SOC1 and three SEP-AP1 tandems. Moreover, we detected three tandem duplication clusters, including eight Mα genes, while no duplication events were observed in the Mβ and Mγ subfamilies. This may explain why there are more Mα genes than Mβ and Mγ genes in mulberry. An analysis of the evolution of A. thaliana and Oryza sativa MADS-box genes suggested that type I genes underwent a faster birth-and-death evolution than type II genes, in part because of a higher frequency of segmental gene duplication and weaker purifying selection in type I than in type II genes 63 . Generally, Ka/Ks = 1 indicates a neutral selection, Ka/Ks <1 corresponds to a purifying selection, and Ka/Ks >1 suggests an accelerated evolution with positive selection 64 . In our study, the Ka/Ks ratios of three whole genome/segmental duplication paralogs from type II were <0.1, consistent with a strong purifying selection pressure. The paralogous pairs of the type I Mα subfamily genes had a higher Ka/Ks ratio, including two pairs (MnMADS26/MnMADS27 and MnMADS2/MnMADS3) with Ka/Ks ratios >1, suggesting the Mα genes evolved under a positive selection pressure. There were three pairs of MIKC c genes distributed in a tandem duplication pattern. According to the duplication types of the paralogous pairs (Table S2), the expansion of the type I gene group was due to tandem duplications, whereas segmental duplications played a leading role in the expansion of the type II gene group. The expansion patterns of the mulberry type I and type II MADS-box genes are similar to those reported for A. thaliana 3 , O. sativa 65 , and P. mume 58 .

Identification of mulberry FLOWERING LOCUS C-like genes. A previous study concluded that
FLC functions as a flowering repressor, which delays flowering by inhibiting the expression of the activators of flowering such as FT and SOC1 37 . A recent transcriptome study in mulberry determined that MnMADS33, MnMADS50, and MnMADS7 (i.e., c75631_g1, c78979_g1, and c79593_g1, respectively) are AGL6 homologs  (Table S1). Our data was further verified by completing a multi-species phylogenetic analysis (Fig. 5). The mulberry FT promoter has been cloned 66 . Four CArG-box like elements were identified in the MnFT promoter (Supplementary File 1). MnFLC-like proteins can directly or indirectly suppress the MnFT promoter activity and suppressed the expression of the luciferase gene in the dual luciferase assay (Fig. 6A). Moreover, the interactions of MnMADS33/MnMADS1 and MnMADS50/MnMADS1 were detected. All these results demonstrated that MnMADS33, MnMADS50, and MnMADS7 might be FLC-like genes. Although, ectopically expressed MnMADS33, MnMADS50, or MnMADS7 in A. thaliana slightly delayed flowering, the heterologous system results in a weaker phenotype than expected potentially because the binding affinity to the interacting partner or target DNA might be weaker than the original one.

MnMADS33 may down-regulate the expression of mulberry FLOWERING LOCUS T in dormant buds to inhibit dormancy release. In A. thaliana, FT is expressed in leaves and the resulting protein
migrates to the shoot apex where it induces the formation of floral primordial 67 . In P. persica, DAM5 and DAM6 expression levels are up-regulated during dormancy development and down-regulated during the winter chilling period 68 . A recent study suggested that DAM5 and DAM6 are key regulators of dormancy 35 . Similarly, we observed that MnMADS33 expression in mulberry floral buds is dramatically up-regulated after endodormancy and then down-regulated, while MnFT exhibits the opposite expression pattern (Fig. 8), suggesting that MnMADS33 affects dormancy and negatively regulates the expression of MnFT during the dormancy period.
For plants in the indirect flowering group, such as mulberry, the formation and development of the inflorescence primordia occur in the buds in summer and then cease growth during endodormancy until next bud burst in spring. Additionally, dormant mulberry buds contain floral and leaf primordia 50 . For the relatively high MnMADS33 and MnMADS50 expression levels in differentiating inflorescence primordial, we propose that the encoded proteins help to maintain the inflorescence meristems in buds and inhibit bud burst and flowering by repressing the localized expression of MnFT. Before endodormancy, the MnFT and MnMADS33 expression levels fluctuate in floral buds, and the balance between these levels likely prevents the bursting of the floral buds. MnMADS50, MnMADS7, and MnMADS33 are three FLC-like genes. However, the low expression level in floral buds during endodormancy gathered with the down-regulation by ABA indicated that MnMADS50 and MnMADS7 had little relationship to endodormancy.
A recent study suggested that P. pyrifolia PpDAM1 up-regulates PpNCED3 expression and forms a feedback regulatory loop with ABA metabolism and signaling pathways to control endodormancy 69 . In A. thaliana, ABSCISIC ACID-INSENSITIVE 4 (ABI4), which is a key component of the ABA signaling pathway, directly promotes FLC transcription and negatively regulates the floral transition 70 . In our study, MnNCED1 expression was up-regulated which is consistent with the hypothesis that high MnMADS33 expression level during the induction of endodormancy up-regulates ABA signaling and promotes endodormancy induction (Supplementary Fig. S4). Moreover, an ABA treatment up-regulated MnMADS33 and MnFT expression, although the MnFT expression level was subsequently down-regulated at 48 and 96 h after treatment (Fig. 9). Our results support the hypothesis that dormancy and flowering share overlapping pathways 46 . We also provide evidence supporting the hypothesis that FLC-like genes affect dormancy by directly or indirectly down-regulating the expression of FT-like genes. However, whether MnMADS33 and MnFT function in the same regulatory pathway or act independently during dormancy remains to be determined.

Conclusions
In summary, 54 putative MADS-box family members were identified in the mulberry genome. We investigated the MnMADS phylogenetic relationships, gene structures, and encoded protein domains. A comparison of the expression patterns of paralogous pairs provided insights into the functional conservation of MADS-box proteins from mulberry and other species. Mulberry FLC-like proteins (i.e., MnMADS33, MnMADS50, and MnMADS7) were mainly localized to the nucleus, although MnMADS33 was also detected in the epicyte and organelles.

Methods
Plant material. Mulberry JQ63 was used for gene expression analysis and vector construction. Mulberry materials were obtained from the Southwestern University mulberry germplasm for field sampling. Floral buds were collected before, during, and after the differentiation of inflorescence primordia (termed bud I, bud II, and bud III) in April 2016. Paraffin sections were prepared to ensure correct stages (Supplementary Fig. S7) 71 . Floral buds during endodormancy (November and December 2015) (as bud IV and bud V) and dormancy releasing stages (January and February 2016) (termed bud VI and bud VII) were gathered. Mulberry catkins were collected in three development stages (before, during, and after pollinated stages; termed catkin I, catkin II, and catkin III) in May 2016. Mulberry cuttings for chilling and abscisic acid (ABA) treatment were obtained on September 28, 2016 and November 4, 2016, respectively. For chilling treatment, cuttings were exposed in 4 °C and were sampled once per 10 days for total RNA extraction. For ABA treatment, cuttings were immersed in vases with 10 μM ABA (Sangon Biotech, Shanghai, China) solution in a growth chamber at 24 °C under a 14/10 h light/dark regime. After incubation for 96 h, cuttings were transferred to tap water. The control was treated with water all the time.  73 . The HMMER (version v3.1b2) 74 was used to search MADS-box genes in the mulberry database using the HMM profile. All sequences obtained using above two methods were used as queries to BLAST the de novo transcriptome assembly data 50 . The truncated sequences were manually corrected based on the de novo transcriptome assembly data. We used the GENSCAN web server (http://genes.mit.edu/GENSCAN.html) 75 to analyze the nucleotide sequences. The gene sequences were uploaded to the website to predict the coding regions, with A. thaliana as the reference organism. Furthermore, the protein sequences were analyzed on the SMART website (http://smart.embl-heidelberg.de/) 76 for domain prediction. Proteins with MADS domains were considered candidate MnMADS proteins.
Analyses of gene and protein properties. Schematic diagrams representing MADS-box gene structures (e.g., introns and exons) were generated using the Gene Structure Display Server (http://gsds.cbi.pku.edu.cn) 77 . We calculated coding sequence lengths and determined the number of amino acids, molecular weights, and isoelectric points of the encoded proteins using the Compute pI/Mw tool on the ExPASy website (http://web. expasy.org/compute_pi/) 78 . The MCScanX program 79 was used to classify duplication types of MnMADS genes based on genome sequences encoding proteins and alignments with default parameters. The ratio of the number of nonsynonymous substitutions per non-synonymous site to the number of synonymous substitutions per synonymous site (Ka/Ks ratio) of each duplicated gene pair was calculated using DnaSP (version 5.10.01) 80 on the full-length coding sequences. The STRING program (version 10.0) 81 was used to construct possible proteinprotein interaction networks based on experimental evidence with the highest confidence values (i.e., >0.900) in A. thaliana MADS-box networks. Furthermore, homologous genes encoding these interacting proteins were identified according to phylogenetic analyses. Their expression patterns in the differentiating inflorescence primordia stage were investigated by qRT-PCR.
Phylogenetic reconstruction. Full-length protein sequences of all A. thaliana, grape, and mulberry MADS-box proteins were used to analyze their phylogenetic relationships. The amino acid sequences were aligned with ClustalW and an unrooted neighbor-joining (NJ) phylogenetic tree was constructed using a bootstrapping method (1,000 replicates) in MEGA 7.0 82 with p-distance and complete deletion parameters. All MnMADS protein sequences were used to build maximum likelihood (ML) tree to ensure the reliability of the constructed tree. The unalignable regions at N and C terminals were manually trimmed to decrease phylogenetic noise. ProtTest 3 83 was used to determine the best-fit model of amino acid sequence evolution, which took the JTT + I + G model under the Akaike information criterion. We used PhyML (version 3.1) 84 to constructed ML tree with 100 bootstrap replicates. The sequence alignment analysis of A. thaliana FLC, MnMADS33, MnMADS50, and MnMADS7 was performed using the pairwise alignment program within the GeneDoc software 85 . Phylogenetic relationship of FLC homologs (Table S5) was generated by ML method using GTR + I + G model. The gene structures, protein structures, and all the phylogenetic trees were visualized using EvolView software on the website (http://www. evolgenius.info/evolview/) 86 .

Gene expression analyses.
We used reads per kilobase of exon per million mapped reads (RPKM) to compare the gene expression levels among samples (Table S3). Values for the roots, bark, winter floral buds, male flowers, and leaves of MnMADSs were obtained from RNA sequencing data (http://morus.swu.edu.cn/morusdb/), which were from single samples. A heatmap was generated based on the log2-transformed (RPKM+1) values using R project software. The expression patterns of MnMADSs during flowering were evaluated using qRT-PCR. Total RNA was extracted from the collected samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Contaminating genomic DNA in the total RNA was removed using RQ1 RNase-free DNase (Promega, Madison, WI, USA). RNA (1 μg) was used for cDNA synthesis with Moloney Murine Leukemia Virus (M-MLV) Reverse Transcriptase (Promega, Madison, WI, USA). Fourfold-diluted cDNA was used for qRT-PCR. The primers for qRT-PCR were designed with Primer Premier 5 (Table S6) and then checked by melting curves. The qRT-PCR was conducted using SYBR Premix Ex Taq II (Takara, Dalian, PR China) and the StepOnePlus Real Time PCR system (Applied Biosystems, Foster City, CA, USA), initiated by 30s at 95 °C and followed by 40 cycles of 95 °C for 5s, 60 °C for 30s, and completed with a melting-curve analysis program. The PCR mixture (20 μl total volume) comprised 10 μl of 2 × SYBR ® Premix Ex Taq, 0.4 μl of each primer (10 μM), 0.4 μl of 50 × ROX Reference Dye II, 2 μl of diluted cDNA and 6.8 μl PCR-grade H 2 O. The mulberry RPL15 gene (Morus024083) was used as an internal control gene. The relative expression of genes for each sample was calculated using the formula 2 − [Ct (target gene) − Ct (control gene)] . The log2-transformed (the relative fold changes+1) values of the 39 MnMADSs were visualized in heat maps prepared with R project software. The expression data of 39 MnMADSs were from duplicates of two biological replicates. Three biological replicates were performed for gene expression analyses in anniversary monthly floral buds, chilling-treated floral buds, and ABA-treated floral buds. The qRT-PCR analyses were performed three times for each biological replicate.

Cellular localization.
The open reading frames (ORFs) of MnMADS33, MnMADS50, and MnMADS7 without stop codons were obtained via PCR amplification using gene-specific primers (Table S7). The purified PCR products were ligated into the pLGNL 35S-GFP vector after digestion with Kpn I and BamH I enzymes. Linker sequences comprising four glycine and one serine residue with three repetitions (GGAGGAGGAGGATCAGGAGGAGGAGGATCAGGAGGAGGAGGATCA) were added between target genes and GFP. The fusion vectors were introduced into Agrobacterium tumefaciens strain GV3101 (pMP90) 87 . The A. tumefaciens cells, containing the pLGNL 35S-GFP, pLGNL 35S-MnMADS33-GFP, pLGNL 35S-MnMADS50-GFP, or pLGNL 35S-MnMADS7-GFP vectors, were grown to OD 600 = 1.0 in LB liquid medium (10 g/L tryptone, 5 g/L yeast extract, and 10 g/L NaCl) containing 35 μg/mL rifampicin, 50 μg/mL kanamycin, and 30 ng/mL acetosyringone. Cells were collected by centrifugation at 5,000 g at room temperature for 10 min and placed for 3-5 hours at 25 °C after resuspending to OD 600 = 0.4 using infiltration medium (2 g/L 2-morpholinoethanesulfonic acid, 1 g/L Magnesium chloride, and 30 ng/mL acetosyringone, pH 5.8). Infiltration was carried out by injecting N. benthamiana epidermal cells. The N. benthamiana samples were kept in darkness for 12 hours and then grown in a chamber at 24 °C under a 16/8 h light/dark regime for 2 days. The subcellular localization results were scanned by FV1200 laser scanning confocal microscope (Olympus, Tokyo, Japan).
Dual luciferase assay. The promoter of the MnFT gene was amplified with the primers described in Table S7. Full-length MnMADS33/50/7 sequences were inserted into the pGreenII 62-SK vector (SK). The promoter of MnFT was inserted into the pGreenII 0800-LUC vector 88 . N. benthamiana leaves were used in dual luciferase assays. The Agrobacterium-mediated method was performed as mentioned in cellular localization. The firefly luciferase (LUC) and Renilla luciferase (REN) were assayed using dual luciferase assay reagents (Promega, Madison, WI, USA) after infiltration for 3 days. Three biological replicates were performed (three technical replicates in each independent experiments).
Yeast two-hybrid and β-galactosidase activity assays. pGBKT7 (bait) and pGADT7 (prey) vectors were used for a yeast two-hybrid assay. Primers used for vector construction are listed in Table S7. The autoactivation and toxicity of both the bait and prey were tested. Positive controls were pGADT7-T and pGBKT7-53, and negative controls were pGBKT7-Lam and pGADT7-T. For each co-transformation pair, 200 ng of bait recombinant plasmid and 100 ng of prey recombinant plasmid were transformed into Saccharomyces cerevisiae Y2HGold, and then the cells were grown at 30 °C on synthetic medium lacking both leucine (Leu) and tryptophan (Trp) (Double Dropout, DDO). The positive clones were grown on synthetic medium without Trp, Leu, histidine (His), and adenine (Ade) (Quadruple Dropout, QDO) after diluting 1, 10, and 100 times. Three replications were conducted for each pair of proteins. The yeast strain AH109 was used for β-galactosidase activity analysis using Yeast β-galactosidase Assay Kit (Thermo Fisher Scientific, Hudson, NH, USA) following the manufacturer's protocol. Yeast transformation and culture methods were the same as mentioned above. Three biological replicates were performed for each combination.

Generation of transgenic Arabidopsis thaliana plants and expression analyses of downstream target genes.
The full-length cDNA of MnMADS33, MnMADS50, and MnMADS7 were amplified by PCR and were inserted into the Kpn I/BamH I cloning sites of a binary vector driven by the cauliflower mosaic virus (CaMV) 35S promoter and followed at the 3′ end by the nopaline synthase gene (NOS) terminator (pLGNL 35s-NOS). The primers used for plasmid construction are listed in Table S7. Overexpression of MnMADS33, MnMADS50, and MnMADS7 in A. thaliana was carried out in wild-type (WT) ecotype Columbia (Col-0) by A. tumefaciens-mediated plant transformation, which was performed by the floral dipping method 89 . Transformed seedlings were obtained from selective medium containing 50 mg/L kanamycin, and then transferred to soil. Homozygous lines of transgenic A. thaliana plants from T3 generation were used for further analyses. Plants were grown in chamber at 24 °C under a 16/8 h light/dark regime. Flowering time and the number of rosette leaves were recorded when the primary inflorescence was 0.5 cm long 90 . To understand the function of mulberry FLC-like genes in transgenic A. thaliana, the expressions of the flowering time genes, FT and SOC1, were investigated using qRT-PCR analyses. A. thaliana Actin7 (At5g09810) was used as the internal control to normalize the expression levels. Conditions for total RNA isolation, reverse transcription, and qRT-PCR are the same as mentioned above. Primers used for expression analyses of downstream target genes in transgenic plants were listed in Table S6.