In this study, we selected fourteen species based on the degree of parasitism, comprising seven orders (Supplementary Table 1), including holoparasitic B. subcupularis, B. fungosa var. globosa and B. fungosa (Santalales) (Leebens-Mack et al., 2019), S. himalayana (Malpighiales) (Cai et al., 2021), C. australis (Solanales) (Sun et al., 2018), hemiparasitic S. asiatica (Lamiales) (Yoshida et al., 2019), S. parasitica var. graciliflora (Santalales), S. album (Santalales) (Dasgupta et al., 2019), M. oleifera (Santalales) (Xu et al., 2019) and their close relatives. The genomes of Balanophora plants were obtained from China National GeneBank DataBase (CNGBdb) (CNP0003054). We employed iTAK¹ and TFplantdb² to search the MADS-box genes from genomes analyzed in this study. Then, all obtained MADS-box gene sequences of B. subcupularis and B. fungosa var. globosa were further analyzed in SMART,³ Pfam,⁴ CDD⁵ databases to verify the existence of MADS domain.

A total of 829 MADS-box protein sequences from a total of 14 species and one MADS-box gene of Aquilegia coerulea as an outgroup were aligned using MAFFT software,⁶ and the phylogenetic tree was constructed by IQ-TREE⁷ program based on the maximum likelihood (ML) method. The parameters were set to −m GTR + R, −bb 1000, −alrt 1000, and the best alternative model is determined by model testing software. The tree was visualized by Evolview⁸ software.

According to the annotation information of B. subcupularis and B. fungosa var. globosa, the intron and exon composition of MADS-box gene were mapped by GSDS 2.0 (Gene Structure Display Server).⁹ The software MEME¹⁰ was used to search for the conservative motifs in the MADS-box genes of Balanophora. The parameter was as follows: the maximum number of motifs was 20, and the motifs width was set to 6–200 amino acids. The motifs were annotated by the Pfam database (see text footnote 4). The number of introns of species analyzed in this study was calculated according to the annotated information. The box plots of the intron length of B. subcupularis, B. fungosa var. globosa, and reference species were plotted using R packages.

We identified MADS-box gene loss primarily according to the phylogenetic tree at the ortholog level with Arabidopsis thaliana and Vitis vinifera. We divided the phylogenetic tree into ortholog groups and allowed the existence of paralogs of Arabidopsis in the same group, each ortholog group was assigned according to Supplementary Figure 2, and the primary identification of gene loss was identified based on it. For those absent genes in Balanohpora, orthologs from Arabidopsis and V. vinifera were used to predict against the assembly again (e-value threshold of 1e-5, coverage threshold of 0.25), the new annotations were checked manually including conserved domains, premature stop codons, and the transcriptome evidences. Then the manually confirmed gene models were used to construct a new phylogenetic tree, and gene loss was identified at the ortholog level again, and the absent genes in each orthologs group were considered as gene loss.

For those important flowering regulation pathways, we further identified the gene loss through orthologs modified from Sun et al. (2018). First, BLASTP was used for all-to-all proteins alignment with e-value of 1e-5, after that, we employed OrthoFinder (v.2.3.7) to cluster genes with a Markov inflation index of 1.5 and a maximum e-value of 1e-5; next, the absent genes in parasitic plants will be searched against their genomes with homologs from Arabidopsis and V. vinifera (e-value threshold of 1e-5, coverage threshold of 0.25), and new annotated gene models will be included for the next step; then, phylogenetic trees were constructed with IQ-TREE in each orthogroup generated in the last step; each tree was divided into subgroups, and each subgroup must be included at least one Arabidopsis gene; finally, the gene numbers for every species in each subgroup were counted as the ortholog gene number.

