Genome-wide identi ﬁ cation and expression analysis of MADS-box transcription factors reveal their involvement in sex determination of hardy rubber tree ( Eucommia ulmoides oliv

Eucommia ulmoides is a famous rubber-producing and medicinal tree species that produces unisexual ﬂ owers on separate individuals from the earliest stage of stamen/pistil primordium formation. To explore the genetic regulation pathway of sex in E. ulmoides , comprehensive genome-wide analyses and tissue-/sex-speci ﬁ c transcriptome comparisons of MADS-box transcription factors were performed for the ﬁ rst time in this work. Quantitative real-time PCR technique was employed to further validate the expression of genes that are assigned to ﬂ oral organ ABCDE model. A total of 66 non-redundant E. ulmoides MADS-box (EuMADS) genes were identi ﬁ ed, they were classi ﬁ ed into Type I (M-type, 17 genes) and Type II (MIKC, 49 genes). Complex protein-motif composition, exon-intron structure and phytohormone-response cis-elements were detected in MIKC-EuMADS genes. Furthermore, 24 differentially-expressed EuMADS genes (DEGs) between male and female ﬂ owers, and two DEGs between male and female leaves were revealed. Amongst the 14 ﬂ oral organ ABCDE model-related genes, there were 6 (A/B/C/E-class) and 5 (A/D/E-class) genes displayed male-and female-biased expression respectively. In particular, one B-class gene EuMADS39 and one A-class gene EuMADS65 were almost exclusively expressed in male trees, no matter in ﬂ ower or leaf tissues. Collectively, these results suggested a critical role of MADS-box transcription factors in sex determination of E. ulmoides , which is conducive to decoding the molecular regulation mechanism of sex in E. ulmoides .

Recent studies revealed that plant sex is highly associated with MADS-box genes . For example, LcMADS42/46/ 47/51/75/93/100 are probably involved in the unisexual flower formation of litchi (Litchi chinensis Sonn.) (Guan et al., 2021). One MADS-box gene QsPISTILLATA is also exclusively expressed in male flowers of Quercus suber L., and is proved to be functional for stamen development (Sobral and Costa, 2017). Remarkably, B-class MADS-box TFs are essential for sex determination of several dioecious plants . In spinach (Spinacia oleracea L.), B-class genes SpAP3 and SpPI are expressed specifically in the floral primordial of male. Suppression of SpAP3 or SpPI turned male individuals into female (Pfent et al., 2005;Sather et al., 2010). Similarly, in Thalictrum dioicum L., all B-class genes (ThdAP3-1, ThdAP3-2a, ThdAP3-2b, ThdPI-1 and ThdPI-2) are only expressed in floral buds of male but not female at very initial stage. Knocking-down the expression of ThdPI-1 or ThdPI-2 also converted male individuals into female (Di Stilio et al., 2005;Larue et al., 2013).
Eucommia ulmoides Oliv. is known as the hardy rubber tree that can synthesize gutta-rubber (trans-1,4-polyisoprene, TPI) almost in the whole plant, especially in fruits . It is also the resource of well-known traditional Chinese medicine 'DuZhong' (Ouyang et al., 2021). Meanwhile, this species is a strict dioecious plant with stamen primordium or pistil primordium initiating in separate individuals, suggesting that its flowers are unisexual from inception Qing et al., 2021). Some MADS-box TFs, especially them involved in 'ABCDE model' of flower development, have been suggested to take part in sex determination of dioecious plants with such type of unisexual flowers (Diggle et al., 2011;Zhang et al., 2022). However, the characteristics and roles of MADS-box TFs in sex determination of E. ulmoides remain largely unknown by far.
In this work, MADS-box TFs of E. ulmoides were identified throughout the genome and comprehensively characterized, e.g., their evolutionary relationships, protein motif compositions, gene structures, and phytohormone response cis-elements. Tissue-and sex-specific expression profiling of all the MADS-box genes were investigated based on transcriptome data. Quantitative real-time PCR (qRT-PCR) analyses were further conducted to uncover the sex-biased expression patterns of A/B/C/D/E-class genes involved in flower development. The results obtained herein will on the one hand shed lights on our understanding of the evolution and function of MADS-box TFs in E. ulmoides, on the other, will promote decoding the molecular regulation mechanism of sex in this dioecious plant.

