Genome-Wide Identification and Characterization of MADS-box Family Genes Related to Floral Organ Development and Stress Resistance in Hevea brasiliensis Müll . Arg .

Elucidating the genetic mechanisms associated with the transition from the vegetative to reproductive phase in the rubber tree has great importance for both theoretical guidance and practical application to yield genetic improvement. At present, many transcription factors, including those that belong to the MADS-box gene family, have been revealed to have roles in regulating the transition from vegetative growth to reproductive growth. However, to the best of our knowledge, the Mad-box gene family from H. brasiliensis Müll. Arg. has not been characterized in detail. To investigate members of the HbMADS-box gene family associated with floral organ and inflorescence development in H. brasiliensis, we performed genome-wide identification and analysis of the MADS-box gene family related to flower development in H. brasiliensis, and a total of 20 MADS-box genes were newly identified in the H. brasiliensis genome. Expression profiling revealed that HbMad-box genes were differentially expressed in various tissues, which indicated that HbMad-box genes may exert different functions throughout the life cycle. Additionally, 12 genes (HbSEP, HbAGL9.1, HbAGL9.2, HbCMB1, HbCMB1-L, HbAGL6, HbAGL8, HbAP1, HbAG, HbDEFL, HbTT16, and HbPADS2) were found to be associated with the differentiation of flower buds and may be involved in flower development in H. brasiliensis. All of these floral-enriched HbMADS-box genes were regulated by hormone, salt, cold, high-temperature, and drought stresses. The present study is the first to carry out the genome-wide identification and analysis of the MADS-box gene family related to flower development in H. brasiliensis, and 20 new HbMad-box genes were identified in H. brasiliensis. Most of the newly identified HbMad-box genes were found to be associated with the differentiation of flower buds and may be involved in flower development in H. brasiliensis. Our results demonstrated that HbMad-box genes may be multifunctional regulators that have roles in distinct aspects of development, and are mainly involved in the maintenance of floral organ and inflorescence development.