The expression profiles of flowering regulate genes in different tissues during the development of Balanophora were explored using transcriptome data downloaded from CNGBdb (CNP0003054). The expression data included three stages, stage 1: samples containing different developmental stages, LC21-YT1-3 (tuber sizes < 8 mm, duplicate: LC21-YT1, LC21-YT2, and LC21-YT3), LC21-YT4-6 (tuber sizes: 8∼15 mm, duplicate: LC21-YT4, LC21-YT5, and LC21-YT6), stage 3 samples without visible inflorescence tissues, tuber and host root were collected and named as LC24-YT1-3 (duplicate: LC24-YT1, LC24-YT2, and LC24-YT3), stage 3 samples with grown inflorescence were separated into different tissues, including tuber LC24-T3-32 (duplicate: LC24-T3 and LC24-T32), male inflorescence LC22-MF1-3 (duplicate: LC22-MF1, LC22-MF2, LC22-MF3), female inflorescence LC24-FF1-3 (LC24-FF1, LC24-FF2, LC24-FF3), inflorescence stem LC24-S1-2 (duplicate: LC24-S1, LC24-S2), and bracts (LC23-BR1). The heatmap was generated by taking the average of these duplicates.

Because the tuber of Balanophora also contains tissues from the host, we carefully classified the reads into two distinct species to avoid possible contaminations. First, the gene set of the Balanophora and its host were decontaminated (removal of genes potentially from fungi, bacteria), and then combined. Next, bowtie2 was used to map the high-quality reads from each tissue with the combined gene models of Balanophora and its host. Those reads which could be mapped into both Balanophora and its host were removed. Finally, the mapped reads were divided into two species.

Using the Cufflinks pipeline¹¹ the FPKM (Fragments Per Kilobase of transcript per Million mapped reads) value of genes in each tissue was calculated, and low-expressed and non-expressed genes were filtered out (FPKM < 10). There were three replicates for each tissue, and we took the average FPKM of the three replicates for the expression analysis. The expression levels were visualized using heatmap tools with FPKM data in each tissue in Hiplot.¹² DESeq2 was used to calculate modest estimates of folding changes and dispersion of RNA-seq data, and log2 (fold change) >1, FDR < 0.05 was considered to be a differentially expressed gene.

To study the evolution of MADS-box gene family in parasitic plants, we identified the MADS-box genes from parasitic plants and their autotrophic relatives (Figure 1B). Totally, 25 and 23 candidate genes were identified as MADS-box genes of B. subcupularis and B. fungosa var. globosa, respectively (Figure 1B). Then we named them BsubMADS1-25 and BgloMADS1-23 (Supplementary Table 2). Compared with Arabidopsis (108) and Vitis (81), we found that the gene number of this family in free-living plants was about twice that in hemiparasitic plants and three-five times that in holoparasitic plants, indicating that the number of MADS-box genes showed a stepwise loss trend with the increasing levels of parasitism. We found that the MADS-box gene family contracted sharply in holoparasitic B. subcupularis, B. fungosa var. globosa, S. himalayana, and their number of MADS-box genes is even fewer than Amborella trichopoda (36) (Albert et al., 2013). So, we inferred that the MADS-box genes of holoparasitic plants are drastically lost during the evolutionary process.

To investigate the evolutionary relationship of MADS-box genes between parasitic plants, we constructed a phylogenetic tree that included autotrophic and parasitic plants (Figure 2). They were categorized into two types: type I (Mα, Mβ, and Mγ) and type II (MIKC^C and MIKC*) (Supplementary Table 3), which were further subdivided into sixteen subgroups following the previous classification (Smaczniak et al., 2012; Fatima et al., 2020). Based on the phylogenetic trees, ten and eight genes were classified as type I in B. subcupularis and B. fungosa var. globosa, respectively, and 15 genes were classified as type II in B. subcupularis and B. fungosa var. globosa, respectively. There are at least 21 MADS-box clades in the Most Recent Common Ancestor (MRCA) of extant angiosperms (Albert et al., 2013). While in our study we found that six clades were totally absent in Balanophora and seven clades were absent in Sapria, during which four clades were convergently lost (Supplementary Table 4).

We noticed several losses of whole subfamilies in B. subcupularis and B. fungosa var. globosa, such as the AGL17, AGL15, AGL12, AGL6, and GNETUM GNEMON MADS13 (GGM13) clades, among which AGL6/12/15/17 function in the flowering transition, and TRANSPARENT TESTA16 (TT16) gene of GGM13 subfamily involved in seed pigmentation and embryo development (Erdmann et al., 2010). Most genes of AGL12 and AGL17 subfamilies are mainly involved in root development in addition to flowering (Tapia-López et al., 2008; Puig et al., 2013), which may be consistent with the lack of typical roots in Balanophora. It is worth noting that several genes classified in type II only contained the MADS domain confirmed by domain identification against Pfam database (Supplementary Figures 1, 2). Phylogenetic analysis showed that these non-K domain genes were more closely related to type II, suggesting that they may have lost their K domain during evolution. Further evidence in the following part such as intron number also supported the classification of these genes in type II.