Genome-wide identification and synteny analysis of MADS-Box TFs in E. ulmoides
A total of 66 non-redundant MADS-box TFs were identified in E. ulmoides for the first time according to BLASTP and HMMER analyses, named as EuMADS01-66 (Table 1). Wherein 64 of them were derived from the haploid E. ulmoide genome , the rest two were discovered in the male genome , with 58 members being commonly detected from both the haploid and male genomes. The identified 66 EuMADS TFs were further verified in SMART and CDD databases for the presence of conserved MADS domain. These MADS-box TFs were sized from 95 aa (EuMADS04) to 471 aa (EuMADS62) with an average of 232 aa. Their isoelectric points (pl) varied from 4.7 (EuMADS02) to 11.5 (EuMADS05), molecular weight (MW) were between 10.7 kDa (EuMADS04) and 53.4 kDa (EuMADS62), number of phosphorylation sites (Ser, Tyr and Thr sites) were from 11 (EuMADS04) to 72 (EuMADS62). As expected, subcellular locations of all the 66 EuMADS TFs were predicted to be in the nucleus. Information of the nucleotide and amino-acid sequences of these 66 EuMADS genes were documented in Additional Data1.
The 66 EuMADS genes were then mapped to the assembled 17 chromosomes of haploid E. ulmoides   (Figure 1A). As a result, a highly uneven distribution pattern of EuMADS genes was revealed, with Chromosome 8 (Chr8) harboring the maximum number of EuMADS genes (11, 16.7%), followed by Chr9 and Chr15 (both contain 6, 9.1%), and Chr4 and Chr17 (both contain 5, 7.6%). While only one was localized on Chr5, none was seen on Chr3, Chr11 and Chr14 ( Figure 1B). Besides, five EuMADS genes were unable to be assigned to any of the 17 assembled chromosomes (Li Frontiers in Genetics frontiersin.org    . #TM8 gene s present in Solanum lycopersicum (Wang et al., 2019) but not in Arabidopsis thaliana (Parenicová et al., 2003).  Table S1). Meanwhile, synteny analyses of MADS-box genes between Eucommia and Arabidopsis/rice showed more collinear regions with Arabidopsis (27) than that with rice (11) ( Figure 2B). This may suggest closer evolutionary relationship between the two eudicot plants, i.e., Eucommia and Arabidopsis.

Phylogenetic relationship and classification of EuMADS TFs
In order to correctly classify the EuMADS genes, phylogenetic relationship among the EuMADS TFs (66) and MADS-box TFs of A. thaliana (AtMADS, 103) were analyzed. Maximum likelihood (ML) tree showed that 17 EuMADS TFs were clustered in Type I (M-type) subfamily, while the rest 49 belonged to Type II (MIKC) (Figure 3). In M-type subfamily, nine EuMADS TFs were clustered in Mα lineage while the remaining eight genes were distributed in Mγ lineage. The Mβ lineage was absent in E. ulmoides.

Protein domain, gene structure, and promotor cis-element of EuMADS TFs
To obtain more information about functions of EuMADS TFs, the protein domain/motif composition, exon-intron organization, and cis-elements in promotor region were analyzed for each EuMADS gene. The conserved domains of 66 EuMADS TFs were identified using Batch CD-Search in NCBI. In summary, EuMADS TFs contained seven types of conserved domains in two categories,  Frontiers in Genetics frontiersin.org between the two EuMADS subfamilies ( Figure 5C). Motif 1 was the most conserved motif in the EuMADS TF family and distributed in the vast majority of members (60/66, 90.9%). Motif 2 was present in most MIKC-EuMADS members (41/49, 83.7%), but not found in any M-type EuMADS TFs. Interestingly, motif composition of EuMADS TFs from the same phylogenetic clade was similar, indicating possible resemblance in biological functions. For example, protein motifs in the AP3/PI subgroup (B-class) all arranged conservatively as 1-3-2, while the AG/SHP subgroup (C-/D-class) all usually arrayed as 1-3-7-2-4.
Given that some EuMADS TFs may participate together in the same pathway, e.g., floral organ development, protein-interactive analysis was performed by STRING. The results showed that 30 out of 66 EuMADS proteins could interact with each other ( Figure 8A, Supplementary Table S4). Among these putatively interacting proteins, 10 were involved in the floral organ 'ABCDE model' ( Figure 8B). The two B-class genes (EuMADS28 and EuMADS39) were shown to interact with each other, and both of them could interact with the C-class gene (EuMADS13). Some E-class genes, e.g., EuMADS64 and EuMADS37 could interact with all the members of B/C/D-class (EuMADS06, 10, 13, 28, and 39), and two A-class genes (EuMADS65 and EuMADS66).

Sex-biased expression pattern of floral organ ABCDE model-related genes
All the floral organ ABCDE model-related genes that were revealed active in flowers by transcriptome data (Figure 7) were further selected for qRT-PCR analysis. The result revealed most of  Table S3).
Remarkably, the expression patterns of most A/B/C/D/E-class genes between the male and female flowers were in accordance with their identities in determining floral organs. For example, B-class genes, i.e., EuMADS28 and EuMADS39 both highly expressed in the male flowers, but hardly expressed in the female flowers. Two D-class genes, EuMADS06 and EuMADS10 both showed female-biased expression patterns. One A-class gene (EuMADS66) and two E-class genes (EuMADS03 and EuMADS33) displayed similar expression level between the male and female flowers. Nevertheless, there were some exceptions, for instance, EuMADS49 and EuMADS65 of A-class were highly expressed in male flowers, while their closely relative EuMADS38 expressed preferably in female. C-class gene EuMADS13 was expressed in both male and female flowers, but with a higher expression level in male. In E-class, EuMADS34 showed malebiased expression, while EuMADS37 and EuMADS64 had malebiased expressions.
In addition, two B-class genes (EuMADS28 and EuMADS39) were also actively transcribed in the male leaves ( Figure 6B). EuMADS39 that was homologous to AP3 of Arabidopsis and DEF of snapdragon (Antirrhinum majus L.) (Table 2) showed a male-biased expression pattern consistently in both flowers and leaves (Supplementary Tables S1; Figures 6, 7). Notably, EuMADS65 of A-class that shared homology with gene FUL in Arabidopsis and gene DEFH28 in snapdragon was also expressed specifically in both flower and leaf of male individuals. These two genes were therefore most probably involved in sex determination of E. ulmoides.

