HaMADS3, HaMADS7, and HaMADS8 are involved in petal prolongation and floret symmetry establishment in sunflower (Helianthus annuus L.)

Plant materials

The sunflower GuangWu Mountain Wild Type (WT) was treated by space mutation via “Shenzhou No.4” (Wu et al., 2020). Space mutation is an efficient technology used for changing the genetic characteristics of living organisms through pace-related factors, such as space radiation, space microgravity, cosmic irradiation, and space magnetic fields. The lpm plants were screened from the offspring of space mutants after several years of inbreeding. These lpm plants were divided into five types: type V, IV, III, II, and I. From type V to I, the mutant traits gradually became obvious. In type V plants, the phenotype was similar to that of WT plants. The type I plants possessed the extreme phenotype, with prolongated petals and degenerated stamens in the most disc floret (He et al., 2022). In our experiment, only the type I lpm plants were used for analyzing the effect of MADS-box genes on petal prolongation and floret symmetry establishment. The experimental materials were planted in the field of the College of Life Sciences, China West Normal University (Nanchong, Sichuan Province, P. R. China, 30°49′N, 106°4′E). All plants were sown in March and collected in July and August 2022 under natural conditions (an average temperature of 25 ± 5 °C and average rainfall 85.8 mm). Fertilization and weeding were performed once a month and the plants were irrigated twice a month.

Morphological and histological analysis

When blooming, the floret materials from WT and lpm, including ray and disc florets, were harvested to compare their morphogenesis, and the floret petals of WT and lpm plants (within the 1^st (ray floret), 5^th, 15^th, and 19^th parastichy) were collected to analyze the petal length using a vernier caliper with 20 biological repeats. The petal length was determined by measuring the petal naturally straightened on a flat surface from bottom to top. Subsequently, the disc florets at different floral developmental stages (25 days prior to blooming (DPB), 15, 5, and 0 DPB) were harvested to analyze the morphological differences of disc florets, stamens, and pistils between WT and lpm via stereomicroscope (Olympus, Tokyo, Japan). The disc florets (at 30, 25, 20, and 15 DPB) were collected and fixed in FAA solution (70% ethanol: acetic acid: formaldehyde, 90:5:5, v/v) for histological analysis. Subsequently, all samples were dehydrated in graded ethanol (30%, 50%, 70%, 80%, 90%, 100%, and 100% ethanol), each for 1 h at room temperature. The dehydrated materials were then incubated in the ethanol with xylene as follows: ethanol/xylene (1:3, v/v) for 2 h, ethanol/xylene (2:3, v/v) for 2 h, ethanol/xylene (1:1, v/v) for 2 h, three changes of xylene for 1 h. The materials were then set in xylene/paraffin (1:1, v/v) for 12–20 h at 40 °C and in paraffin at 60 °C with three changes for 1 h as interval, and then embedded in pure molten wax. Samples were sectioned into pieces with 7 μm using microtome (Leica, Wetzlar, Germany). The sections were deparaffinized 10 min in xylene two times and rehydrated as follows: 100%, 100%, 95%, 80%, and 70% ethanol and distilled water, each for 5 min. Sections were then stained with toluidine Blue O (G3668, Solarbio, Beijing, China) for 5 min and washed in water. After that, sections were dehydrated for 1 min through the ethanol series (70%, 80%, 95%, 100%, and 100% ethanol) and set in xylene for 15 min. Finally, the sections were covered with glass cover and observed under optical microscope (Olympus, Tokyo, Japan). Petals from different positions (the 1^st, 5^th, 15^th, and 19^th parastichy) on the capitulum were selected to analyze the cells located in the upper, middle, and lower part of the petal in order to record cell size and number data.

RNA-sequencing and data analysis