To understand the structural diversity of MADS-box genes, we analyzed the intron and exon structure of these genes. We found that both Balanophora showed bimodal distribution between type I and II genes like other species (Wang et al., 2017, 2019; Fatima et al., 2020), significant differences in the number of introns were observed in the MADS-box genes of type I and type II in both Balanophora species (Figure 3). The number of the type II gene varied greatly, containing introns ranging from 0 to 8, among which 67% of the type II genes had at least five introns in both Balanophora. In addition, we also analyzed the number of introns of MADS-box genes of autotrophic plants and found that the intron number of type I genes was generally less than that of type II genes (Supplementary Table 5), suggesting that intron number in two types of the MADS-box gene family are conserved during evolution in flowering plants.

Intriguingly, there were very long introns inserted into the MADS-box genes of Balanophora, for example, the first intron of FLC homologue genes BsubMADS15 and BgloMADS12 are 64 kb and 56 kb, respectively, which are longer than other species (Supplementary Figure 3A). Moreover, the maximum intron length of MADS-box gene in B. subcupularis, B. fungosa var. globosa, and S. himalayana were longer than in other species (Supplementary Figure 3B). Furthermore, we analyzed the eight non-K domain genes of B. subcupularis and B. fungosa var. globosa, which were classified into type II on the phylogenetic tree. Previous studies showed that the length of the 1–6 exons of type II genes are conserved (Johansen et al., 2002). The average length of the first intron of the eight non-K domain genes in two Balanophora (197bp) was highly similar to that of type II (188bp), which was markedly smaller than that of type I genes (647bp) (Supplementary Table 6). Therefore, the results further proved the reliability of the classification of MADS-box non-K domain genes in the phylogenetic tree.

Next, we verified the conserved motifs of MADS-box genes in B. subcupularis and B. fungosa var. globosa using MEME program, then annotated the obtained motifs employing Pfam database. As a result, a total of 15 conserved motifs were identified and named motifs 1–15 (Supplementary Table 7). As shown in Figure 3, genes in the same family tend to have common patterns of conserved motifs, especially the type II gene whose domains are more conserved than type I, such as subgroups SOC1, SEP, AP1/FUL, AP3/PI. The MADS-box genes of B. subcupularis and B. fungosa var. globosa both contained motifs 1, 2, and 3, which were the most typical MADS domain according to Pfam database search. Motif 4, 6, and 10 were verified to be K-box domain, which was another conserved domain, and all MIKC^C type genes contained the K domain except BsubMADS14, BsubMADS16. We also observed that motif 7 represented the I domain that existed in all subgroups except for the subgroups AP3/PI, Mγ, Mβ of B. subcupularis and B. fungosa var. globosa. Besides, motif 8, 9, and 14 could be found in Arabidopsis, but absent in two Balanophora. Though the number and lengths of introns are variable in Balanophora, the motifs are conserved between Balanophora and Arabidopsis.

The MADS-box gene family showed stepwise contraction from hemiparasitic to holoparasitic plants (Figure 1B). We analyzed the gene loss of nine parasitic plants, most MADS-box genes of parasitic plants showed a decreasing trend in each subfamily compared with autotrophic plants (Figure 4A). In addition, there was convergent genes loss among different lineages of parasitic plants, during which the most remarkable is the loss of subfamilies AGL6, AGL12, AGL15, AGL17 in Balanophora and S. himalayana, AGL12 and AGL17 were also lost in C. australis. There was some species-specific gene loss, such as the loss of SVP in C. australis, and the loss of SOC1 in S. himalayana (Figure 4A). Hemiparasitic plants showed mild gene loss and most of them were also lost in holoparasites.