Genome-wide characterization of MADS-box TFs in E. ulmoides
MADS-box TFs play foundational and indispensable roles in floral organogenesis and flowering of plants (Theien, 2001;Theissen et al., 2018). In present work, 66 MADS-box TFs of E. ulmoides were identified (Table 1), they were unevenly distributed across the genome (mainly on 14 chromosomes, Figure 1). The number of MADS-box genes identified in present work (66) is fewer than that in previous study (100) (Wang and Zhang, 2017). This inconsistency may be largely attributed to the alternative splicing of genes (Parenicová et al., 2003;Arora et al., 2007). In E. ulmoides, 49 MADS-box genes of the Type II (MIKC) subfamily averagely contained 7.8 exons per gene (Table 1; Figure 5), which likely would generate >49 transcripts due to different exon combinants during mRNA processing. Under such circumstance, prediction of more putative EuMADS unigenes solely based on transcriptome data is probable. For example, two putative MADS-box unigenes (Cluster- 47702.5450 and Cluster-47702.5456) predicted previously (Wang and Zhang, 2017) were practically derived from the same gene locus of EuMADS65 (Supplementary Table S3) in present work. In addition, the quality of data assembling (both transcriptome and genome) and approaches employed are certainly essential for MADS-box family screening (Wang et al., 2019;Guan et al., 2021). The MADS-box TF family has been studied in certain species with bisexual flowers, e.g., tomato (137 genes) (Wang et al., 2019), foxtail millet (89 genes) (Lai et al., 2022) and American beautyberry (78 genes) (Alhindi and Al-Abdallat, 2021), and a few species with unisexual flowers, such as litchi (L. chinensis, 101 genes) (Guan et al., 2021), hop (H. lupulus L., 65 genes) (Gutierrez et al., 2022) and Populus trichocarpa Torr. & A. Gray ex Hook (105 genes) (Leseberg et al., 2006). Interestingly, the size of MADS-box TF family identified in H. lupulus (65) (Gutierrez et al., 2022) is very close to the numbers of EuMADS TFs we found here. Neither male nor female flowers in E. ulmoides comprise sepal and petal, only have bracts encompassing 8-12 stamens (male) or one pistil (female) . Likewise, in H. lupulus the female flowers entirely lack sepal and petal, while male flowers comprise five sepals with five stamens inside (Shephard et al., 2000). Their resembling dioecious sexual system and unisexual floral composition may reflect similar family size of MADS-box TFs in E. ulmoides and H. lupulus, or vice versa.

