Gene networks orchestrated by MeGI: a single-factor mechanism underlying sex determination in persimmon.

Separating male and female sex organs is one of the main strategies used to maintain genetic diversity within a species. However, the genetic determinants and their regulatory mechanisms have been identified in only a few species. In dioecious persimmons, the homeodomain transcription factor, MeGI, which is the target of a Y chromosome-encoded small-RNA, OGI, can determine floral sexuality. The basic features of this system are conserved in the monoecious hexaploid Oriental persimmon, in which an additional epigenetic regulation of MeGI determines floral sexuality. The downstream regulatory pathways of MeGI remain uncharacterized. In this study, we examined transcriptomic data for male and female flowers from monoecious persimmon cultivars to unveil the gene networks orchestrated by MeGI. A network visualization and cistrome assessment suggested that class-1 KNOTTED-like homeobox (KNOX)/ovate family protein (OFP)/growth regulating factors (GRFs) and short vegetative phase (SVP) genes mediate the differences in gynoecium and androecium development between male and female flowers, respectively. The expression of these genes is directly controlled by MeGI. The gene networks also suggested that some cytokinin, auxin, and gibberellin signaling genes function cooperatively in the KNOX/OFP/GRF pathway during gynoecium differentiation. Meanwhile, SVP may repress PI expression in developing androecia. Overall, our results suggest that MeGI evolved the ability to promote gynoecium development and suppress androecium development by regulating KNOX/OFP/GRF and SVP expression levels, respectively. These insights may help to clarify the molecular mechanism underlying the production of unisexual flowers, while also elucidating the physiological background enabling a single-factor system to establish dioecy in plants.

Introduction

Sexuality is the main strategy used to maintain genetic diversity within a species. By contrast with animals, most angiosperms are hermaphroditic, which is the ancestral state of flowering plants. However, some angiosperms have evolved as monoecious species, in which each individual produces distinct male and female flowers, or as dioecious species, which produce male and female flowers on different individuals (Renner, 2014). Dioecy is thought to be often associated with the presence of sex chromosomes, which include genetic determinants of sex, although only a small part of the whole dioecious species have been assessed so far (Ming et al., 2011; Renner, 2014; Charlesworth, 2016). Recent advances in genomics‐based research have helped to reveal the structure of sex chromosomes in some dioecious plants (Liu et al., 2004; Ming et al., 2011; Wang et al., 2012; Charlesworth, 2016; Kazama et al., 2016; Harkess et al., 2017; Muyle et al., 2017). Additionally, a few genetic determinants of sex have been detected on sex chromosomes in some species, including persimmons (Diospyros spp.) (Akagi et al., 2014), garden asparagus (Asparagus officinalis L.) (Harkess et al., 2017), and kiwifruit (Actinidia spp.) (Akagi et al., 2018). The recent discovery of two Y chromosome‐encoded sex‐determining genes in garden asparagus, aspTDF1 as the male‐promoting factor (M) and SOFF as the female‐suppressing factor (SuF), directly supports the ‘two‐mutation model’, which is a representative framework for the evolution of dioecy (Charlesworth and Charlesworth, 1978). The expression of a Y chromosome‐encoded sex‐determining gene identified in kiwifruit (Akagi et al., 2018), Shy Girl, can suppress female functions, which is also consistent with the two‐mutation model. On the other hand, in persimmons, a single gene located on the Y chromosome might be sufficient for determining sexuality. The OGI gene on the Y chromosome is a non‐coding RNA gene that produces a small‐RNA, and is a genetic determinant of sex in persimmons, while its autosomal counterpart, MeGI, is targeted by the OGI small‐RNA, and is thought to be the integrator of sex expression (Akagi et al., 2014). As Oriental persimmon (Diospyros kaki) evolved into a hexaploid species, its dioecious sex determination system was transferred into a more plastic system, including monoecious individuals, which are genetically male because they carry a Y chromosome (Akagi et al., 2016; Henry et al., 2018). Although the fundamental regulatory pathways in the monoecious individuals are likely to be identical to those of diploid dioecious Diospyros species, OGI is substantially silenced by a SINE‐like insertion in the promoter region (Akagi et al., 2016). By contrast, the epigenetic conditions of the MeGI promoter region and the resulting MeGI expression level are sufficient for determining the sex of each flower on monoecious trees. This implies that MeGI is the single integrator of sexuality in persimmons (Henry et al., 2018). Nevertheless, the molecular pathways underlying this integration by MeGI that is essential for androecia and gynoecia development remain uncharacterized.

Regarding the factors affecting plant sex expression, phytohormones have traditionally been considered to play important roles, although the effects are likely to differ across plant species (Grant et al., 1994). In particular, cytokinin signals are thought to be important for gynoecium development, for which the responsible molecular mechanisms have been well characterized in hermaphroditic model plant species such as Arabidopsis thaliana (Marsch‐Martínez et al., 2012). The treatment of male flowers with exogenous cytokinins often induces the development of gynoecia in some dioecious or monoecious plants, such as wild grape (Vitis amurensis) (Wang et al., 2013), kiwifruit (Actinidia spp.) (Akagi et al., 2018), and Oriental persimmon (D. kaki) (Yonemori et al., 1993). Genes encoding regulators of sex expression have recently been gradually unveiled. In Silene latifolia, components of the CLV‐WUS and CUC‐STM pathways are reportedly upregulated in a bisexual mutant that was putatively derived from a SuF‐disrupted male plant, suggesting that the Y chromosome‐encoded SuF in this species can regulate these pathways during the repression of gynoecium development (Koizumi et al., 2010). Similar pathways are also likely to affect gynoecium development in kiwifruit (Akagi et al., 2018). In garden asparagus, AMS, MS2, LAP3 and LAP5, which are genes involved in a late pollen development stage, are male‐biased genes whose expression can be influenced by aspTDF1 (Harkess et al., 2015, 2017). Additionally, the ABCDE model is well known for its role in the specification of floral organs, in which B type genes, APETALA3 (AP3) and PISTILLATA (PI), can specify androecia development along with a C type gene, AGAMOUS (AG), and E type SEP genes (Yanofsky et al., 1990; Bowman et al., 1991; Weigel and Meyerowitz, 1994; Mizukami and Ma, 1997; Rijpkema et al., 2010). During early flower development in A. thaliana, short vegetative phase (SVP) and suppressor of overexpression of constans1 (SOC1), which were originally identified as flowering time‐related genes, may encode repressors of class B genes, AP3 and PI, and a class C gene, AG (Wagner, 2008; Gregis et al., 2009). Differential B and C class gene expression levels are reportedly associated with the production of unisexual flowers in dioecious Spinacia oleracea and monoecious Quercus suber L. (Pfent et al., 2005; Sobral and Costa, 2017).