Based on the gene function of their homologies in Arabidopsis, we speculated that AGL12 and AGL17 genes may be related to the degradation of root and photoperiod regulation of flower transition in Balanophora and S. himalayana (Han et al., 2008; Tapia-López et al., 2008; Puig et al., 2013; Shu et al., 2020), they also convergently lost AGL15 and AGL6 genes, among which AGL15 genes play roles in regulating flowering through photoperiod and somatic embryo development (Perry et al., 1999) and AGL6 genes are involved in floral meristem regulation, development of floral organs, and ovule and seed, and AGL6 genes also have possible roles in the development of male and female germline and gametophyte (Dreni and Zhang, 2016). Two Balanophora also lost AGL63 and TRANSPARENT TESTA16 (TT16) genes in GGM13 subfamily, which are involved in fruit development, seed formation, and embryo development (Erdmann et al., 2010). In conclusion, parasitic plants showed convergent and functional-biased gene loss in MADS-box gene families, which may be related to the parasitic lifestyle.

There are several genes related to floral meristem identity, including AP1, FUL (AGL 8), CAL (CAULIFLOWER), AGL24, SVP, SOC1, and LEAFY (LFY) (Gregis et al., 2009; Grandi et al., 2012). We observed that two Balanophora species only lost FUL gene, however, AP1/CAL, and SOC1 genes were lost in S. himalayana (Supplementary Table 8). FUL gene promotes the identity of early flower meristem in coordination with AP1 and appears to be partially redundant to the function of AP1 (Ferrándiz et al., 2000). AP1, AGL24, and SVP control floral meristem identity redundantly by inhibiting the expression of class B, C, and E genes (Gregis et al., 2009). SOC1 regulates the expression of LFY, and LFY is a non-MADS-box gene that links floral induction and floral development (Lee et al., 2008). These results showed that Balanophora and S. himalayana showed different levels of gene loss related to floral meristem, and they may lose those functional redundant genes which were unnecessary for them.

ABCE model genes coordinate together to determine floral organ identity (Soltis et al., 2007; Murai, 2013). Two Balanophora species retained all ABCE model genes (Figure 4C). While no AP1 and no SEP3 homologs were found in S. himalayana, suggesting that the A function does not rely on the function of AP1 homologs or that the perianth has an entirely different contribution. Another holoparasite C. australis also lost AP1 in class A, and PI in class B. The result suggested that holoparasites retained different levels of the ABCE genes, these may be a consequence of the highly modified floral structure or employing different flower identity genes in them.

Flowering is essential to plant reproduction and can be regulated by a variety of pathways, such as photoperiod, circadian clock, vernalization, temperature, and autonomous, etc. (Liu et al., 2015). Previous studies have illustrated that there were at least 16 MADS-box genes involved in the regulation of flowering transition in Arabidopsis, but in our study, only three flowering signal integrators SOC1, SVP, FLC were found to be retained in Balanophora, and most genes were lost, for example, AGL15, AGL17, AGL18, MAF1, MAF3/4/5 genes involved in photoperiod pathway (Adamczyk et al., 2007; Han et al., 2008; Kim and Sung, 2010), AGL19 and MAF2 genes associated with the vernalization pathway (Ratcliffe et al., 2003; Schönrock et al., 2006), AGL28 gene involved in autonomic pathway (Yoo et al., 2006), the SOC1-like genes AGL42, AGL71, AGL72 appeared to act through gibberellin-dependent pathway in Arabidopsis (Dorca-Fornell et al., 2011; Figure 4B and Supplementary Figure 4). We observed that the subfamilies AGL12, AGL15, AGL17, and AGL6 were convergently lost in B. subcupularis, B. fungosa var. globosa, and S. himalayana, while the subfamilies of SOC1 and FLC were specifically lost in S. himalayana, and GGM13 was specifically lost in Balanophora (Figure 4A). Among them, SOC1 gene is an integrator of flowering signals involved in photoperiod, temperature, age, and gibberellin regulatory pathways in Arabidopsis to promote flowering transition, while FLC gene is a repressor of flowering, in part, it retards flowering by inhibiting the expression of SOC1. Another integrator FLC are involved in the vernalization and autonomic pathways (Michaels and Amasino, 1999, 2001; Michaels et al., 2005).