Introduction
Natural rubber (NR) is an important industrial and strategic raw material, and has been applied to many aspects of social production [1].Although more than 2000 plant species in the world are considered to be latex producers, the rubber tree (Hevea brasiliensis Müll.Arg.) is the only commercial source of NR because of its high yield and the excellent physical properties of its rubber products [2], and it supplied 92% of the 10.2 million tons of NR consumed worldwide in 2016 [3].
With the rapid development of the world economy, the consumption of NR in major industrial countries is increasing year by year.To meet the growing demand for NR, it is necessary to expand the planting area of rubber trees.However, the rubber tree originates from the Amazon rainforest, so the planting area must be located in sub-tropical to tropical zones [4].At present, 95% of rubber trees in the world are now mainly cultivated in South-East Asia, so the regions suitable for planting rubber trees are very limited.Therefore, there is an urgent need to improve the rubber yield per area.
To the best of our knowledge, breeding new varieties is one of the most effective approaches to increasing rubber yield per hectare.However, breeding experiments to yield genetic improvement of rubber trees are very inefficient and time-consuming, mainly because of the rubber tree's long life cycle of more than 30-35 years; it is immature for five to eight years until the rubber tree reaches the age of commercial productivity [5,6], and takes more than three decades to breed and select new clones for commercial production [7].Furthermore, due to the low rate of success for controlled pollinations, genetic improvements of H. brasiliensis are very difficult and slow [8].Therefore, research on the genetic mechanisms that affect the transition from the vegetative to reproductive phase in H. brasiliensis can provide insight for producing advantageous genetic improvement methods for controlling reproduction by genetic engineering and accelerate Hevea breeding.
In flowering plants, the transition from vegetative growth to reproductive growth is an important developmental process that involves many gene regulatory processes [9].To date, many researchers have attempted to elucidate the functional genes in association with the transition from the vegetative growth to reproductive growth of plants.It is worth noting that many transcription factors (TFs), including those that belong to the MADS-box gene family, have been demonstrated to have roles in regulating the transition from vegetative growth to reproductive growth [10].However, to the best our knowledge, the rubber tree Mad-box gene family has not been characterized in detail.
As a floral homeotic gene family, the Mad-box gene family was previously identified and investigated in the model plants Arabidopsis Heynh. in Holl & Heynh.and Nicotiana tabaccum L. [11][12][13], and has evolutionarily conserved DNA-binding domains, called the MADS-box [14].Typically, the MADS-box protein sequences can be divided into four characteristic domains from the N to the C terminus: the MADS-box (M), intervening (I), keratin-like (K), and C-terminal (C) domains [15].In plants, based on the structural features, MADS-box TFs usually contain two main groups-type I (M-type) and type II (MIKC-type) genes [16]; whereas, the type II genes can be categorized into MIKCc-and MIKC*-type [17].
In previous studies, it has been reported that this superfamily encodes transcriptional regulators that are involved in various processes, including floral organ development [18,19], root development [20], leaf development [21], fruit development and maturation [22][23][24], and embryonic development [25,26].In addition to growth and development-related functions, some MADS-box genes also play important roles in response to stress stimuli [27,28].For instance, MADS-box genes have already been proved to play important roles under low temperature stress in tomato plants [29], while several MADS-box genes have been demonstrated to take part in cold, salt, and drought responses in rice [30].Furthermore, a few MADS-box genes have been shown to be affected by the application of hormones and they exhibited differential expression in response to cytokinin, gibberellin [31], ethylene [32], and auxin [33] application in other plants.
Despite the fact that the MADS-box gene plays a great role in plant growth and development, only a few MADS-box genes have been identified and characterized in the rubber tree to date [34].For example, previous studies have found that HbAGL62 is a specific expression in flowers and embryos, and highly expressed in the flower bud differentiation stage, which indicated that HbAGL62 might play an important role in flowering regulation of the rubber tree [35]; MADS27 is highly expressed in the flower buds of rubber trees and may participate in the flowering process of rubber trees [36]; HbMADS1 and HbMADS3 have highly frequent transcriptions in the laticifer cells and somatic embryogenesis, and their transcriptions are induced in the laticifer cells by jamonic acid, which indicates that HbMADS1 and HbMADS3 may be important in natural rubber biosynthesis and somatic embryogenesis in the rubber tree [37].In the present research, we newly identified 20 MADS-box genes in the rubber tree, and analyzed their gene structure and phylogenetic relationships.To identify differentially expressed patterns of Mad-box genes in various tissues, the 20 MADS-box genes of H. brasiliensis were detected and analyzed using real-time quantitative PCR (RT-qPCR).Furthermore, to understand the responses of HbMADS-box genes to various stresses, the expression profiles of 12 floral organ-specific HbMADS-box genes were examined in leaves of Hevea seedlings after hormone, salt, cold, high-temperature, and drought stress treatments.