The total RNA of root, stem, leaf, and flower were extracted using an E.Z.N.A®Plant RNA Kit (R6827; Omega, Norcross, GA, USA) and sent to Shanghai Majorbio Biopharm Technology Co., Ltd. (Shanghai, China) for RNA-sequencing. The cDNA library was constructed using the Illumina NovaSeq Reagent Kit (Illumina, San Diego, CA, USA) and 1 μg of RNA, and RNA-sequencing was performed through the Illumina NovaSeq 6000 platform. The raw paired end reads were trimmed, and the quality of the raw paired end reads was controlled by fastp (Chen et al., 2018). The clean reads were aligned to the sunflower genome (https://www.ncbi.nlm.nih.gov/genome/?term=txid4232[orgn]) using HISAT2 software (Kim, Langmead & Salzberg, 2015), and the mapped reads of each sample were assembled via StringTie (Kovaka et al., 2019). The number of raw and clean reads, percentage of total-mapped reads, and percentage of unique-mapped reads obtained from each sample are listed in Table S1. The expression levels of the transcripts were calculated using the transcripts per million reads (TPM) method. The genes were annotated through web sites, including GenBank of National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/genbank/) and the European Molecular Biology Laboratory (EMBL) (https://www.embl.de/).

The identification of MADS-box genes and analysis of conserved motifs and phylogenetic tree

The putative MADS-box genes were identified based on the sunflower genome (https://www.ncbi.nlm.nih.gov/genome/?term=txid4232[orgn]) and our RNA-seq data from a collection of root, stem, leaf, and flower tissues. Using all MADS-box protein sequences of A. thaliana obtained from the TAIR database (http://www.arabidopsis.org) as queries, the MADS-box sequences of sunflower were identified through BLASTP. Subsequently, based on the annotations in our RNA-seq data, the members of the MADS-box family were finally identified and analyzed in this study. Notably, the MADS-box members that could be effectively detected in our RNA-seq data were identified, but not all the members of the MADS-box family in sunflower. The detailed information of these MADS-box genes is provided in Tables S2–S4, including transcript sequences, detailed annotations, and TPMs.

The protein sequences of the 11 reported MADS-box genes in sunflower (Dezar et al., 2003; Shulga et al., 2008) were then aligned with all candidate MADS-box genes in our study using DNAMAN 8.0 (https://dnaman.software.informer.com/) and annotated in Table S3. Subsequently, the conserved motifs of protein sequences encoded by the candidate genes were analyzed and visualized via MEME (http://memesuite.org/tools/meme). The multiple sequence alignment was implemented by Clustal 2.0 (http://www.clustal.org/) using the protein sequences of the MADS-box genes of sunflower and A. thaliana (see NCBI IDs and TAIR IDs in Table S5). A phylogenetic tree was constructed using MEGA 11.0.13 (Tamura, Stecher & Kumar, 2021), using neighbor-joining (NJ) with 1,000 bootstrap replications, p-distance model, and the pairwise deletion of gaps.

RNA extraction and qPCR analysis

The roots, stems, leaves, and flowers from WT were collected for analyzing the tissue expression pattern of the MADS-box genes. The disc florets of WT and lpm at different stages were collected to explore the temporal expression pattern of MADS-box genes. The petals, pistils, and bracts in the WT disc florets and abnormal disc (ray-like) florets of lpm were used to investigate the spatial expression patterns of these genes. The petals from ray and disc florets in WT were applied to explore the function of candidate MADS-box genes during the process of petal development, and the floret petals of WT and lpm plants within the 1^st, 5^th, 15^th, and 19^th parastichy were collected to analyze the expression levels of MADS-box genes in petals across different positions. The detailed information of materials used for quantitative real-time PCR (qPCR) are listed in Table S6. All materials used in this study for qPCR were gathered from 15 individual plants that were grown in the same conditions and then mixed for RNA extraction. At least 0.5 g of materials were harvested, immediately frozen in liquid nitrogen, and then stored at −80 °C for no more than 3 months.

Total RNA was extracted using the E.Z.N.A^®Plant RNA Kit. NanoDrop2000 (ThermoFisher, New York, NY, USA) was used to detect quantification, OD260/280, OD260/230, and contamination assessment of RNA. Then 1% agarose gel electrophoresis was used to analyze RNA integrity via electrophoresis apparatus (Bio-Rad, Hercules, CA, USA) and gel imaging system (Bio-Rad, Hercules, CA, USA), using Tris-Borate-EDTA (BL540A, Biosharp, Hefei, Anhui, China) as buffer and SuperRed (BS354A, Biosharp, Hefei, Anhui, China) as dye. See detailed results in Table S7 and Fig. S1.

The cDNA was obtained using the PrimeScript^™ RT reagent Kit with gDNA Eraser (RR047A, TaKaRa, Shiga, Japan). The reaction for removing gDNA was performed using 10 μL mix including: 2.0 μL 5 × gDNA eraser buffer, 1.0 μL gDNA eraser, 1.0 μg total RNA, and adding RNase free dH₂O up to 10 μL. The 10 μL mix was then incubated at 42 °C for 10 min. RNA reverse transcription was performed using 20 μL mix including: 10 μL mix after removing gDNA reaction, 1.0 μL PrimeScript RT enzyme mix I, 1.0 μL RT primer mix, 4.0 μL 5×PrimeScript buffer 2, and 4.0 μL RNase free dH₂O. The 20 μL mix was incubated at 37 °C for 15 min and 85 °C for 5 s. Finally, the 20 μL cDNA solution was diluted 10 times and stored at −20 °C.

qPCR was performed using TB Green^® Premix Ex Taq^TM II (RR820A, TaKaRa, Shiga, Japan) via CFX96^TM Real-Time System (Bio-Rad, Hercules, CA, USA). qPCR reaction was performed using 20 μL mix including: 10 μL TB Green Premix Ex Taq II (Tli RnaseH Plus), 1.0 μL forward primer (10 μM), 1.0 μL reverse primer (10 μM), 7.0 μL Rnase free dH₂O, and 1.0 μL cDNA. qPCR was performed with following program: 30 s at 95 °C, then 40 cycles of 5 s at 95 °C, and 30 s at 56 °C. Data analysis for qPCR was performed by Bio-Rad CFX maestro 1.1 software (version 4.1.2433.1219, https://www.bio-rad.com/). For each sample, 15 biological (pooled) and three technical replicates were conducted. Gene expression levels were calculated using the 2^−∆∆Ct method, with elongation factor 1-alpha (EF-1α; GenBank accession: XM_022117715.2; also known as EF1A) as the reference gene (De Paolis et al., 2022; He et al., 2022; Najafabadi & Amirbakhtiar, 2023). All primers are listed in Table S8. The MIQE checklist for qPCR is listed in Table S9.

Tissue expression pattern and correlation analysis

The tissue expression pattern was analyzed based on the TPMs in our RNA-seq data (Table S4) and visualized via the pheatmap package of R 4.3.1 (R Core Team, 2023). The relationship between petal length and expression levels of HaMADS3, HaMADS7, and HaMADS8 were analyzed using the Pearson correlation coefficient (PCC). Pearson correlation analysis was performed via R WGCNA package. PCC > 0.6 and P-value < 0.05, or PCC < −0.6 and P-value < 0.05 were considered positively or negatively correlated, respectively (Zhang et al., 2015).

Analysis of cis-acting elements

The 2,000 bp sequences (−2,000 bp) upstream of the transcription start site of HaCYC2c, HaMADS3, HaMADS7, and HaMADS8 were obtained from NCBI (https://www.ncbi.nlm.nih.gov/). The cis-acting elements in the promoter region of HaCYC2c bound by MADS-box transcriptional factors, as well as in the promoter region of HaMADS3, HaMADS7, and HaMADS8 bound by HaCYC2c, were analyzed by JASPAR (Khan et al., 2018). The results were visualized using TBtools (Chen et al., 2020).

Statistical analysis

Statistical analysis was performed by one-way ANOVA and Tukey multiple comparison using SPSS 22.0 (IBM, Armonk, NY, USA).

Morphogenesis of florets in WT and lpm

To analyze the morphogenetic difference between WT and lpm, the capitulum structure and composition, floret symmetry, and morphogenesis of floral organs were investigated. The results showed that the capitulum of WT was composed of multiple whorls of actinomorphic disc florets with shorter petals surrounded by a single whorl of zygomorphic ray florets with longer petals (Fig. 1A). In WT plants, the zygomorphic ray floret consists of zygomorphic petals and a pseudo-ovary, while the actinomorphic disc floret possesses a single bract, a pair of pappus, one actinomorphic corolla fused by five petals, five synantherous stamens, and one pistil (Figs. 1B and 1C). In lpm, an insignificant difference was observed in ray florets, compared with WT (Fig. 1C). Comparatively, in lpm, the petals of abnormal disc florets were significantly elongated and the stamens were partially degenerated, while the structure of other floral organs remained normal. The symmetry of disc florets in lpm was transformed from actinomorphic to zygomorphic, their outline was similar to that of ray florets in WT, and these florets were thus defined as ray-like florets (Figs. 1A–1C). Notably, the ray-like florets in lpm were not completely uniform in morphology, and the transformation of petal length and symmetry gradually weakened from outside to center (Fig. 1D).

Analyzing the development schedule of floral organs of disc florets via the morphological and histological method indicated that the difference between WT and lpm at 25 DPB was insignificant, while the divergency initially emerged at 20 DPB. In particular, the petals in lpm started to elongate faster than those of the WT plants, and the discrepancy of petal length gradually increased afterwards (Fig. 1E). Coincidentally, until 25 DPB, stamen development was almost synchronized between lpm and WT, indicating that they shared similar stamen morphology. Subsequently, stamens developed following the normal schedule in WT, while stamen development in lpm was suppressed, causing only stamen traces to be present in the disc florets when blooming (Figs. 1E and 1F). Regarding pistils, an insignificant difference was observed between WT and lpm except in color. The pistil stigma was light yellow in disc florets of WT, then gradually darkened at 10 DPB, and finally transformed into purple-red, but the stigma color remained light yellow in lpm (Fig. 1E).

To determine the cause of petal prolongation in lpm, the cell size and number were analyzed in petals. The results showed an insignificant difference in the cell size presented between WT and lpm (Fig. 1G). The cell number of petals in disc florets of lpm was significantly higher than that of WT, despite no obvious difference in ray floret petals between WT and lpm (Fig. 1G). These results revealed that the increase in cell number was the main contribution to petal prolongation of disc florets in lpm.

Identification of MADS-box genes in sunflower

Based on our RNA-seq data and sunflower genome (https://www.ncbi.nlm.nih.gov/genome/?term=txid4232[orgn]), 29 MADS-box candidate genes were identified: 11 previously reported genes and 18 unreported genes (Dezar et al., 2003; Shulga et al., 2008). They were labeled HaMADS1 to HaMADS29 (see Table S3 for details).

In order to classify the obtained MADS-box genes, 13 conserved motifs of the proteins encoded by these genes were further analyzed (Table 1). The results showed that among these genes, HaMADS28 and HaMADS29 belonged to the Type I class because their proteins only contained motifs 1 and 3, while the remaining 27 genes belonged to the Type II class because their proteins contained motifs 1, 2, 3, and 5 simultaneously (Fig. 2A). The results of the phylogenetic analysis indicated that these 29 MADS-box genes were distributed in 11 clades, and the detailed information is listed in Fig. 2B.

Identifying the MADS-box genes with high expression in flower

To predict the function of the 29 MADS-box genes, their tissue expression patterns were analyzed based on our RNA-seq data (Table S4). Nineteen genes presented their maximum expression levels only in flower: HaMADS1, HaMADS3, HaMADS7, HaMADS8, HaMADS9, HaMADS11, HaMADS12, HaMADS15, HaMADS16, HaMADS18, HaMADS19, HaMADS20, HaMADS21, HaMADS22, HaMADS23, HaMADS24, HaMADS25, HaMADS26, and HaMADS27 (Fig. 3A). The other 10 genes that had high expression levels both in flower and other organs, or were hardly expressed in flower, were set aside for further investigation (Fig. 3A). Therefore, only 19 MADS-box genes were considered to play important roles in flower development.

Differential expression of MADS-box genes between WT and lpm during flower development

To analyze the precise action period of the 19 MADS-box genes during flower development, their expressions were analyzed at four developmental stages in WT and lpm using qPCR. In WT sunflower, 11 genes exhibited higher expression at 25 DPB than other stages: HaMADS7, HaMADS8, HaMADS9, HaMADS11, HaMADS12, HaMADS15, HaMADS21, HaMADS23, HaMADS24, HaMADS26, and HaMADS27. Among these genes, seven (HaMADS9, HaMADS11, HaMADS15, HaMADS23, HaMADS24, HaMADS26, and HaMADS27) were hardly expressed or only presented with low levels in the other stages (Fig. 3B). The other eight genes (HaMADS1, HaMADS3, HaMADS16, HaMADS18, HaMADS19, HaMADS20, HaMADS22, and HaMADS25) showed higher levels at 5 DPB (Fig. 3B).

Comparing the expression levels of these genes in WT and lpm during floral development, the results indicated that six genes showed a reduced expression trend in lpm compared with WT: HaMADS1, HaMADS3, HaMADS9, HaMADS11, HaMADS22, and HaMADS15 (Fig. 3B). Eight genes were down-regulated in the early stages but up-regulated later, including HaMADS7, HaMADS12, HaMADS18, HaMADS19, HaMADS21, HaMADS25, HaMADS26, and HaMADS27 (Fig. 3B). Only one gene, HaMADS18, was up-regulation at all stages in lpm. In addition, four genes, HaMADS7, HaMADS8, HaMADS20, and HaMADS21, showed a steadily high level at all stages in lpm (Fig. 3B).

Notably, among these genes, eighteen were found strongly inhibited in lpm at the stage when they presented with peak value in WT. Eleven of these genes were strongly inhibited in lpm at 25 DPB: HaMADS7, HaMADS8, HaMADS9, HaMADS11, HaMADS12, HaMADS15, HaMADS21, HaMADS23, HaMADS24, HaMADS26, and HaMADS27 (Fig. 3B). Similarly, the other seven genes were inhibited in lpm at 5 DPB: HaMADS1, HaMADS3, HaMADS16, HaMADS19, HaMADS20, HaMADS22, and HaMADS25 (Fig. 3B). Perhaps these genes are related to the development of abnormal flower organs, such as petal and stamen, in lpm.

Identified MADS-box genes associated with petal development

To further identify the MADS-box genes related to the formation of abnormal floral organs in lpm, the difference of their expression in floral organs of disc floret were analyzed between WT and lpm. The results showed that 12 genes were involved in the development of floral organs in ray-like floret of lpm (HaMADS3, HaMADS7, HaMADS8, HaMADS11, HaMADS15, HaMADS16, HaMADS18, HaMADS19, HaMADS20, HaMADS21, HaMADS24, and HaMADS26) due to their regular expression discrepancy in floral organs between WT and lpm (Figs. 4 and S2). Among these genes, only five displayed expression variation in petals between WT and lpm: HaMADS21, HaMADS3, and HaMADS11, which were suppressed seriously in lpm petals from 5 DPB, and HaMADS7 and HaMADS8, which were significantly up-regulated (Fig. 4). And these five genes were considered to be associated with the formation of abnormal petal in disc floret of lpm.

Taken together, we found that HaMADS3 was strongly inhibited not only in various organs, including bract, petal, and pistil, but also in the whole process of flower development in lpm. HaMADS7 was significantly inhibited in the early stage (25 DPB) and became gradually up-regulated as the flower developed, when compared with WT (Figs. 3B, S2 and S3). However, comparatively, it was strongly up-regulated in all floral organs of lpm, particularly in petals (Figs. 4, S2 and S3). The expression of HaMADS8 exhibited strong stage-specificity in WT, with the most optimal value at 25 DPB followed by 5 DPB, but stage-specificity was weakened in lpm (Figs. 3B, S2 and S3). The expression of HaMADS8 was significantly higher in lpm than in WT, although no obvious difference was observed among the bract, petal, and pistil in WT, nor in lpm (Figs. 4, S2 and S3). HaMADS11 was strongly inhibited from 25 DPB to 5 DPB in lpm, and then slightly up-regulated compared with WT (Figs. 3B, S2 and S3). Its expression significantly decreased from 5 DPB to blooming only in the petal and pistil, but not in the bract of lpm plants (Figs. 4, S2 and S3). The expression of HaMADS21 was repressed originally at 25 DPB, but up-regulated from 15 DPB to blooming in lpm when compared with WT (Figs. 3B, S2 and S3). Its expression was inhibited in all floral organs except for bract in lpm (Figs. 4, S2 and S3). As the spatiotemporal expression results of HaMADS21 were contradictory, it was excluded for its contribution to the formation of abnormal floral organs of lpm. Therefore, HaMADS3, HaMADS7, HaMADS8, and HaMADS11 were considered to play important roles during the development of floral organs, such as petals.

In order to further identify the effects of HaMADS3, HaMADS7, HaMADS8, and HaMADS11 on the elongation and symmetry determinacy of petals, their expression levels were detected at various stages in ray and disc floret petals of WT. The results showed that HaMADS3 was strongly inhibited in the ray floret petals of WT (Fig. 5), similar to the ray-like floret petals in lpm (Fig. 4). HaMADS7 and HaMADS8 were significantly up-regulated in ray floret petals compared to disc floret petals in WT and showed a similar up-regulation trend in the ray-like floret petals of lpm (Figs. 4 and 5). It could be deduced that the lower expression of HaMADS3 and higher expression of HaMADS7 and HaMADS8 coordinately contributed to the formation of long and zygomorphic petals. On the other hand, the other gene HaMADS11 presented with up-regulation in ray floret petals when compared with disc florets in WT (Fig. 5), but was inhibited in ray-like floret petals in lpm (Fig. 4). For the logically conflicting exhibition of HaMADS11 in WT and lpm, only HaMADS3, HaMADS7, and HaMADS8 were considered to close association with the formation of long and zygomorphic petals in lpm.

Subsequently, the expression levels of HaMADS3, HaMADS7, and HaMADS8 were investigated in floret petals within the 1^st (ray floret), 5^th, 15^th, and 19^th parastichy in WT and lpm plants, and the correlation coefficients between the expression levels of the three MADS-box genes and petal length were analyzed. The results showed that the expression levels of HaMADS7 and HaMADS8 were significantly higher in the petals of ray-like florets from lpm and ray florets from WT compared to disc florets from WT. The opposite expression of HaMADS3 was detected (Fig. 6). Further correlation analysis suggested that the correlation coefficient between petal length and the expression levels of HaMADS3, HaMADS7, and HaMADS8 was −0.99, 0.96, and 0.97, respectively, suggesting a positive correlation between petal length and expression level of HaMADS7 and HaMADS8, and a negative correlation between petal length and expression level of HaMADS3 (Table 2). These results indicated the involvement of HaMADS3, HaMADS7, and HaMADS8 in petal prolongation and floral symmetry establishment in sunflower.

Interacting analysis between MADS-box genes and HaCYC2c

A previous study (He et al., 2022) found that the abnormally higher expression of HaCYC2c (LOC110936862) could result in the symmetrical transformation of disc floret from actinomorphic to zygomorphic in lpm sunflower, via recognizing the cis-elements in the promoter region of HaNDUA2. To further understand the interaction between these three MADS-box genes and HaCYC2c in regulating long and zygomorphic petal formation, their interaction relationship was analyzed. The results indicated that 14 cis-acting elements existed in the promoter region (−2,000 bp) of HaCYC2c as the potential transcriptional recognition binding sites for MADS-box transcription factors (Fig. 7 and Table 3). However, no HaCYC2c binding cis-acting element was found within the promoter regions (−2,000 bp) of HaMADS3, HaMADS7, or HaMADS8 (Fig. 7 and Table 3). These results implied that HaMADS3, HaMADS7, and HaMADS8 may regulate the expression activity of HaCYC2c, resulting in petal prolongation and symmetry transformation from actinomorphic to zygomorphic in lpm sunflower.