Cuscuta australis is another holoparasitic plant without roots and leaves that also lost many MADS-box genes. We observed that the subfamily SVP was specifically absent in C. australis (Figure 4A), which are involved in flowering time regulation in response to temperature changes by controlling the expression of the FT gene (Lee et al., 2007). Previous research has manifested that C. australis eavesdrops FT signals of the host to control flowering (Shen et al., 2020), taken together, we predicted that Cuscuta may have totally lost this SVP-FT signal pathway. In addition, it also lost subfamilies such as AGL17, AGL12, FLC, and GGM13 (Figure 4A), which showed convergent gene loss with holoparasitic plants. These indicated that holoparasites dramatically lost the flowering transition genes in MADS-box, Balanophora and S. himalayana showed a more severe loss.

Furthermore, we found that the gene loss in semi-parasitic plants with relatively intact body structures was not as pronounced as in holoparasitic plants. For example, S. album only lost subfamilies AGL12 and GGM13, S. parasitica lost AGL17 and GGM13 (Figure 4A). Loss of function in these subfamilies may be compensated by other genes, such as AGL17, a gene associated with root development, was functionally redundant with AGL12 (Tapia-López et al., 2008; Puig et al., 2013).

We further explored the conserved subfamilies in B. subcupularis and B. fungosa var. globosa and discovered that the number of retained genes was decreased and only one gene was kept in some subgroups. For example, the subgroups SVP, SOC1, AP1/FUL, and FLC all retained only one gene in each Balanophora species, while there were more than two copies in Arabidopsis in each subgroup (Figure 4A). We concluded that the redundancy of MADS-box genes in holoparasitic plants was greatly reduced, and they would abandon sequence similar genes and retain the diversity of protein families.

We further investigated genes involved in flower regulation pathways, and found that two Balanophora plants and S. himalayana lost the circadian clock and photoperiod pathways (Figure 5A), comprising the important genes CONSTANS (CO), TOC1, LHY, CCA1, EARLY FLOWERING 4 (ELF4), Two-component response regulator-like APRR9 (PRR9) which are not MADS-box genes (Matsushika et al., 2000; Mizoguchi et al., 2002; McWatters et al., 2007; Nakamichi et al., 2010). However, S. himalayana also some species-specific loss, including FT/FD, SOC1, and FLC functioned as the floral integrator, the FRIGIDA (FRI) complex functioned in vernalization pathways (Choi et al., 2011; Figure 5A). Compared with holoparasitic plants, hemiparasitic plants kept all crucial flowering regulatory genes (Figure 5B). Based on these results, we inferred that Balanophora and Sapria may have lost the photoperiod and circadian clock pathways to control flowering.

To further confirm the functions of MADS-box genes in Balanophora, the transcriptome data from eight different tissues at different developmental stages were analyzed. In the MADS-box genes of B. fungosa var. globosa, the expression of eleven genes could be detected. The other genes cannot be identified because their expression levels were too low, or they have specific expression patterns, which were not collected in this study, it is also possible that they have lost their functions as they adapt to parasitic life. Several genes showed specific expression patterns (Figure 5C), in B. fungosa var. globosa, genes AP1 (BgloMADS23), SEP1 (BgloMADS21), AG (BgloMADS19), STK (BgloMADS18) were highly expressed in the female inflorescence, and AG (BgloMADS19) may be involved in the development of carpel, and STK (BgloMADS18) in controlling the development of ovules of the female flower. Furthermore, PI (BgloMADS15), AP3 (BgloMADS14), AGL33 (BgloMADS11), and SVP (BgloMADS17) were highly expressed in the male inflorescence, and AG (BgloMADS19) was also expressed in male flower, suggesting that PI (BgloMADS15), AP3 (BgloMADS14), and AG (BgloMADS19) may be involved in the development of stamen. The expression pattern of these genes was consistent with the expression pattern of the ABCE model gene. In addition, we noticed that gene FLC (BgloMADS12) with extremely long introns showed high expression levels in the tuber and low expression levels in inflorescence (Figure 5C), which was similar to the expression pattern of FLC in Arabidopsis (highly expressed in the vegetative apex and root tissue). These findings revealed that the expression patterns of MADS-box genes in B. fungosa var. globosa were similar to other species, indicating that the function of these genes may be still conserved.