Plant Materials and Treatments
12-year-old rubber tree clones CATAS 7-33-97 were used as the experimental material in this study.The rubber trees were grown under normal field conditions at the Experimental Station of the Rubber Research Institute, the Chinese Academy of Tropical Agricultural Sciences (Danzhou, Hainan, China).The fresh tissues and organs (including roots, stem, stem tips, leaves, labeled bark, xylem, latex, fruits, inflorescence, male flowers, and female flowers) were collected from 12-year-old mature trees of CATAS 7-33-97 during the flowering period.Each sample was harvested on average from five trees, and three biological replicates were taken from each sample (the image of some samples shown in Figure S1).Then, the prepared samples were frozen in liquid nitrogen and transferred to a −80 • C refrigerator for RNA separation.
The tissue culture seedlings of CATAS 7-33-97 were treated with cold, high temperature, and drought stress.Each treatment was set up with three replicates, and each replicate consisted of five seedlings.Under the cold stress condition, the seedlings were grown in a culture incubator set at 5 • C and continuous illumination.For high-temperature stress treatment, the tissue-cultured seedlings were planted at 40 • C and maintained at a relative humidity of 80% in the incubator.Leaf samples of 0, 3, 6, 12, and 24 h treated with low and high temperature stresses were collected for RNA extraction.For drought stress treatments, the tissue-cultured plants were grown in Hoagland nutrient solutions [38] containing 20% PEG6000, and incubated at different times (0 h, 3 h, 6 h, 12 h, 1 day, 3 day, 4 day, and 7 day).Then, the leaf samples of each drought-treated time point were collected for RNA extraction, and samples from untreated plants were used as controls.
For hormone and salt treatments, the tissue-cultured seedlings were treated with abscisic acid (ABA) (200 µmol/L), gibberellin (GA) (100 µmol/L), and high salt (1 M NaCl), respectively.Among these chemicals, ABA and GA were diluted in distilled water that contained 0.05% (v/v) ethanol.The diluted ABA, GA, and NaCl solutions were sprayed on the leaves and stems of seedlings until the runoff occurred.For control plants, the distilled water that contained 0.05% (v/v) ethanol was sprayed on the leaves and stems of seedlings.Leaf samples were harvested 0, 0.5, 2, 6, 12, 24, and 48 h after treatments.In all treatments, one leaf from each of the five plants was taken and mixed together for RNA extraction.

RNA Isolation and First-Strand cDNA Synthesis
Total RNA was isolated from the collected samples by the described methods [39], and the extracted RNA was digested with DNase (Promega, Madison, WI, USA) to remove genomic DNA contamination.The integrity and concentration of the RNA samples was detected by 1.5% agarose gel electrophoresis and NanoDrop 2000 (Thermo Scientific Inc., Waltham, MA, USA), respectively.Then, the RNA samples were reverse transcribed into First-strand cDNA with the RevertAid™ First Strand cDNA Synthesis Kit (TaKaRa, Shiga, Japan).