Essential position of floral organ ABCDE model-related genes in sex determination pathway of E. ulmoides
Sex determination of dioecious plants are usually reflected in floral phenotype, which are related with the expression of floral organ ABCDE model-related genes (Cronk and Muller, 2020;Zhang et al., 2022). In fig (Ficus hispida L.), a MADS-box gene FhAG2 was found locating at the sex determination region of Y chromosome as potential male activator . Likewise, studies on cycad (Cycas panzhihuaensis L. Zhou and S. Y. Yang) (Liu et al., 2022), papaya (Carica papaya L.) (Yu et al., 2007), ginkgo (Ginkgo biloba L.) (Liao et al., 2020;Gong and Filatov, 2022), Nepenthes pitcher species (Scharmann et al., 2019), and Silene latifolia Poir. (Matsunaga et al., 2003). also suggested floral MADS-box genes were putative sex determination genes to exhibit sex-linked inheritance. In P. tremula L., ARR17 was a feminizing factor (Müller et al., 2020), CRISPR/Cas9-induced arr17 mutants upregulated B-class gene popPI to turn female individuals into male (Leite Montalvão et al., 2022). Similarly, in Diospyros lotus L., B-class gene DlPI was negatively controlled by the feminizing gene MeGI and a regulation cascade of OGI-MeGI-DlSVP-DlPI was proposed for sex determination in persimmon (Akagi et al., 2014;Yang et al., 2019).
Herein we found that 24 out of the 66 EuMADS genes were DEGs between the male and female flowers, meanwhile two EuMADS genes were also DEGs between the male and female leaves (Figure 7, Supplementary Table S3). It is worth noting that most (11/14, 78.6%) of the floral organ ABCDE model-related genes showed sex-biased expression patterns in E. ulmoides (Figures 7A,9). Unisexual flowers can be produced when the spatiotemporal expression domain of B-and C-class MADS-box genes are overlapped (Irish, 2017;Jabbour et al., 2022). As expected, two B-class genes we identified in E. ulmoides, i.e., EuMADS28 and EuMADS39 (Table 1; Figure 4) were both highly expressed in the male flowers based on transcriptome analysis and qRT-PCR validation ( Figures 7A, 9). Protein interactive analysis also predicted protein-protein interactions among B/C/E-class genes (Figure 8). In particular, EuMADS39 was expressed specifically in the male leaves too ( Figure 7B), suggesting its constitutiveexpression pattern in male individuals. This gene was previously described as EuAP3 (Wang and Zhang, 2017) by our group or EuDEF (Zhu, 2019) by another team, consistent with its functional annotation here (Table 2). More importantly, EuMADS39 was also physically adjacent to a male-linked molecular marker (MSL4) on Chr10  and exhibited a persistently male-biased expression pattern during the whole period of floral bud differentiation (Qing et al., 2021). All these above clues lead us to hypothesize that B-class EuMADS genes, particularly EuMADS39,