In this study, we aimed to identify the gene networks involved in the MeGI‐mediated differentiation of female and male flowers by analyzing transcriptomic data sets from various developing flowers in diverse monoecious D. kaki cultivars. Co‐expression networks have recently been commonly applied to integrate the information in large transcriptional data sets (Li et al., 2015; Liseron‐Monfils and Ware, 2015). The ‘guide‐gene approach’ (Serin et al., 2016) represents one of the effective strategies for co‐expression analyses, especially when key components of a specific pathway have been identified (Itkin et al., 2013; Serin et al., 2016). In this study, we used MeGI as the guide gene (or ‘bait gene’) to analyze the co‐expression network. We also revealed the candidate gene networks directly controlled by MeGI, in which two independent core paths regulate gynoecium and androecium development. The data presented herein may be useful for elucidating the molecular mechanisms underlying the production of unisexual flowers, while also clarifying the physiological background that enables a single‐gene system to establish dioecy.

Results

Transcriptome profiles in developing flowers

The developmental stages of D. kaki androecia/gynoecia from primordia initiation to maturation were morphologically divided into four stages (Figures 1a and S1). During these development stages, MeGI expression was substantially repressed by the methylation of the MeGI promoter and the accumulation of small RNA, which occurred in a male‐specific manner (Akagi et al., 2016). We sequenced the mRNA‐Seq Illumina libraries of each male and female flower collected in stage 1 for seven cultivars, and in stage 3 for 10 cultivars (Table S1). The reads were mapped onto the reference gene sequences of the diploid Caucasian persimmon, Diospyros lotus (Dlo_r1.0, http://persimmon.kazusa.or.jp/index.html), to calculate the expression levels as reads per kilobase of transcript per million mapped reads (RPKM). A principle component analysis (PCA) was conducted to profile the expression patterns of all genes that were substantially expressed (RPKM > 1.0) from stage 1 to stage 3 in male and female flowers (Figure 1b). PC1 and PC2 represented 42.9 and 13.2% of the total variance, respectively. The PCA clearly separated stages 1 and 3. Additionally, there were no significant differences in PC1 between female and male flowers in each stage, while significant differences were observed in PC2 between the male and female flowers in stage 3 (P = 0.038). To further investigate the relationships between male and female flowers, Pearson's distance matrix was examined. The matrix revealed a strong correlation (r = 0.87 in average) between female and male flowers in each cultivar (N = 7) in stage 1, regardless of the genotype (Figure 1c). However, the matrix indicated the correlation between female and male flowers in each cultivar (N = 10) was weaker in stage 3 (r = 0.68 in average) than in stage 1. These results suggested that dynamic changes in the transcriptomes of male and female flowers occurred during the transition from stage 1 to stage 3. This tendency was consistent with the morphological characterization, in which only slight differences between male and female organs were detected in stage 1, while there were considerable variations in stage 3 (Figure 1a).

Figure 1

Transcriptomic profiles during male and female flower development.

(a) Morphological observation and definition of the developmental stages of male and female flowers. At, anther; APr, androecium primordia; Cp, carpel; GPr, gynoecium primordia; Ova, ovary; Ovu, ovule; Pe, petal; Pi, pistil; RC, rudimentary carpel; RS, rudimentary stamen; St, stamen. Scale bars: 0.5 mm for stages 1 and 2, and 1.0 mm for stages 3 and 4. (b) Characterization of gene expression dynamics in male and female flowers by a PCA of expression data for each cultivar in stages 1 and 3. Male and female flowers were not substantially separated in stage 1, while significant differences were observed for the second component in stage 3. (c) Pearson correlation matrix of the gene expression data sets for 17 male and female flower samples collected in stages 1 and 3. For the PCA (b) and Pearson correlation analysis (c), normalized expression data for 15 644 genes with RPKM > 1 were used.

Differentially expressed genes in male and female flowers

To analyze sex‐biased gene expression, some of which might be closely associated with MeGI expression, we attempted to identify the differentially expressed genes (DEGs) between female and male flowers in stages 1 and 3. We identified 1115 and 4720 DEGs [RPKM > 1, P < 0.1; edgeR test with paired option (biological replicates N = 7 × 2 and 10 × 2 for stage 1 and 3, respectively, see Materials and Methods)], and assigned putative functions according to the annotated A. thaliana genome (TAIR10) (Dataset S1). To simplify the analysis, each persimmon gene was called based on the putative orthologous genes or functions annotated in the TAIR10 database. The persimmon gene IDs are provided in Dataset S1.

In stage 1, MeGI was identified as a female‐biased gene (Figure 2a). Moreover, genes related to meristem and gynoecium development were highly expressed in female flowers (Table S2a). For example, genes in the class‐1 KNOTTED1‐like homeobox (KNOX) subfamily, including shoot meristemless (STM), BREVIPEDICELLUS/KNAT1 (BP/KNAT1), and KNAT6, which influence flower meristem and carpel development (Arnaud and Pautot, 2014), were more highly expressed in female flowers than in male flowers (Table S2a). Additionally, ovate family protein (OFP) genes, which affect ovule development, fruit shape, and secondary wall formation (Pagnussat et al., 2007; Tsaballa et al., 2011; Wang et al., 2011, 2016), were also biased toward female flowers. Meanwhile, some male‐biased genes reportedly influence stamen development. For example, the class B genes, AP3 and PI, which are indispensable for androecium development (Weigel and Meyerowitz, 1994), were categorized as highly expressed male‐biased genes (Figure 2a, Table S2a). One of the class C genes, AG, which potentially contributes to the development of androecia and gynoecia (Bowman et al., 1991), was also identified as a male‐biased gene.

Figure 2

Detection of differentially expressed genes (DEGs) between male and female flowers, and correlations with the Me expression pattern.(a, b) Distribution of the expression patterns of the DEGs between male and female flowers in stage 1 (a) and stage 3 (b). The X and Y axes correspond to the normalized expression level (RPKM) and female/male expression ratio, respectively. The DEGs (P < 0.01) are highlighted in red. The DEGs annotated with representative gynoecium‐related, androecium‐related, or (flowering‐related) meristematic functions are indicated with pink or blue circles or a green triangle, respectively. (c) Pearson correlation coefficients of functionally annotated DEGs against the MeGI expression pattern in stage 1 were calculated. Putative gynoecium‐related, androecium‐related, or meristematic genes are indicated with pink, blue, or green bars, respectively.

We expected to detect specific genes under the direct control of MeGI in stage 1, during which there were no morphological or dynamic gene expression differences between male and female flowers (Figure 1). Pearson's product‐moment correlation test between MeGI expression patterns and all transcripts revealed that some gynoecium‐related or meristematic‐related genes biased toward female flowers, such as SVP, SOC1, AGL6, class‐1 KNOXs and OFPs, were positively correlated with MeGI (Figure 2c, Table S2a). Conversely, genes biased toward male flowers exhibited a weaker correlation than the female‐biased genes (Figure 2c). One of the most negative correlations was observed between the expression levels of a representative androecium‐related gene, PI (r = −0.5).

In stage 3, we detected 4720 DEGs between female and male flowers (Figure 2b). Moreover, there was a substantial decrease in the MeGI expression level and no significant differences between male and female flowers (P > 0.1). The 2212 female‐biased genes included representative genes related to gynoecium development (Figure 2b, Table S2b). Annotations of the 2508 male‐biased genes indicated pollen development, pollen wall assembly, and stamen development were enriched functions (Table S2b).

Core gene networks correlated with MeGI expression

Next, we attempted to establish networks reflecting the relationships between the DEGs and MeGI by applying a weighted correlation network analysis (WGCNA) to identify the module (cluster) (Langfelder and Horvath, 2008) correlated with the MeGI expression pattern. To consider the genes expressed in developing flowers, we applied all transcriptomic data from stage 1 to stage 3 (N = 44; Table S1) to construct co‐expression networks.

The DEGs between male and female flowers in stage 1 (N = 1115) were first clustered into six modules, M1–M5 and ‘unclassified’ (Figure 3a,b). The M1 module included MeGI and 366 genes, which were mostly female‐biased DEGs (Dataset S2). Thus, we defined this module as the ‘female module’. Meanwhile, male‐biased genes were enriched in the M2 and M4 modules. Two genes likely involved in androecium formation, PI and AG, were nested in the M4 module, which we designated the ‘male module’ (Figure S2). A gene‐guide approach was used to assess the correlation between the expression levels of each module and MeGI (Serin et al., 2016). As expected, the highest correlation was observed between MeGI and the female module (r = 0.76). A correlation between MeGI and the male module was detected (r = 0.082) even though genes negatively correlated with MeGI, such as AG and PI, were included in this module. These results were likely due to the dilution of the specific genes significantly correlated with MeGI (Table S2, Dataset S2).

Figure 3

Clustering and networking of the DEGs in male and female flowers.(a) Clustering of the DEGs between male and female flowers in stage 1 with the WGCNA package revealed five main modules (M1–M5). For the heat map, 44 transcriptomic samples were divided into six categories: male (M) and female (F) flowers in stages 1 and 3, with an annual replication for stage 3. (b) Correlation coefficients of each module against the MeGI expression pattern were calculated with WGCNA. The M1 module (‘female module’), which includes MeGI, exhibited the highest correlation. (c) Visualization of the female module network. Genes are clustered according to their putative functions and are presented in different nodes. Additionally, MeGI and the first‐degree genes directly connected to MeGI are indicated in orange circles. The gene clusters putatively related to plant hormones and floral organ development are presented in green and pink circles, respectively. The genes annotated with other functions are indicated in blue (see Dataset S2). The size of the nodes reflect the number of edges connected to other nodes.

A co‐expression network analysis revealed that MeGI was directly connected to 18 of 366 genes in the female module. These 18 genes were divided into the following five groups: SVP, class‐1 KNOX family (KNAT1, KNAT6, and STM), OFP (OFP7 and OFP4), growth regulating factors (GRFs), and unclassified genes (Figure 3c, Table S3). We set a topological overlap measure threshold for co‐expression at 0.17 based on the distribution of the number of edges. These 18 genes were defined as first‐degree genes directly connected to MeGI (Table S3). Consistent with our Pearson correlation test results (Figure 2c, Table S2), some meristematic‐related and/or gynoecium‐related genes, such as SEPALLATA 1 (SEP1), were identified as second‐degree genes (Figure 3c). Genes related to fruit or embryo development, such as PEROXISOME BIOGENESIS FACTOR 10 (PEX10) and ARABIDOPSIS THALIANA NUCLEOSOME ASSEMBLY PROTEIN 1 (ATNAP1), were also detected in the female module, although they were not connected to MeGI (Dataset S2). Importantly, we identified genes biased toward the female flowers in stage 1 that were nested in the female module (Dataset S2), and were related to cytokinin, auxin, and gibberellin biosynthesis/signaling, including Arabidopsis response regulator 15 (ARR15) and cytokinin oxidase 7 (CKX7) (for cytokinin); carboxyl‐terminal domain (CTD), PHOSPHATASE‐LIKE 2 (CPL2), late elongated hypocotyl (LHY), and ETTIN (ETT) (for auxin); and gamma tonoplast intrinsic protein (GAMMA‐TIP), gamma tonoplast intrinsic protein 3 (GAMMA‐TIP3), and short internodes (SHI) (for gibberellin) (Nemhauser et al., 2000; Fridborg et al., 2001; Hutchison and Kieber, 2002; Yanhui et al., 2006; Ueda et al., 2008; Köllmer et al., 2014).

The SVP gene, which encodes a regulator of the floral transition from vegetative to reproductive growth (Hartmann et al., 2000), was directly connected to MeGI in the female module (Figure 3c). Additionally, SOC1 was included in the female module, but it was not directly connected to MeGI. As mentioned earlier in the manuscript, the SVP and SOC1 genes mediate early flower development in A. thaliana by repressing class B genes (AP3 and PI) and a class C gene (AG) (Wagner, 2008; Gregis et al., 2009), which affects androecium fertility. Although SVP/SOC1 and PI/AG were nested in the female and male modules, respectively, the Pearson correlation matrix based on the expression level in stage 1 indicated that SVP, SOC1, and MeGI were negatively correlated with AG and PI (Figure S3). Consistent with this result, the pairwise correlation matrix for the female and male modules revealed a distinct negative correlation between these two modules (Figure S4), suggesting the modules may be functionally connected.

In the male module, in addition to PI and AG, LIPOXYGENASE 4 (LOX4) and SWEET7 may also affect androecium development (Bock et al., 2006, Rijpkema et al., 2010; Caldelari et al., 2011). Some genes related to auxin transport, which is critical for many aspects of androecium development (Cecchetti et al., 2008), were identified in the male module (Dataset S2), including auxin transporter protein 1 (AUX1) (Yang et al., 2006) and BIG (encoding a calossin‐like protein) (Gil et al., 2001). We also detected Arabidopsis response regulator 2 (ARR2), which is related to cytokinin signaling, and GAST1 protein homolog 1 (GASA1), which influences gibberellin signaling (Aubert et al., 1998; Borner et al., 2000) in the male module.

Organ‐specificity in candidate gene expression

To examine the possibility that MeGI promotes gynoecium development and represses androecium development, RNA in situ hybridization was applied to localize MeGI expression. Substantial MeGI expression was observed in the center of the floral meristems of female flowers in stage 1 (Figure S5a), including both gynoecium and androecium primordia. However, we did not detect substantial MeGI expression in the androecium and gynoecium primordia of male flowers in stage 1 (Figure S5a), presumably because of gene silencing induced by DNA methylation and the accumulated small RNA (Akagi et al., 2016). Furthermore, we assessed the expression bias between the gynoecium and androecium regarding the genes putatively under the direct control of MeGI as described above. We also conducted an mRNA‐Seq analysis on developing gynoecia and (rudimentary) androecia in female flowers in stage 2. The expression of class‐1 KNOX and OFP genes was considerably biased toward the gynoecium, whereas GRF and SVP genes exhibited no bias toward either the androecium or gynoecium (Figure S5b). The PI expression level was significantly lower in female flowers than in male flowers as described, and was higher in the androecium than in the gynoecium of female flowers.

Cistrome assessment to identify genes directly targeted by MeGI

We applied DNA affinity purification sequencing (DAP‐Seq) (Bartlett et al., 2017) using MeGI fused to a Halo‐Tag to examine the genes and/or motifs directly targeted by MeGI. An Illumina gDNA library of the diploid Caucasian persimmon (D. lotus) cv. Kunsenshi‐male was filtered based on the affinity to the MeGI fusion protein. On the basis of DAP‐Seq data, we identified 72 746 MeGI‐ binding sites (peaks) in vitro. The DAP‐Seq reads were mapped to the D. lotus genome to characterize the accumulated recognition motifs by using MACS2 (Zhang et al., 2008). We also identified motifs using the top 2000 high‐confidence peaks, and determined that AATWATT was enriched in MeGI‐binding regions, by using MEME‐ChIP (Machanick and Bailey, 2011). This motif is similar to the binding motifs of closely related A. thaliana HD‐ZIPs, such as ATHB21, ATHB40, and ATHB53 (Figure 4a) (Khan et al., 2018). Among the genes directly connected to MeGI in the female module (Figure 3c), SVP, KNOX genes, and OFP genes were targeted by MeGI. Specifically, MeGI binds to the promoter regions and/or the introns of these genes at sites overlapping the AATWATT motif in vitro (Figure 4b). However, there was no significant peak detected in the promoter region and only one in the second intron of STM, which is also connected to MeGI in the female module (Figure 4b). These results suggest that MeGI can bind directly to the promoters of SVP, KNOX genes, and OFP genes to regulate their expression.

Figure 4

Cistrome analysis to identify the direct target motifs/genes of Me(a) Nucleotide motifs recognized by MeGI compared with the motifs bound by the other Arabidopsis thaliana HD‐ZIP1 transcription factors nested in the MeGI clade (Akagi et al., 2014). The enriched motifs were mostly identical for MeGI and three MeGI orthologs in A. thaliana, AtHB20, AtHB40, and AtHB53. (b) Ability of MeGI to bind to the genes directly connected to MeGI in the gene co‐expression network (see Figure 3c). The full‐length sequences (5′–3′) with introns and the 2‐kb regions 5′ upstream (left side) and 3′ downstream (right side) from the start and stop codons, respectively, were analyzed. For each gene, mapped read coverages in a random genomic DNA library (‘Cont’, yellow peaks) and in the DAP‐Seq library using a MeGI fusion protein (‘DAP‐MeGI’, green peaks), a MeGI‐recognizing motif (‘AATWATT’, blue bars), and gene structures (exons, red boxes) are provided. The DAP‐Seq peaks were almost consistent with the recognized motifs. Although MeGI could bind to the promoter sequences of all first‐degree genes, except for STM, the recognized motifs were especially enriched in SVP.

Gene activation by ectopic expression of MeGI in Arabidopsis thaliana

We previously proved that the overexpression of MeGI driven by the CaMV 35S promoter represses anther and petal development (Figure 5a–d). To investigate whether MeGI regulates the expression of conserved target genes in persimmon and A. thaliana, we completed an mRNA‐Seq analysis of CaMV35S‐MeGI and CaMV35S‐empty control lines (N = 3 × 2 for biological replicates, Table S4) using developing flowers (or inflorescences) collected around stage 8 (Smyth et al., 1990). We observed that MeGI was specifically expressed in the CaMV35S‐MeGI lines (Table S4). Of the A. thaliana orthologs of the candidate target genes of MeGI (Figure 3c), the expression of SVP was significantly upregulated in the CaMV35S‐MeGI lines compared with the control lines (Figure 5e). However, STM, class‐1 KNOX genes (KNAT1 and KNAT6), and OFP genes were not significantly affected by MeGI overexpression (Figure 5e, P > 0.1). This may have been due to the saturated expression of these genes because A. thaliana is originally hermaphroditic with a functional gynoecium that requires high class‐1 KNOX gene (e.g., STM and KNAT1) expression levels (Scofield et al., 2007) to develop. By contrast, PI expression was significantly downregulated in the CaMV35S‐MeGI lines compared with the control (Figure 5e). These results are consistent with our hypothetical path of MeGI‐SVP‐PI, and with the feminized phenotype of the transformed lines, similar to the pi mutant (Bowman et al., 1989; Figure 5a–d). We detected 815 DEGs (P < 0.05) between the CaMV35S‐MeGI and control lines (Dataset S3). A comparison with the DEG orthologs in male and female persimmon flowers revealed that 16.2% (24/148) of the genes in the male module, including PI, overlapped with the CaMV35S‐MeGI A. thaliana DEGs (Figure 5f). Meanwhile, only 2.2% (8/366) of the genes in the female module, but including SVP, overlapped with the CaMV35S‐MeGI A. thaliana DEGs (Figure 5e, Dataset S3). These results were consistent with the observation that SVP expression is upregulated by MeGI in persimmon, resulting in the downregulated expression of PI and other male module genes. Of the 815 DEGs, most did no overlap with persimmon DEGs, presumably because MeGI was constitutively overexpressed in the CaMV35S‐MeGI A. thaliana lines, while MeGI expression in persimmon is limited to the meristematic region (Figure S5). However, although these DEGs were annotated with 12 statistically enriched GO terms (Table S5) (Tian et al., 2017), they were not directly related to the development or differentiation of the gynoecium/androecium.

Figure 5

Transformation of Arabidopsis thaliana with CaMV35S‐Me to characterize the activated genes.(a–d) Representative phenotypes of the transgenic lines with CaMV35S‐MeGI (a–c) and the control plants (d). The MeGI‐overexpressing plants exhibited dwarfism (a) and repressed petal and anther development (b and c). An, anther; Pe, petal; Pis, pistil; Se, sepal; Sg, stigma. (e) Expression levels (RPKM) of the genes directly connected to MeGI in the co‐expression network (see Figure 3c) as well as PI and SOC1. Among these genes, only SVP and PI exhibited significant differences (P < 0.01) between the MeGI‐overexpressing (MeGI‐OX) transgenic and control lines. Bars indicate standard errors. (f) Venn diagram of the DEGs detected between the MeGI‐OX and control A. thaliana plants (gray circle in the center, N = 815), and between male and female persimmon flowers in stage 1. Genes were clustered into six modules (M1–M5 and unclassified, see Figure 3a). The male module (blue circle) had more overlapping orthologs with the A. thaliana DEGs than the other modules, which is consistent with the phenotypic changes in the MeGI‐OX lines, which mostly produced androecia. The SVP and PI genes were included in the overlapping areas for the female and male modules, respectively.

Discussion

In monoecious D. kaki, the DNA methylation pattern in the MeGI promoter, which affects MeGI expression directly or via the accumulated small‐RNA, determines floral sexuality (Akagi et al., 2016). Consequently, this species may be suitable for identifying the sex‐determination pathways controlled by MeGI during comparisons of male and female flowers from genetically identical individuals, similar to the comparative analysis of twins. Our transcriptomic data indicated MeGI was more highly expressed in female flowers than in male flowers, which is consistent with the results of previous studies (Akagi et al., 2014, 2016). Several floral organ identity genes, whose expression was synchronized with that of MeGI, were differentially expressed between female and male flowers. On the basis of a co‐expression network analysis and the identification of candidate target genes of MeGI, we were able to define two separate pathways governing androecium and gynoecium development.

Regarding androecium differentiation in male and female flowers, SVP (and possibly SOC1) encodes one of the main repressors of androecium development. Moreover, SVP expression was synchronized with that of MeGI in the male/female flowers of D. kaki and transgenic A. thaliana flowers overexpressing MeGI. Additionally, the SVP promoter region was directly targeted by MeGI in vitro. In A. thaliana, AP1 forms a dimer mainly with SVP (and/or SOC1) to repress the expression of class B and C genes during early flower developmental stages (Liu et al., 2008; Gregis et al., 2009). Consistent with this observation, our transcriptomic data for D. kaki and transgenic A. thaliana revealed a significantly negative correlation between MeGI/SVP and PI in the androecium. These results suggest that SVP can function as an important intermediate that connects the expression of MeGI and PI during androecium differentiation. Similar mechanisms likely exist in other dioecious species. Unisexual flowers have often been used to study the relationships between ABCDE‐like genes (Sather et al., 2010; Larue et al., 2013). For example, Sobral and Costa (2017) observed that in Quercus suber L., the class B genes, especially QsPISTILLATA, were predominantly expressed in male flowers, suggesting PI is mainly responsible for the sexual differentiation in developing androecia in diverse plant species.

Our results imply that during gynoecium differentiation in male and female flowers, the class‐1 KNOX, OFP, and/or GRF genes encode mediators that are directly controlled by MeGI. Their physiological functions may also be interrelated (Hackbusch et al., 2005), which is reminiscent of the signal‐mediated regulation of gibberellin and cytokinin (Frugis et al., 2001; Jasinski et al., 2005; Yanai et al., 2005). Cytokinin is essential for gynoecium development, from gynoecium initiation, embryonic development, and even fruit development (Müller and Sheen, 2007; Werner and Schmülling, 2009; Marsch‐Martínez et al., 2012). Female gametophyte arrest reportedly occurs in the Arabidopsis histidine kinase (AHP) mutants (ahk2‐7 ahk3‐3 cre1‐12) of A. thaliana (Cheng et al., 2013). Furthermore, a type‐C cytokinin response regulator gene in kiwifruit, Shy Girl, can serve as a Y chromosome‐encoded sex determinant, while the cytokinin signaling pathway genes are differentially expressed between male and female flowers (Akagi et al., 2018). Consistent with these observations, our results imply that the expression levels of three cytokinin metabolism/signaling genes (ARR3, ARR15, and CKX7) are positively correlated with MeGI and KNOX/OFP/GRF gene expression levels (Figure 3c). Furthermore, our co‐expression analysis confirmed that STM and CUP‐SHAPED COTYLEDON 2 (CUC2), which activate cytokinin biosynthesis, are in the same female module. The cytokinin‐related pathway may be important for gynoecium differentiation in male and female persimmon flowers. This hypothesis is supported by the fact that treating male persimmon flowers with cytokinin leads to the production of hermaphroditic flowers (Yonemori et al., 1993).

Auxin and gibberellin signaling genes are also involved in gynoecium differentiation in persimmon, and function cooperatively with cytokinin signaling genes (Figure 3c, Table S2). The coordinated expression of auxin‐ and cytokinin‐related genes in the female module (Figure 3c) may be explained by crosstalk between the two plant hormones and can affect primordium formation and apical‐basal patterning during gynoecium development (Marsch‐Martínez et al., 2012; Besnard et al., 2014; Zuniga‐Mayo et al., 2014). However, auxin signaling is also essential for the late stages of stamen development. A lack of auxin signaling in the androecium causes precocious pollen maturation, anther dehiscence, and irregular filament development (Cecchetti et al., 2008). We observed that auxin‐related genes, such as CPL2, LHY, and ETT, were more highly expressed in male flowers than in female flowers in stage 3 (Dataset S1, Table S2). The gibberellin signal in the shoot apical meristem (SAM) is integrated by class‐1 KNOX genes (Sakamoto et al., 2001; Jasinski et al., 2005), and can repress cytokinin signaling (Greenboim‐Wainberg et al., 2005; Fleishon et al., 2011). Thus, among the female module genes, the lower expression levels of gibberellin‐related genes in female flowers than in male flowers (Dataset S2, Table S2) are consistent with our results described above. The expression of gibberellin‐regulated genes, such as GASA1, is also mediated by auxin (Paponov et al., 2008). Gibberellin signaling is often involved in activities associated with the formation of unisexual flowers, such as the abortion of androecium primordia in female flowers in maize (Lebel‐Hardenack and Grant, 1997) and the promotion of male flower development in cucumber (Cucumis sativus L.) (Zhang et al., 2017). Our results suggest that gibberellin signaling is important for sex determination in persimmon, potentially through cooperative effects with cytokinin signals.

Although the two‐mutation model is a representative framework for the evolution of dioecy, a single‐factor model has been proposed for the dioecious sex determination system in persimmons (Akagi et al., 2014, 2016; Henry et al., 2018). However, the underlying physiological mechanism has not been characterized. The data presented herein may provide some physiological evidence for a single‐factor model, in which MeGI regulates two independent pathways for androecium and gynoecium development (Figure 6). We hypothesize that MeGI‐SVP‐PI and MeGI‐KNOX/OFP/GRF pathways are responsible for the differences in androecium and gynoecium differentiation between male and female flowers, respectively. Although the downstream genes or mechanisms in these core pathways must still be identified, we observed some commonalities in the regulation of gene/plant hormones among dioecious plants or unisexual flowers. Overall, the results of this study imply that persimmons evolved a sex determination system in which MeGI controls the separate androecium and gynoecium pathways, while the downstream mechanisms affecting the development of these organs may be similar in diverse plant species.

Figure 6

Model for the single‐factor sex determination mechanism of persimmon.MeGI can integrate two independent pathways to promote gynoecium development and repress androecium development. The class‐1 KNOX, OFP, and GRF genes can function as intermediates that activate cytokinin (CK)‐dependent pathways [and cooperative auxin (AUX)‐ and gibberellin (GA)‐related genes] to induce gynoecium development. However, SVP (and potentially SOC1) is important for repressing PI expression, thereby leading to an aborted androecium. In hexaploid D. kaki, MeGI expression is regulated by the epigenetic status of MeGI (Akagi et al., 2016). In diploid dioecious persimmon species, male individuals with a Y‐chromosome stably express OGI, resulting in an accumulation of small‐RNA targeting MeGI. In both cases, the MeGI expression level is sufficient for determining floral sexuality.

Experimental procedures

Plant materials

The developmental stages of the Oriental persimmon (Diospyros kaki Thunb.) buds/flowers were defined in previous studies (Yonemori et al., 1993; Akagi et al., 2016). In this study, we further divided the early flower development period into the following four stages based on observations of gynoecia and androecia development: stage 1 (before organ differentiation; 1–4 April 2016), stage 2 (organs are developing; 14–17 April 2016), stage 3 (organs are maturing; 28–30 April 2016, and 27–29 April 2015), and stage 4 (flowering; 10–15 May 2016). To construct mRNA‐Seq libraries, male and female flowers were separately harvested from seven monoecious cultivars in stage 1 (‘Tohachi’, ‘Okugosho’, ‘Amayotsumizo’, ‘Iwasedo’, ‘Egosho’, ‘Meotogaki’ and ‘Taiwanshoshi’), and from 10 monoecious cultivars in stage 3 (the seven cultivars mentioned before and ‘Taishu’, ‘Zenjimaru’, and ‘Kakiyamagaki’). Additionally, the pistil of female flowers and the stamen of male flowers were collected from ‘Taiwanshoshi’ plants in stage 2 for an organ‐specific mRNA‐Seq analysis.

RNA isolation and Illumina sequencing

Total RNA was extracted from flower, pistil, and stamen samples using the PureLink® Plant Reagent (Invitrogen, Carlsbad, CA, USA), after which the mRNA was isolated using the Dynabeads™ mRNA Purification Kit (Ambion, Foster City, CA, USA). Illumina sequencing libraries were prepared as previously described (Akagi et al., 2016, 2018). Briefly, cDNA was synthesized with random primers and Superscript III (Life Technologies, Carlsbad, CA, USA) followed by a heat inactivation for 5 min at 65°C. Second‐strand cDNA was synthesized in second‐strand buffer (200 mm Tris−HCl, pH 7.0, 22 mm MgCl2, and 425 mm KCl) containing DNA polymerase I (NEB, Ipswich, MA, USA) and RNase H (NEB). Samples were incubated at 16°C for 2.5 h. The resulting double‐stranded cDNA was used to prepare libraries with the KAPA Hyperplus Kit (Kapa Biosystems, Wilmington, MA, USA), while AMPureXP (Beckman Coulter Life Science, Brea, CA, USA) was used to remove fragments shorter than 300 bp. The purified libraries were quantified with a Qubit 2.0 fluorometer (Invitrogen) and then analyzed with the Illumina HiSeq 4000 system at the QB3 Genomic Sequencing Laboratory of UC Berkeley (http://qb3.berkeley.edu/gsl/). The resulting 50‐bp single‐end reads were analyzed at the Vincent J. Coates Genomics Sequencing Laboratory at UC Berkeley. Raw sequencing reads were processed using custom Python scripts developed in the Comai laboratory, which are available online (http://comailab.genomecenter.ucdavis.edu/index.php/Barcoded_data_preparation_tools), as previously described (Akagi et al., 2014).

Transcriptomic data profiling

The mRNA‐Seq Illumina reads were aligned to the reference coding sequences (CDSs) of the diploid Caucasian persimmon, D. lotus (http://persimmon.kazusa.or.jp/index.html), using the default parameters of the Burrows‐Wheeler Aligner (version 0.7.12) (Li and Durbin, 2009, https://github.com/Ih3/bwa). The read counts per CDS were determined from the aligned SAM files using a custom R script to calculate the RPKM for each gene. To examine the gene expression dynamics in female and male flowers collected in stages 1 and 3, a PCA was conducted using prcomp in R. Additionally, Pearson's product‐moment correlation coefficients were calculated for each sample, using the cor.test with the ‘pearson’ argument in R, with gene expression levels as the parameters.

Identification of differentially expressed genes

Genes that were differentially expressed between female and male flowers were detected with edgeR (Robinson et al., 2010; McCarthy et al., 2012), using the paired‐test option, an in‐house R script, as well as 7 and 10 biological replicates for stages 1 and 3, respectively. Male and female flowers from the same cultivars were paired and used as biological replicates. The DEGs were filtered according to RPKM and P‐values (RPKM ≥ 1.0, P < 0.1). Putative functions of each gene were determined with a BLASTX search of the TAIR10 database (https://www.arabidopsis.org/index.jsp). A Pearson correlation analysis was also applied for determining the correlation between the expression levels of DEGs and MeGI in each stage, using the cor.test with the ‘pearson’ argument in R.

Construction of the co‐expression network

The DEGs in stage 1 were selected to construct a gene co‐expression network based on the WGCNA package, which is a representative algorithm used for developing co‐expression networks (Langfelder and Horvath, 2008). The soft‐thresholding power for a signed network was set at 6, with a scale‐free model fitting index R 2 > 0.8. A relatively large minimum module size (30) and a medium sensitivity (deepSplit = 2) to cluster splitting were also set. In the co‐expression network, genes were represented by nodes, and the correlation values (weight) between two genes were calculated by raising Pearson's correlation coefficient. The genes in the same module were first visualized with the VisANT program (Hu et al., 2004), and only lines with a weight greater than 0.175 and 0.065 were visualized in the module with MeGI and the module with genes related to the formation of male organs, respectively. The final networks were designed with the igraph and ggplot2 packages (Csardi and Nepusz, 2006; Wickham, 2009).

DAP‐Seq analysis

The DAP genomic DNA library was prepared and the DAP reaction was completed as previously described (O'Malley et al., 2016; Bartlett et al., 2017). Briefly, the Covaris M220 ultrasonicator (with the manufacturer‐recommended setting) was used to fragment gDNA to an average size of 200 bp. The resulting fragmented gDNA was end‐repaired with the End‐It DNA Repair kit (Epicentre, Madison, WI, USA). Next, a dA‐tail was added using the Klenow fragment (3′→5′ exo‐) (NEB). The DAP‐Seq adaptor (i.e., truncated Illumina TruSeq adapter) was ligated to the fragmented gDNA with T4 DNA polymerase (NEB). Full‐length MeGI cDNA was cloned into the pDONR221 vector (Life Technologies) and then transferred to pIX‐Halo using LR clonase II (Life Technologies) to generate pIX‐Halo‐MeGi. The N‐terminally Halo‐tagged MeGI was produced using the TNT SP6 Coupled Wheat Germ Extract System (Promega, Fitchburg, WI, USA) and purified with Magne Halo‐Tag beads (Promega). In total, 50 ng DAP gDNA library was incubated with Halo‐tagged MeGI at room temperature for 1 h. The recovered library was sequenced with the Illumina HiSeq 4000 system at the Vincent J. Coates Genomics Sequencing Laboratory at UC Berkeley.

RNA in situ hybridization

The MeGI cDNA was cloned into the pGEM‐T Easy vector (Promega). Additionally, RNA probes for the sense‐MeGI and antisense‐MeGI sequences were labeled with the DIG RNA Labeling Kit (SP6/T7) (Roche, Mannheim, Germany). The RNA in situ hybridization analysis was completed as described by Akagi et al. (2018), with minor modifications. Specifically, female and male flowers collected in stage 1 were fixed in FAA (1.8% formaldehyde, 5% acetic acid, and 50% ethanol). The FAA was then replaced by a 10–30% sucrose solution series before flower samples were sliced with a CM1520 cryostat (Leica, Wetzlar, Germany) using cryofilm as previously described (Kawamoto, 2003). The tissues were sliced into approximately 10‐μm sections, and mounted on Frontier coated glass slides (Matsunami Glass Ind., Kishiwada, Japan). The tissue sections were rehydrated in an ethanol series and then incubated in a Proteinase K solution (700 U ml−1 Proteinase K, 50 mm EDTA, 0.1 m Tris−HCl, pH 7.5) for 30 min at 37°C, followed by an acetylation with acetic anhydride (0.25% acetic anhydride in 0.1 m triethanolamine solution) for 10 min. The full‐length SyGI cDNA sequence was cloned into the pGEM‐T Easy vector (Promega) to synthesize DIG‐labeled antisense RNA probes using the DIG‐labeling RNA synthesis kit (Roche, Switzerland). The probe solution including RNaseOUT (Thermo Fisher Scientific, Waltham, MA, USA) was applied to the slides, which were then covered with Parafilm. Hybridizations were completed at 48°C for >16 h. For the subsequent detection, 0.1% anti‐digoxigenin‐AP Fab fragments (Sigma‐Aldrich, St. Louis, MO, USA) were used as the secondary antibody, which were visualized with NBT/BCIP solutions.

Transformation and transcriptomic analysis of transgenic Arabidopsis thaliana

The full‐length MeGI sequence was inserted into the pPLV26 vector (Rybel et al., 2011) in an earlier study (Akagi et al., 2014). A. thaliana plants were transformed with MeGI as described by Akagi et al. (2018). Briefly, A. thaliana ecotype Columbia‐0 plants were grown at 21°C under white light (400–750 nm) with a 16‐h light/8‐h dark photoperiod. The pPLV26‐MeGI construct was introduced into Agrobacterium tumefaciens strain EHA105 by electroporation along with the helper vector pSOUP. Wild‐type A. thaliana plants were transformed using a floral‐dip method. The putative transgenic plants were screened on Murashige and Skoog medium containing 30 μg ml−1 kanamycin.

The mRNA‐Seq analysis was completed using developing flowers that were collected from three plants for each pPLV26‐MeGI transgenic and control line around stage 8 (Smyth et al., 1990). The mRNA‐Seq libraries were constructed and sequenced as described above. The Illumina reads were aligned to the reference CDSs of A. thaliana (Araport 11) (https://www.arabidopsis.org/download_files/Genes/Araport11_genome_release/Araport11_blastsets/Araport11_genes.201606.cds.fasta.gz) using the default parameters of the Burrows‐Wheeler Aligner. The read counts per CDS were calculated based on the aligned SAM files using a custom R script. Differences in gene expression levels between the pPLV26‐MeGI transgenic and control plants were detected with the edgeR package in R (version 3.0.1).

Accession numbers

All sequence data generated in this study have been deposited in the appropriate DDBJ database, with the Illumina reads for the mRNA‐Seq analysis deposited in the Short Read Archive database (BioProject ID PRJDB7688).

Author Contributions

T.A., and R.T. conceived the study. Y.H. and T.A. designed the experiments. Y.H., T.A. and T.K. conducted the experiments. Y.H. and T.A. analyzed the data. Y.H. and T.A. drafted the manuscript. All authors approved the manuscript.

Conflict of Interest

The authors declare no conflicts of interest.

Figure S1. Morphological characterization of stage 1 and stage 3.

Figure S2. Visualization of the male module network.

Figure S3. Pearson correlation matrix in the MeGI‐SVP‐SOC1 and PI/AG.

Figure S4. Pearson correlation matrix in the whole genes between the female and male modules.

Figure S5. Organ‐specificity in expression of the MeGI and the downstream candidates.

Click here for additional data file.

Table S1. List of plant materials.

Table S2. List of the representative DEGs putatively related to androecium, gynoecium, and meristem development in stage 1 (a) and stage 3 (b).

Table S3. Genes directly connected to MeGI in the female module.

Table S4. Phenotypes of the MeGI‐overexpressed and control transgenic lines.

Table S5. Enriched GO terms in the DEGs between the MeGI‐overexpressed and control in Arabidopsis.

Click here for additional data file.

Dataset S1. Information of DEGs with persimmon IDs and annotation by TAIR.

Click here for additional data file.

Dataset S2. Genes list in each module.

Click here for additional data file.

Dataset S3. List of DEGs in the MeGI‐OX and control in Arabidopsis.

Click here for additional data file.

Click here for additional data file.