Identification and Isolation of Mad-box Genes in H. brasiliensis
Twenty full-length cDNA sequences of H. brasiliensis Mad-box genes were obtained from RNA sequencing.The cDNA sequences of these genes were compared with the Transcriptome Shotgun Assembly (TSA) and Expressed Sequence Tags (EST) of H. brasiliensis in the NCBI database (http: //www.ncbi.nlm.nih.gov/), or searched against the Hevea genome database.Then, the NCBI ORF Finder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html)and Softberry (http://linux1.softberry.com/)were used to determine open reading frames (ORFs) of candidate mRNA or genome DNA sequences.In addition, in order to confirm the presence of the Mad-box domain, all the candidate HbMad-box genes were further validated by conserved domain searching using CDD (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and InterProScan (http://www.ebi.ac.uk/interpro/scan.html).After the similarity comparison, ORF Finder, and conserved domain searching, redundant sequences were removed and the sequences of the HbMad-box gene were obtained.
The gene-specific primers used to amplify the corresponding full-length cDNA sequences of HbMad-box genes were designed by Primer 3.0 (http://primer3.ut.ee/).The primer pairs for all HbMad-box genes are listed in Table S1.RT-PCR amplification of HbMad-box genes was conducted using Pyrobest™ DNA polymerase (TaKaRa, Japan) according to the instructions.

Protein Properties and Gene Structure Analysis of HbMad-box Genes
We used the ProtParam tool (http://web.expasy.org/protparam/) to predict the theoretical molecular weight (Mw) and isoelectric point (PI) of HbMad-box gene proteins.To further analyze the structural diversity of HbMADS-box genes, the exon-intron structures of HbMad-box genes were identified by comparing the coding sequence with their corresponding genomic sequence using the FGENESH-C tool (http://linux1.softberry.com), as previously described [40].

Multiple Sequence Alignments and Phylogenetic Analysis
In the present study, amino acid sequence identities and multiple alignments of 26 HbMad-box proteins were calculated using DNAMAN6.0.To examine the evolutionary history and phylogenetic relationships of the HbMADS-box genes, we constructed a phylogenetic tree using MEGA6.0(http:// www.megasoftware.net/) based on multiple sequence alignment of HbMADS-box TFs from Arabidopsis, Oryza sativa L., Vitis vinifera L., Jatropha carcas L., and H. brasiliensis.

Quantitative Real-Time PCR (qRT-PCR) Analysis
The real-time quantitative RT-PCR (qRT-PCR) was performed according to the following procedures: 94 • C for 30 s, 94 • C for 5 s, followed by 40 cycles, 60 • C for 15 s, and 72 • C for 10 s.The reaction volume was 20 µL: including 60 ng cDNA per sample, 1×SYBR ® Premix Ex Taq™ (TaKaRa, Shiga, Japan), and 0.4 µM per primer.The reaction was carried out in 96-well plates using the CFX96™ Real-Time System (Bio-Rad, Hercules, CA, USA).The 18S rRNA gene (GenBank accession No.: AB268099) was used as a reference gene in the qRT-PCR reaction [41].After the reaction, we used Bio-Rad CFX Manager Software 3.0 (Bio-Rad, Hercules, CA, USA) to analyze and visualize the data as previously described [42].All primers used for qRT-PCR analysis were designed by Primer3.0 (http://frodo.wi.mit.edu/primer3).The primer sequences are given in Table S2.

Statistical Analysis
Data and graphical analysis was performed with Sigma Plot 12 software (Systat Software Inc., San Jose, CA, USA).The 2 −∆∆CT method was used to calculate the relative expression levels of all HbMad-box genes [43].The data are represented by the mean ± SD (standard deviation) of the three biological repeats.The statistical significance of the values was determined by a t-test.
Pairwise sequence comparisons were performed to check the sequence identities between HbMADS-box proteins.The results showed that the identities between two HbMADS-box genes ranged from 13.67% to 98.57% (Table S4).The average sequence homology between two HbMADS-box genes was 39.99%.The largest sequence identity was observed between HbMADS4 and HbSVP1 (98.5%).HbAGL30 and HbMADS3 showed the least sequence identity (13.67%).

Analysis of Conserved Domains and Structural Features of HbMad-box Proteins
Analysis of conserved domains showed that all of the deduced HbMADS-box proteins contained a highly conserved MADS-box domain of approximately 50 amino acid residues in length, a semi-conservative K domain with an obvious coiled-coil region, a less conservative I domain, and a non-conservative C domain (Figure 1), which indicates that they belonged to the MADS-box TF family.We further analyzed the conserved domain of HbMADS-box proteins, and the results indicated that for amino acids 3 (R), 17 (R), 20 (T), 23-24 (KR), 27 (G), 30-31 (KK), 34 (E), 38-39 (LC), 48 (F), and 52 (G), up to 13 sites were conserved in the MADS-box domains (Figure 1), indicating that the deduced amino acid sequences of HbMADS-box proteins shared significant homology with each other.The conserved motifs and structural characteristics analysis showed that all HbMADS-box genes included four conserved motifs: motifs 1 and 3 belonged to the MADS-box domains, whereas motifs 2 and 4 belonged to the K-box domains.Most HbMADS-box proteins (except HbAGL30 and HbAGL61) contained the MADS-box and K-box domains (Figure S2), which is similar to the previously reported MADS-box conserved motifs in other species [44].

Gene Structures and Sequence Characteristics of HbMad-box Genes
Sequence analysis revealed that all of the deduced H. brasiliensis MADS-box genes belong to the type II (MIKC-type) MADS-box genes.These HbMad-box genes usually contained multiple introns and exons, with a maximum of 11 exons (Figure S3).According to their predicted structures, the 26 candidate HbMADS-box genes could be divided into eight groups (Table S5).The first group of HbMADS-box genes consisted of eight exons interrupted by seven introns, including eight type II HbMADS-box genes (HbSEP, HbAGL9.2,HbCMB1, HbCMB1-L, HbSVP1, HbMADS4, HbAGL65, and HbMADS5).The second and third groups included five type II HbMADS-box genes, which contained six exons interrupted by five introns (HbAGL6, HbMADS2, HbMADS27, HbSVP2, and HbTT16), and seven exons interrupted by six introns (HbAP1, HbMADS3, HbAGL12, HbAG, and HbDEFL), respectively.The fourth and fifth groups consisted of two type II HbMADS-box genes, which contained 10 exons interrupted by nine introns (HbAGL9.1 and HbAGL8), and five exons interrupted by four introns (HbAGL11 and HbPMADS2), respectively.The sixth group included two type II HbMADS-box genes (HbAGL30 and HbMADS1), which contained 11 exons interrupted by 10 introns.The seventh and eighth groups each contained one type II HbMADS-box gene (HbAGL15 and HbAGL61), which consisted of four exons interrupted by three introns, and only one exon, respectively.These results revealed that most of the HbMADS-box genes (69.2% belong to the first, second, and third groups) within the same group shared conserved exon-intron structures, and the other HbMADS-box genes (30.8% of the fourth to eighth groups) with the different exon-intron structures may belong to distinct types of HbMADS-box genes in the rubber tree (Table S5).These characteristics are consistent with the features of MADS-box genes in other flowering plants, such as Arabidopsis, Oryza sativa, and Vitis vinifera [45][46][47].

Phylogenetic Analysis of HbMad-box TFs
Based on the phylogenetic tree analysis, the 26 HbMADS-box genes could be divided into three clades (I, II, and III), which contained 11 main branches: SEP, AGL6, FUL/AP1, AGL12, AG, SOC1, AGL15, AGL17, SVP, BS, and AP3/PI subfamilies (Figure 2).Clade I contained 12 type II HbMADS-box genes, and it could be further categorized into two subclades, which contained nine and three members.Clade II contained six type II HbMADS-box genes, and it could also be divided into two subclades, and one subclade contained one member, whereas the other contained five members.Clade III contained eight type II MADS-box genes, which had two subclades, one of which contained three members and the other contained five members.A phylogenetic tree of MADS-box genes was generated by the neighbor-joining (NJ) method.The multiple sequence alignment and construction of phylogenetic tree were performed with MEGA6.06 using the neighbor joining method with 1000 bootstrap replicates.The proteins were clustered and divided into three distinct sub-families.The three sub-families were further divided into 12 sub-groups.

Expression Profiles of HbMad-box Genes in Various Developmental Stages during Rubber Tree Reproductive Development
In the present study, analyzing the expression patterns of HbMADS-box genes in different tissues and organs (including roots, stem, stem tips, leaves, labeled bark, xylem, latex, fruits, inflorescence, male flowers, and female flowers) of H. brasiliensis by qRT-PCR, the results showed that most of the deduced HbMADS-box genes might be primarily involved in floral organ differentiation and inflorescence development.For example, HbSEP, HbAGL9.1,HbAGL9.2,HbCMB1, HbCMB1-L, HbAGL6, HbAG, HbDEFL, and HbPADS2 had significantly higher expression in stem tips and floral organs than other tissues; HbAP1 showed the highest expression pattern in stem tips compared to other tissues (Figure 3); HbAGL8 was significantly higher expressed in stem tips and leaves than other tissues (Figure 3); the MADS-box genes HbAGL15, HbTT16, HbMADS2, and HbAG were higher expressed in fruits than other tissues (Figure 3); and the MADS-box genes HbAGL9.2 and HbAGL6 were specifically expressed in floral organs and they were higher expressed in male flowers than other tissues (Figure 3).The higher expression of the HbMad-box genes in floral organs and shoot tips suggests that they may play a specific role in the corresponding tissues.However, in contrast, some MADS-box genes displayed different tissue expression patterns, such as HbAGL15, HbAGL30,

Expression Patterns of HbMADS-box Genes in Response to Abiotic Stress
Analyzing the expression patterns of 12 floral-enriched HbMADS-box genes under cold, high-temperature, drought, and salt stresses by qRT-PCR, the results indicated that the expressions of 12 HbMADS-box genes specifically related to the floral organ were affected by low-temperature stress (Figure S4).Among the 12 cold-regulated HbMADS-box genes related to floral organ development, most of them (HbSEP, HbAGL9.1,HbAGL9.2,HbCMB1, HbCMB1-L, HbAGL6, HbAG, HbDEFL, and HbMADS2) were rapidly induced by cold stress, and reached the highest levels at 3 or 6 h after treatment, and then declined at 12 or 24 h after treatment.Only three HbMADS-box genes, HbAGL8, HbAP1, and HbTT16, exhibited down-regulated expression at all analyzed time-points in response to cold stress.
Under polyethylene glycol induced drought stress, the expression patterns of 12 HbMADS-box genes related to floral organ development were markedly affected (Figure S6).

Expression Patterns of HbMADS-box Genes in Response to Phytohormones
After ABA treatment, HbAGL9.1,HbAGL6, HbAP1, and HbPADS2 displayed obvious up-regulation across all time points (Figure S8).Their expression levels exhibited rapid up-regulation under ABA treatment, and peaked at 48, 48, 12, and 48 h after treatment, respectively.HbAGL9.2 and HbCMB1-L showed similar expression patterns; their transcripts exhibited down-regulated expression at the 6-h time point, but they were significantly up-regulated at other time points.HbSEP, HbCMB1, HbAGL8, HbAG, and HbDEFL were down-regulated at the time points of 24, 12, 0.5, 24, and 24 h compared with the control, respectively, but they displayed obvious up-regulation at other time points.In comparison, HbTT16 showed an irregular expression pattern under ABA treatment; it was up-regulated at 48 h, but was down-regulated at 0.5, 2, 12, and 24 h.

Discussion
Understanding the flowering regulation of H. brasiliensis is important for accelerating the breeding process of this species.Based on studies of Arabidopsis and other plants, extensive efforts have been devoted to clarifying the molecular mechanism of reproductive development in plants [48,49], and techniques such as RT-qPCR and high-throughput RNA sequencing have been used in many transcriptional-level studies to identify the genes that regulate the metabolism of plant reproductive development.However, available information about TFs related to reproductive development in H. brasiliensis is still limited.To date, no study has been performed to determine the expression patterns of HbMad-box genes relative to rubber tree reproductive development.
As an important gene family, Mad-box genes are widespread in the plant kingdom.The number of Mad-box genes varies considerably in genomes of different species, from 20 in Physcomitrella patens (Hedw.)Bruch & Schimp.to 167 in Brassica rapa L. [50,51].In this study, we performed genome-wide identification and analysis of the MADS-box gene family related to flower development in H. brasiliensis.A total of 20 MADS-box genes were newly identified in the H. brasiliensis genome, which was lower than the number of MADS-box genes in the woody plant Malus pumila Mill.[52].A domain search using EMBL with the corresponding H. brasiliensis candidate protein sequences confirmed that 26 of the sequences contained a 'MADS' domain.We classified all 26 putative H. brasiliensis MADS-box proteins into three clades, which were consistent with the classifications of MADS-box gene family members in other flowering plants [53].Sequence alignment analysis of 26 MADS-box genes revealed that their ORFs ranged from 609 to 1116 bp, and predicted protein lengths that ranged from 202 to 371 amino acids (Table 1).Subsequent gene structures analysis showed that most of the HbMad-box genes usually contained multiple introns, with a maximum of 11 introns; the exception was HbAGL61, which did not have any introns.Interestingly, we observed that only MIKC-type (no M-type) MADS-box genes existed in H. brasiliensis.In contrast, the number of M-type MADS-box genes in Arabidopsis is more than that of the MIKC C -type [54].Furthermore, we found that the expansion of MIKCC-type and MIKC*-type MADS-box genes was disproportionate, and there were more MIKC C -type (twenty-one members) than MIKC*-type (five members) MADS-box genes presented in H. brasiliensis, which might be related to the fact that MIKC C -type genes conducted as functional genes to perform more complex functions in H. brasiliensis flower organogenesis.
Most studies on the functional identification of MADS-box genes were conducted on model plants [55][56][57], which also substantially contributed to revealing the diverse functions of MADS-box genes.In order to detect the expression of HbMADS-box genes in different tissues, we analyzed the expression patterns of HbMADS-box genes in 11 sample tissues and four stages of inflorescence development by qRT-PCR.The results revealed that for close orthologs of SEP and AGL6 in the rubber tree, HbSEP, HbAGL9.1,HbAGL9.2,HbCMB1, HbCMB1-L, and HbAGL6 showed higher expression during stem tip and floral organ development; for close orthologs of FUL-AP1 and AG in the rubber tree, HbAGL8, HbAP1, and HbAG, also exhibited higher expression during stem tip and floral organ development (Figure 3).The same expression pattern was also observed in subclade AP3-PI; HbDEFL and HbPADS2 were found to be highly expressed during stem tip and floral organ development, which indicates association with the development of reproductive organs.These findings are consistent with those of previous results that revealed the MIKC-type MADS-box gene as the flower homeotic gene that plays a dominant role in floral organ development [58].
Interestingly, HbAGL6 was progressively increased in four stages of inflorescence development, and the expression level of HbAGL6 in the male flowers was almost four-fold more than the female flowers in H. brasiliensis, which indicates that HbAGL6 may be involved in the development of inflorescence and floral organs, and especially participates in the development of male flowers organs (Figure 3).In addition, HbDEFL also displayed strong expression in inflorescence and was almost four-fold more abundant in the male flowers than the female flowers in H. brasiliensis, which indicates that HbDEFL may play an important role in the development of male flowers.In contrast, HbTT16 showed strong expression in the fruits and female flowers, and the expression level of HbTT16 in the female flowers was almost 40-fold more than the male flowers in H. brasiliensis, which indicates association with the development of female organs in the rubber tree.In other plant species, TT16 homologs were also demonstrated to mediate the crosstalk between the endothelium and nucellus for the development of female organs [59,60].
Besides being involved in the essential regulation of the development of floral organs, many studies have shown that MADS-box genes also play an essential role in the response to various stresses and exhibit differential expression patterns under abiotic stress [61,62].Rubber trees originated in the Amazon Basin of South America, but are now widely planted in the northern margin of tropical areas, such as Southeast Asia countries (e.g., Thailand, Vietnam, and southern China), which are often suffer from cold, drought, typhoon, and other abiotic stresses.
To understand the responses of HbMADS-box genes to various stresses, the expression profiles of 12 floral-enriched HbMADS-box genes were examined in leaves of Hevea seedlings after hormone, salt, cold, high-temperature, and drought stress treatments.In the present study, almost all of the HbMADS-box genes specifically related to floral organs exhibited changes in expression patterns when responding to hormone, salt, cold, high-temperature, and drought stresses (Figures S4-S9).Our results show that different HbMADS-box genes responded differently to abiotic stress, which indicates that stress responses of H. brasiliensis are regulated by many factors, involving a variety of signaling pathways.In addition, the hormone and salt stress responses of 12 floral-enriched HbMADS-box genes were more intense than those of temperature and drought stress responses.Among them, HbAGL9.2,HbAGL6, and HbAP1 were strongly induced by ABA; HbAGL9.1,HbAGL9.2,HbCMB1-L, HbAGL6, HbAGL8, HbAP1, HbTT16, and HbPADS2 were strongly induced by GA3; and HbAGL9.2,HbCMB1-L, HbAGL6, HbAGL8, HbAP1, HbTT16, and HbPADS2 were strongly induced by salt treatments, whereas their expressions were slightly affected by cold, high-temperature, and drought stresses.The expressions of HbAGL9.2,HbAGL6, and HbAP1 were strongly induced by a variety of hormones, which revealed that these genes might play a crucial role in hormone signaling in H. brasiliensis.
Overall, we characterized HbMADS-box genes as multifunctional regulators in H. brasiliensis based on tissue specific expression analysis, as well as various stress responses.Our results revealed that HbSEP, HbAGL9.1,HbAGL9.2,HbCMB1, HbCMB1-L, HbAGL6, HbAGL8, HbAP1, HbAG, HbDEFL, and HbPADS2 may mainly regulate the differentiation of flower buds and could help regulate reproduction in H. brasiliensis; whereas HbTT16 may mainly regulate the development of fruits and female organs in the rubber tree.Expression profiling revealed that different HbMADS-box genes responded differently to abiotic stresses, which indicates that abiotic stresses in H. brasiliensis are regulated by many factors and various signaling pathways.Our findings will promote the development of technology that can control the reproduction of H. brasiliensis by genetic engineering.In the near future, we will be able to verify the role of these HbMADS-box genes in the transition of vegetative to reproductive development using transgenic H. brasiliensis, as well as transgenic Arabidopsis.

Conclusions
In this study, 20 new HbMad-box genes were identified in the rubber tree.Subsequently, the bioinformatics characteristics, sequence identity, conserved domains, gene structure, and phylogenetic relationship of these genes were systematically analyzed.Expression profiling revealed that HbMad-box genes were differentially expressed in various tissues, which indicated that HbMad-box genes may exert different functions throughout the life cycle.Additionally, HbSEP, HbAGL9.1,HbAGL9.2,HbCMB1, HbCMB1-L, HbAGL6, HbAGL8, HbAP1, HbAG, HbDEFL, HbTT16, and HbPADS2 were found to be associated with the differentiation of flower buds and may be involved in flower development.All of these floral-enriched HbMADS-box genes were regulated by hormone, salt, cold, high-temperature, and drought stresses, which revealed that abiotic stresses in

Figure 1 .
Figure 1.Phylogenetic analysis of HbMADS-box genes related proteins.MADS-box, I domain, K domain, and C domain are marked, respectively; coiled-coil domain is boxed.

Figure 2 .
Figure 2. Phylogenetic analysis of HbMADS-box genes with other MADS-box genes by MEGA version 6.0.A phylogenetic tree of MADS-box genes was generated by the neighbor-joining (NJ) method.The multiple sequence alignment and construction of phylogenetic tree were performed with MEGA6.06 using the neighbor joining method with 1000 bootstrap replicates.The proteins were clustered and divided into three distinct sub-families.The three sub-families were further divided into 12 sub-groups.

Figure 3 .
Figure 3. Relative expression levels of HbMad-box genes were determined by qRT-PCR and normalized by 18S rRNA gene expression.Relative expression levels of HbMad-box genes were determined by qRT-PCR and normalized by the 18S rRNA gene expression.For each gene, the transcript level in the root was used to normalize the transcript levels in other tissues.Values are means ± SD (standard deviation) of three biological replicates.1-14 represent roots; stem; stem tips; leaves; labeled bark; xylem; latex; fruits; 3 cm inflorescence; 6 cm inflorescence; 9 cm inflorescence; 12 cm inflorescence; male flowers; and female flowers, respectively.