Class
Gene Name Ortholog in 1 Arabidopsis / 2 snapdragon ( 3 petunia)  Frontiers in Genetics frontiersin.org seem to participate in sex determination pathway of E. ulmoides as (in) direct players. Gene function verification of EuMADS39 is urgently needed in the future. According to the classical 'ABCDE model', we know that A/C/ D/E-class MADS-box genes direct the organogenesis of (sepal and petal)/(stamen and carpel)/ovule/the whole floral organ, respectively (Theien, 2001;Theissen and Saedler, 2001;Soltis et al., 2007). As expected, D-class genes (EuMADS06 and EuMADS10) both showed female-biased expression pattern ( Figures 7A, 9), in line with its role in ovary development. However, expression patterns of certain A/C/ E-class EuMADS genes unfit their identity as described in 'ABCDE model'. For example, EuMADS65 (A-class gene), homologous to FUL in Arabidopsis and DEFH28 in snapdragon (Table 2), was specifically expressed in male individuals, both in the flowers and leaves (Figures 7, 9). In A. thaliana, FUL is paralog of AP1 and plays multiply functions in controlling flowering time, meristem differentiation, and flower development, etc. (Ng and Yanofsky, 2001;Parenicová et al., 2003;Soltis et al., 2007). There are different flower arrangements between male and female individuals of E. ulmoides, with the male flowers forming a capitulum inflorescence while the female flowers being solitary on branch Qing et al., 2021), similar to that in persimmon (Akagi et al., 2014;Yang et al., 2019) and kiwifruit (Akagi et al., 2018;Varkonyi-Gasic et al., 2021). It is thus possible for A-class EuMADS65 functioning in determining the male flower inflorescence.