Because the vascular systems in the tuber of Balanophora contained both host and Balanophora tissues, the intimate connection may facilitate the usage of host signals, especially the long-distance mRNA movement between host and parasites (Kim et al., 2014). We investigated the expression levels of host genes related to regulating flowering in the tuber of Balanophora, and the results showed that several flowering regulation genes from the host were upregulated in the tuber compared to its own root, including (i) genes that can regulate FT or FLC levels transcriptionally or epigenetically, for example, AHL22, EBS, and bHLH63 function in the regulation of FT expression (Piñeiro et al., 2003; Yun et al., 2012; Liu et al., 2013); (ii) CO depression genes including MIP1A and CDF3 (Fornara et al., 2009; Graeff et al., 2016); (iii) Serine/threonine-protein kinase WNK1 regulates flowering time by modulating the photoperiod pathway (Wang et al., 2008; Supplementary Figure 5). Based on these lines of evidence, we speculated one possibility that Balanophora may take advantage of transcripts or proteins from the host plant to retain functionality in pathways of photoperiod regulation of flowering and has not lost associated regulatory pathways completely.

A phylogenetic tree of the fourteen species was built and divided into sixteen subgroups based on the taxonomy of Arabidopsis (Smaczniak et al., 2012). Phylogenetic analysis revealed that there were eight non-K domain genes clustered into type II. This phenomenon was also found in MADS-box genes of rice, apple, and pear (Arora et al., 2007; Tian et al., 2015; Wang et al., 2017), indicating that it might be common in angiosperms. Studies have shown that MIKC genes are conserved in the length of the 1–6 exons (Johansen et al., 2002), so we analyzed the length of the first exon of these eight genes and found that the length was closer to that of type II genes but different from type I genes, among them, BgloMADS11 was highly expressed in inflorescence, indicated that this gene might be functional. This result illustrated that some MADS-box genes are probably experiencing K domain loss, and at least several genes without K domain are still functional.

Reports on the evolutionary mechanism of phenotypic adaptation indicated gene loss can be beneficial by providing an evolutionary mechanism for phenotypic adaptation (Albalat and Cañestro, 2016; Sahu et al., 2019; Wang et al., 2021; Li et al., 2022). Studies on the flowering mechanism of C. australis show that in order to ensure flowering, dodder eavesdropped on the flowering signal of the host and kept the synchronization of flowering (Sun et al., 2018; Shen et al., 2020). In this study, severe gene loss of flowering regulation genes was also found in B. subcupularis, B. fungosa var. globosa, and S. himalayana, especially the circadian clock and photoperiodic flowering pathways, but there were several photoperiodic flowering response genes from the host plant are upregulated in the chimeric tuber, speculating that Balanophora may take advantage of transcripts or proteins from the host to retain functionality in pathways of photoperiod regulation of flowering, instead of losing the regulatory pathway completely. It is probably a convergent strategy for parasitic plants to steal flowering signals from the host (Shen et al., 2020). However, in hemiparasitic plants, almost all the flowering regulation pathways were retained. Dramatic loss of flower regulation genes in holoparasites analyzed in this study indicated reduced robustness in parasites with extremely degraded structure and high levels of dependence on the host.

Redundant genes generally have a genetic compensation function (i.e., the loss of one gene can be compensated by another with overlapping functions and expression patterns), and probably establish a genetically robust system for adaption (Wagner, 2000; Dean et al., 2008; Rutter et al., 2017; Yang et al., 2022). Studies have indicated that the reduction of gene expression after replication promotes the long-term maintenance of duplicate genes and functional redundancy (Qian et al., 2010). Several subgroups showed genetic redundancy in MADS-box, for example, the AP1/FUL subgroup, including AP1/CAL, FUL, and AGL79, among which AP1/CAL and FUL played important roles in the transition to floral meristem, floral meristem development, perianth, or fruit development (Gu et al., 1998; Ferrándiz et al., 2000; Acri-Nunes-Miranda and Mondragón Palomino, 2014), only one member could be identified in Balanophora (retained AP1/CAL) and Sapria (retained AGL79), respectively (Figure 2). In FLC and SOC1 subgroups, six members were confined to Arabidopsis, respectively, and only one member in each subgroup was identified in Balanophora, with a complete absence in Sapria. The reduction of redundancy in holoparasites, together with the loss of other regulatory genes, indicated that holoparasites may employ a relatively simpler flowering regulation system compared to autotrophs. The simplicity may be due to the utilization of signals directly from the host, or the highly modified structure, or the reduced demand for additional robustness when hosts offer stable environments.