Putative Function
Notably, neither sepals nor petals were present in E. ulmoides female flowers Qing et al., 2021). Given that A + B + E genes determine the petal identity (Theien, 2001;Theissen and Saedler, 2001;Soltis et al., 2007), the silence of B-class genes (EuMADS28 and EuMADS39) in E. ulmoides female flowers ( Figures 7A, 9) may partly contribute to the absence of petals. Meanwhile, the four A-class genes (EuMADS38, EuMADS49, EuMADS65 and EuMADS66) showed varied expression patterns between the male and female flowers of E. ulmoides (Figures 7A,9) and different protein-protein interactions (Figure 8), suggesting their functional differentiation during evolution. This may lead to the loss of ancestral function of A-class genes in determining sepal and petal identity, resuling in the absence of both sepals and petals in E. ulmoides female flowers. Additionally, the number of stamens (8-12) in each male flower is much larger than that of carpel (only 1) in each female flower of E. ulmoides (Zhu, 2019;Wang et al., 2020), as in poplar (Xue et al., 2020) and willow . This probably leades to higher expression level of C-class gene (EuMADS13) in male flowers than that in female flowers ( Figures 7A, 9). Further functional analysis of A-class genes (e.g., EuMADS65) and the C-class gene (EuMADS13) is required to understand their roles in sex determination pathway of E. ulmoides.
Based on GFF3 file of genome annotation (https://www.ncbi. nlm.nih.gov/data-hub/genome/GCA_016647705.1/), the physical positions of 66 EuMADS genes were mapped against the assembled 17 chromosomes of E. ulmoides . TBtools software  was used for visualization of genes on chromosomes. The number of EuMADS genes per chromosome was then counted and plotted. Moreover, the information of length and gene-density of 17 chromosomes of E. ulmoides were obtained from the reference genome . The multiple collinearity scan toolkit X (MCScanX) program (Wang et al., 2012) was performed with default parameters for collinearity analysis among the 66 EuMADS genes. Then the circos map was drawn using TBtools software . Similarly, the homology of these EuMADS genes among E. ulmoides, A. thaliana, and O. sativa were also analyzed using MCScanX (Wang et al., 2012). The synteny genome blocks between species were then visualized in TBtools .

Protein domain, gene structure, and ciselement analyses of EuMADS TFs
Conserved protein domains of the 66 EuMADS TFs were blasted in the NCBI-CDD website (https://www.ncbi.nlm.nih.gov/Structure/ bwrpsb/bwrpsb.cgi) (Marchler-Bauer et al., 2015). The online Multiple Em for Motif Elicitation (MEME) program (http://meme-suite.org/tools/ meme) (Bailey et al., 2009) was also employed to analyze the protein motifs of EuMADS TFs. Parameters for MEME search were set as follows: maximum number of motifs = 20, motif width = 6-60, number of repetitions = 0/1, according to our knowledge of plant MADS-box proteins (Wang et al., 2019;Alhindi and Al-Abdallat, 2021;Dong et al., 2021;Gutierrez et al., 2022). Moreover, exon-intron structures of the EuMADS genes were constructed based on the alignment between the full-length coding sequences (CDS) and the corresponding genomic sequences Li et al., 2020). The annotation information for each EuMADS gene was also extracted from the GFF3 files (PRJNA599775 and PRJCA000677) of the genome data to verify the constructed gene structure. The protein domains, motif composition and gene structure of EuMADS TFs were then visualized using TBtools .
The 2,000 bp promotor sequences upstream the initiation codon (ATG) of each EuMADS gene were extracted from the genome data Li et al., 2020) in TBtools . The online PlantCARE tool (http://bioinformatics.psb.ugent.be/webtools/ plantcare/html/) was then performed for the cis-acting regulatory element (CARE) prediction. The cis-elements related to phytohormone responsiveness, defense and stress responsiveness, light responsiveness, low-temperature responsiveness, wound responsiveness, anaerobic induction, and circadian control were analyzed. In the phytohormone-response cis-elements, ABRE was involved in ABAresponsiveness, AuxRR-core and TGA-element were involved in IAAresponsiveness, O2-site was involved in zein-metabolism, GARE-motif, TATC-box and P-box were involved in GA-responsiveness, TGACGmotif and CGTCA-motif were involved in MeJA-responsiveness, and TCA-element was involved in SA-responsiveness (Alhindi and Al-Abdallat, 2021;Ye et al., 2022) The location of cis-elements in promoter region of each EuMADS gene was visualized in TBtools . Cis-elements associated with phytohormones (ABA, IAA, zein, GA, MeJA, and SA) were also counted and plotted in TBtools .

Tissue-and sex-specific expression analysis of EuMADS TFs by RNA-seq and protein-interaction prediction
To reveal the expression patterns of MADS-box genes in both vegetative and reproductive tissues of E. ulmoides, leaves and flowers from multiple individuals of male and female were sampled for transcriptome analysis. The RNA-seq data of male and female leaves with three biological replicates were from our previous work (Wang and Zhang, 2017). Male and female flowers from the same corresponding trees that offered leaf tissues were collected in April 2021. Flower samples were named as Eu_Male1, Eu_Male2, Eu_Male3 and Eu_Female1, Eu_Female2, Eu_Female3. Total RNA was extracted using RNeasy Plant Mini Kit (74904, Qiagen, German). mRNA library construction and Illumina sequencing were conducted according to manufacturer's procedures at BGI Technologies Corporation (Shenzhen, China).
About 6 Gb clean reads were generated for each sample after quality control. These reads were then used to quantify the expression level of each EuMADS gene using RSEM software (Li and Dewey, 2011). FPKM (fragments per kilobase of transcript per million mapped) values of EuMADS genes were calculated from the uniquely mapped reads (Zhao et al., 2021). For each tissue, i.e., flower and leaf, reads from male and female individuals were batched into one dataset, respectively for differential expression analysis with DESeq R package v1.10.1 (Anders, 2010). Benjamini and Hochberg's approach was employed to adjust the resulting p-values for false discovery rate (FDR) control (Benjamini and Hochberg, 1995). Genes with log 2 (fold change value) ≥ 1 or ≤ −1, and padj < 0.05 were considered as differentially expressed genes (DEGs) (Wang and Zhang, 2017;Guan et al., 2021), i.e., male-or female-biased expression. Flower RNA-seq data newly generated in this study were deposited in SRA database (accession number: PRJNA399774).
Given that MADS-box proteins often form dimers or polymers with each other to play roles in plant growth and development (Theissen and Saedler, 2001;Theißen et al., 2016;Ruelens et al., 2017). The protein sequences of 66 EuMADS TFs were submitted to the STRING website (https://cn.string-db.org/) (Szklarczyk et al., 2015) for protein interactions prediction. Then the protein-protein interaction network was exported and beautified by the Adobe Illustrator tool (McLean, 2002). The protein interactions of the identified 14 A/B/C/D/E-class EuMADS genes was also predicted in the same way.  (Wang and Zhang, 2017;Wuyun et al., 2018;Zhu, 2019), only stamens and pistils were isolated and immediately immersed into liquid nitrogen before stored at −80°C until use. Total RNA extraction and cDNA synthesization were performed using RNeasy Plant Mini Kit (74904, Qiagen, German) and SuperScript ™ IV VILO ™ Master Mix (Thermo Fisher, United States of) respectively following the manufacturer's instructions.
qRT-PCR analysis was conducted according to the manuals of Hieff qRT-PCR SYBR Green Master Mix (YEASEN, Shanghai, China) on a LightCycler 480 II Real-Time PCR Platform (Roche, Germany). The reaction system (20 μL) consisted of 10.0 μL of qRT-PCR Mix, 1.0 μL cDNA, 0.4 μL forward primer, 0.4 μL reverse primer, and 8.2 μL of ddH 2 O. EuGAPDH gene was used as internal control for data normalization . 2 −ΔΔCT method was applied to calculate the relative gene expression of 14 A/ B/C/D/E-class EuMADS genes (Livak and Schmittgen, 2001). Primers used for EuMADS and EuGAPDH genes (Supplementary Table S6) were designed using the online program Primer3 (http:// bioinfo.ut.ee/primer3-0.4.0/primer3/). Five biological replicates and two technical replicates were carried on for each gene. Differences of gene expression level between male and female flowers were analyzed via ANOVA followed by Student's t-test in SPSS software (v.24, IBM).
The identified 14 A/B/C/D/E-class EuMADS genes were also functionally annotated by the online BLAST in NCBI database (https://blast.ncbi.nlm.nih.gov/Blast.cgi). The most similar genes in model plants (Arabidopsis, snapdragon/petunia) of flower development studies (Theißen et al., 2016;Irish, 2017;Ruelens et al., 2017) were screened as orthologous genes to predict gene function.

Conclusion
In this study, genome-wide comprehensive analyses of the MADS-box TF family in E. ulmoides were performed for the first time. A total of 66 EuMADS TFs were identified, including 49 Type II (MIKC) members and 17 Type I (M-type) members. All the previously recorded 13 MIKC C subgroups were detected in E. ulmoides, but the Mβ lineage was missing. More complicated protein-motif composition, exon-intron architecture, and phytohormone-response cis-elements were found in MICK-EuMADS TFs, suggesting more diverse biological functions in these genes. Furthermore, tissue-and sex-specific transcriptome analyses revealed 24 DEGs between male and female flowers, and two DEGs between male and female leaves. qRT-PCR validation of the 14 EuMADS genes involved in the floral organ 'ABCDE model' showed six male-biased A/B/C/E-class genes, five female-biased A/D/E-class genes and three non-biased A/E-class genes in E. ulmoides flowers. Notably, the B-class gene EuMADS39 (ortholog of AtAP3) and the A-class gene EuMADS65 (ortholog of AtFUL) were significantly expressed in the male individuals, no matter in flower or leaf tissues. In short, these results suggested essential roles of floral organ ABCDE model-related MADS-box TFs in sex determination of E. ulmoides, as demonstrated in dioecious poplar (Leite Montalvão et al., 2022), persimmon (Yang et al., 2019) and kiwifruit (Ye et al., 2022) recently.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions
XZ: Data analysis and original draft; XW, LP, WG, and YL: sample collection and data analysis; WW: Inputs and revision. All authors read and approved the submitted version.

Funding
Rural Revitalization Project of Guangdong Province (KTP20210181) and Innovation Project for Forestry Science and Technology in Guangdong (2021KJCX015).