Fine-mapping and candidate gene analysis of the Mcgy1 locus responsible for gynoecy in bitter gourd (Momordica spp.)

doi:10.21203/rs.3.rs-2103453/v1

Download PDF

Research Article

Fine-mapping and candidate gene analysis of the Mcgy1 locus responsible for gynoecy in bitter gourd (Momordica spp.)

https://doi.org/10.21203/rs.3.rs-2103453/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 23 Mar, 2023

Read the published version in Theoretical and Applied Genetics →

You are reading this latest preprint version

Gynoecy plays an important role in high-efficiency hybrid seed production, and gynoecious plants are excellent materials for dissecting sex differentiation in Cucurbitaceae crop species, including bitter gourd. However, the gene responsible for gynoecy in bitter gourd is unknown. Here, we first identified a gynoecy locus designated Mcgy1 using the F₂ population (n=291) crossed from the gynoecious line S156G and the monoecious line K8-201 via bulked segregant analysis with whole-genome resequencing and molecular marker linkage analysis. Then, a large S156G×K8-201 F₂ population (n=5,656) was used for fine-mapping to delimit the Mcgy1 locus into a 292.70-kb physical region on pseudochromosome MC01, which did not include any homologous gynoecy genes previously reported in Cucurbitaceae species. Furthermore, a novel gene, MC01g1681, which encodes a cytidine triphosphate synthase, was considered the underlying gene of Mcgy1 on the basis of the results of RNA sequencing, multiple genomic sequence variation and expression analyses. In addition, transcriptome analysis of S156G and its monoecious near-isogenic line S156 revealed the potential molecular mechanisms of the formation of gynoecy in bitter gourd. Our findings provide a new insight into inheritance modes of gynoecious plants in Cucurbitaceae species and a theoretical basis for breeding gynoecious bitter gourd lines by the use of molecular markers.

bitter gourd

gynoecy

bulked segregant analysis

fine-mapping

transcriptome analysis

The Mcgy1 locus responsible for gynoecy was fine-mapped. MC01g1681, encoding a cytidine triphosphate synthase, was considered as the underlying gene for Mcgy1 according to RNA sequencing, multiple genomic sequence variation and expression analyses.

In angiosperms, there are three basic types of flower sexes, namely, hermaphroditic, male, and female flowers, as well as combinations of them. Thus, there are a total of seven types of sexually reproducing plants, namely, androecious (only male flowers per plant), gynoecious (only female flowers per plant), hermaphroditic (only hermaphroditic flowers per plant), monoecious (male and female flowers per plant), andromonoecious (male and hermaphroditic flowers per plant), gynomonoecious (female and hermaphroditic flowers per plant), and trimonoecious (male, female, and hermaphroditic flowers per plant) plants (Dellaporta & Calderon-Urrea, 1993; Lewis, 1942). Generally, hermaphroditic flowers are considered the ancestors of unisexual (female or male) flowers (Barrett, 2002), and the latter are capable of cross-pollination and thus more easily promote abundant genetic diversity (Hamrick & Godt, 1996). Among these different types of traits, gynoecy has become a popular topic and focus of research because it plays an important role in high-efficiency hybrid seed production by avoiding the possibility of pollen contamination and the need for artificial emasculation, and gynoecious plants serve as excellent materials for dissecting sex differentiation, especially in cucurbit crop species (Martin et al., 2009; Zhang et al., 2020a).

Cucurbit crop species include many important vegetable and fruit species, such as cucumber (Cucumis sativus), melon (Cucumis melo), watermelon (Citrullus lanatus), and bitter gourd (Momordica charantia), and are usually considered ideal model plant species for dissecting sex differentiation because plants of all seven types of sexual reproduction can be found (Schaefer & Renner, 2011). Many researchers consider that ethylene is a main factor regulating sex differentiation in cucurbits, and the endogenous biosynthesis of ethylene is mainly determined by two genes, ACS and ACO (Byers et al., 1972; Manzano et al., 2011; Manzano et al., 2014; Zhang et al., 2017). ACS encodes a 1-aminocyclopropane-1-carboxylic acid (ACC) synthase catalysing the conversion of S-adenosylmethionine (SAM) to form ACC, and ethylene is then synthesized from ACC, which is catalysed by ACC oxidase encoded by ACO (Boualem et al., 2015). For example, a unique F locus associated with gynoecy was identified in cucumber, and previous works have suggested that monoecious plants carry only one CsACS1 gene, while gynoecious plants have not only the CsACS1 gene but also an additional copy of CsACS1 named CsACS1G (F) (Mibus & Tatlioglu, 2004; Li et al., 2020). The loss of function of the androecy (a) gene CsACS11/CmACS11, which is the hypostatic gene of F and controls female flower development, leads to androecious cucumber and melon plants (Boualem et al., 2015). CsACO2/CmACO3, the downstream gene of CsACS11/CmACS11, was identified as the second androecy-related gene, and the loss of function of this gene resulted in monoecious cucumber and melon mutants becoming androecious (Chen et al., 2016). Unlike in unisexual flowers, a monoecious (m) gene (CmACS7/CsACS2/CitACS4) is believed to arrest the development of stamens in female flowers, and in cucumber, melon, and watermelon mutants harbouring a nonfunctional m gene, the development of stamens cannot be prevented, leading to hermaphroditic flowers (Boualem et al., 2008; Boualem et al., 2009; Li et al., 2009; Tan et al., 2015; Ji et al., 2016; Manzano et al., 2016). In addition, the interaction of F and m (FFmm) results in the formation of hermaphroditic cucumber plants (Li et al., 2009).

Some studies have suggested that, in addition to the genes involved in ethylene biosynthesis, the C₂H₂ zinc finger transcription factor WIP plays a significant role in the formation of gynoecy in cucurbits. The expression of CmWIP1 (g) leads to carpel abortion and suppresses the expression of CmACS7 in male flowers, and epigenetic changes of CmWIP1 induced by the insertion of a transposon was shown to induce the expression of CmACS7 in melon, resulting in the conversion of male flowers to female thus leading to gynoecy (ggMM) (Martin et al., 2009). Similarly, chromosome translocation disrupts the expression of ClWIP1, leading to gynoecy in watermelon (Zhang et al., 2020a). The function of CsWIP1 has been shown to be the same as that of CmWIP1/ClWIP1 in cucumber via clustered regularly interspaced short palindromic repeats/CRISPR-associated system 9 (CRISPR/Cas9) technology, as no natural mutants of CsWIP1 are available (Hu et al., 2017). These studies indicated that these WIP orthologues are highly conserved in three cucurbit crop species: cucumber, melon, and watermelon.

Bitter gourd, an important member of the Cucurbitaceae family, originated in the tropical regions of Africa (Schaefer & Renner, 2010). The fruit of this species is rich in nutrients and has medicinal properties, including antidiabetic, anticancer, hypotensive, anti-obesity, antimicrobial, and antihyperlipidaemic properties (Wang et al., 2017). It is currently widely cultivated in Africa, Asia, South America, the Caribbean, parts of the Amazon Basin, and so forth (Basch et al., 2003). Like other cucurbit crop plants such as cucumber, bitter gourd plants can display many types of sexual reproduction, but only two types, monoecy and gynoecy, have been reported thus far (Kole, 2020). Previous genetic analyses have suggested that gynoecy in bitter gourd is controlled by a single recessive gene named Mcgy, gy, or gy-1 (Matsumura et al., 2014; Cui et al., 2018; Gangadhara Rao et al., 2018). Based on the genetic map, Mcgy and gy were delimited via restriction site-associated DNA sequencing (RAD-seq) to genetic regions of 12.7 cM and 42.9 cM, respectively (Matsumura et al., 2014; Cui et al. 2018), and gy-1 was delimited via genotyping by sequencing (GBS) technology to a 10.10 cM genetic interval (Gangadhara Rao et al., 2018). However, the gene responsible for gynoecy in bitter gourd remains unknown.

In this study, the candidate locus controlling gynoecy was first identified via bulked segregant analysis with whole-genome resequencing (BSA-seq), and further narrowed by molecular marker linkage analysis and fine-mapping. Through RNA sequencing (RNA-seq), multiple genomic sequence variation and expression analyses, we suggested that MC01g1681 is an important candidate gene responsible for gynoecy in bitter gourd. Our findings provide a basis for breeding gynoecious lines and help elucidate the molecular formation mechanism of gynoecy in bitter gourd.

Plant materials

Three gynoecious lines, S156G, K20-95, and K20-506, whose causal gene of gynoecy was transferred from the spontaneous mutant K44 (Cui et al., 2018), and three monoecious lines, K8-201, S060, and K20-507, were used as female and male parents, respectively, to construct three F₁ populations, namely, S156G×K8-201 F₁, K20-95×S060 F₁, and K20-506×K20-507 F₁ populations, by artificial hybridization. Three F₂ populations, S156G×K8-201 F₂, K20-95×S060 F₂, and K20-506×K20-507 F₂ populations, were then produced by artificial self-pollination using the corresponding F₁ materials. These three sets of independent genetic materials, which included three gynoecious lines, three monoecious lines, three F₁ populations, and three F₂ populations, were used to verify the inheritance of gynoecy. Furthermore, a set of S156G, K8-201, and their resulting F₂ plants was used for BSA-seq, molecular marker linkage analysis, and fine-mapping to delimit the candidate region of the gynoecy-related gene. The above six parental lines and six additional monoecious inbred lines, namely, K20-508, THMC170, FOLI112, K7-359, S121, and TR, were subjected to whole-genome resequencing and genomic sequence variation analysis to identify the candidate genes of gynoecy, and 79 monoecious inbred lines of bitter gourd germplasms (Supplementary Table S1) with distinct genetic backgrounds were used to verify the candidate genes. In addition, S156G and its monoecious near-isogenic line, S156, were used for RNA-seq and real-time quantitative polymerase chain reaction (qRT-PCR). All the plants were cultivated in the experimental field at the SCAU Teaching & Research Base in Zengcheng District, Guangzhou (23.24N, 113.64E).

Whole-genome resequencing and variant calling

DNA of all the samples was extracted from young leaves by using the modified cetyl-trimethylammonium bromide method (Porebski et al., 1997). Two DNA pools, namely, monoecy-associated pool and gynoecy-related pool, were constructed by mixing equal amounts of DNA from 30 monoecious and 30 gynoecious plants in the S156G×K8-201 F₂ population for use in BSA-seq.

A total of 14 DNA libraries for whole-genome resequencing, namely, those corresponding to two DNA pools, six parental lines, and six additional monoecious inbred lines, were prepared as follows. The qualified DNA samples were randomly fragmented by Covaris, followed by fragment selection to an average size of 200 ~ 400 bp by an Agencourt AMPure XP-Medium Kit (Thermo Fisher Scientific, USA). The selected fragments were end repaired and 3’adenylated, and the adaptors were then ligated to the ends of the 3’adenylated fragments. The products were amplified by PCR and purified by the Agencourt AMPure XP-Medium Kit. The purified double-stranded PCR products were denatured (via heating) to a single strand and then circularized by splint oligo sequences. Single-strand circular DNA was used to constitute the final DNA library, which was then sequenced on an Illumina X-ten system.

The raw data generated from the sequencing pipeline were filtered using SOAPnuke (Chen et al., 2018) by discarding reads with more than 50% adaptor sequences, low-quality reads with more than 50% bases and Phred values less than 20, and reads with 2% or more “N” bases to obtain clean data. The clean reads were aligned to the Dali-11 reference genome using Burrows‒Wheeler Aligner (Li & Durbin, 2009; Cui et al., 2020), and variant calling of single-nucleotide polymorphisms (SNPs) and insertions-deletions (InDels) was performed via the Genome Analysis Toolkit (Mckenna et al., 2010). All the variants were annotated by ANNOVAR (Wang et al., 2010).

BSA-seq

Based on the SNPs generated from the S156G, K8-201, monoecy-associated pool, and gynoecy-related pool data, the Euclidean distance (ED) and SNP index algorithms were used for BSA. The ED^2, which represented the association value, was calculated by the ED algorithm, and the candidate region controlling gynoecy was delimited where the fitted value was greater than the median + 3SD (Hill et al., 2013). The ΔSNP index, which represented the association value, was the difference value of the SNP index of two pools calculated by the SNP index algorithms by sliding window analysis with a sliding window size of 200 kb and a sliding step of 100 kb, and the fitted value within the 99% confidence interval was taken as the lower limit to delimit the candidate region controlling gynoecy. After overlapping the results from the ED and SNP index algorithms, we preliminarily mapped the final candidate region responsible for gynoecy.

Molecular marker linkage analysis and fine-mapping

To genotype each plant within the S156G×K8-201 F₂ population (n = 291), seven polymorphic molecular markers, GY_1 ~ GY_7, that were polymorphic between S156G and K8-201 and targeted the candidate region delimited by BSA-seq were developed using Primer3 Plus (https://www.primer3plus.com/index.html) according to the variations obtained from the resequencing data; the primer sequences are listed in Supplementary Table S2. PCR was performed in a 10 µL reaction mixture consisting of 5 µL of Green Taq Mix (Vazyme, Nanjing, China), 0.2 µL of forward and reverse primers (10 µmol L^− 1), 1 µL of template DNA (50–100 ng µL^− 1), and 3.6 µL of nuclease-free water. The PCR procedure was as follows: an initial denaturation cycle of 3 min at 94°C; 34 cycles of denaturation of 15 s at 94°C, annealing of 15 s at 55°C, and extension for 30 s at 72°C; and a final extension cycle of 5 min at 72°C. The SNP markers were further transformed into cleaved amplified polymorphic sequences (CAPS) markers or derived cleaved amplified polymorphic sequences (dCAPS) markers through a stationary temperature condition of 30 min at 37°C followed by digestion with the corresponding restriction enzyme (Supplementary Table S2). The PCR products or the digested products were then visualized using 6% polyacrylamide gel electrophoresis (PAGE). A molecular marker linkage map was constructed based on the genotype and phenotype of each plant within the F₂ population using JoinMap 4.0 software; the logarithm of odds (LOD) score threshold was set to 3.0, and the recombination frequency was set to less than 0.4 (Van Ooijen, 2006).

Based on the molecular marker linkage map, two flanking markers, GY_3 and GY_7 (Supplementary Table S2), of the gynoecy locus were used for screening recombinant plants among the 5,656 individuals composing the large S156G×K8-201 F₂ population. Moreover, four new polymorphic molecular markers, GY_8 ~ GY_11 (Supplementary Table S2), between the two flanking markers were developed and used to genotype each recombinant in the same manner as that described above. Then, after combining the genotype and phenotype data of all the recombinants, we finely mapped the candidate region to an even narrower region. A collinearity analysis within the candidate region between the Dali-11 genome (Cui et al., 2020) and OHB3-1 genome (Matsumura et al., 2020) was performed using MUMmer and LINKVIEW 2 (Delcher et al., 2003) to ensure the reliability of the genome assembly in this candidate region.

Multiple genomic sequence variation analysis

Three gynoecious inbred lines, namely, S156G, K20-95, and K20-506, and nine monoecious inbred lines, namely, K8-201, S060, K20-507, K20-508, THMC170, FOLI112, K7-359, S121, and TR, were used for genomic sequence variation analysis within the candidate region via TBtools (Chen et al., 2020), to identify the unique or few variations belonging to only the gynoecious lines. The candidate variation responsible for gynoecy was further genotyped by the use of the molecular marker GY_10 (Supplementary Table S2) in 79 bitter gourd inbred lines with diverse genetic backgrounds.

RNA-seq and comparative transcriptome analysis

The stem apexes of S156G and S156 at the two-leaf, six-leaf, 10-leaf, and 14-leaf developmental stages were collected and frozen in liquid nitrogen for total RNA isolation using an Eastep Super Total RNA Extraction Kit (Promega, China). Three biological replicates were included for each sample. The qualified RNA was prepared for library construction using an NEBNext Ultra™ RNA Library Prep Kit (NEB, USA) following the manufacturer's instructions. The accuracy of the library concentrations were evaluated by the qPCR method, and after which the libraries were then sequenced on an Illumina sequencing platform.

Based on the raw data generated by RNA-seq, clean data were obtained by removing the low-quality reads containing adaptor and poly-N (N ≥ 10%) sequences and then mapped to the Dali-11 reference genome using HISAT2 (Kim et al., 2015). The mapped reads were subsequently assembled by StringTie (Pertea et al., 2015) and then used for comparison with the original gene annotation information from the Dali-11 genome (Cui et al., 2020) and to optimize annotated genes or search for new genes. The Genome Analysis Toolkit and ASprofile were used for calling variants and predicting alternative splicing events, respectively (Mckenna et al., 2010; Florea et al., 2013). Differential alternative splicing analysis was then performed by rMATS (Khan et al., 2018). Fragments per kilobase of transcript per million fragments mapped (FPKM) values were applied to measure the expression level of transcripts by StringTie using the maximum flow algorithm (Florea et al., 2013). Differential expression analysis was performed using DESeq2 (Love et al., 2014), and the genes were considered differentially expressed genes (DEGs) only when the resulting P value was < 0.01, which was adjusted based on Benjamini and Hochberg’s approach. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the DEGs was performed using KOBAS software (Mao et al., 2005). To verify these DEGs, qRT-PCR was performed on a CFX384 Real-Time System (Bio-Rad, CA, United States) using Eastep qPCR Master Mix (Promega, China) according to the manufacturer’s protocols. The GAPDH gene (MC08g1155) was used as the internal control, and the sequences of the primers used are listed in Supplementary Table S3. The relative expression levels were calculated via the 2^−ΔΔCt method (Livak & Schmittgen, 2001).

Characteristics and inheritance of gynoecy

Gynoecious bitter gourd plants produce only female flowers on both the main stems and the lateral branches, although there is the occasional appearance of one or two male flowers at low nodes along the main stem (Fig. 1). Monoecious bitter gourd plants produce both female and male flowers along the main stems and lateral branches (Fig. 1). We evaluated each plant among the three sets of independent genetic materials. We found that all plants resulting from female S156G, K20-95, and K20-506 lines were gynoecious, and all plants from male K8-201, S060, and K20-507 lines were monoecious (Table 1). Like their male parents, the plants composing these F₁ generations were monoecious (Table 1). Among the plants of the three F₂ populations, 228 and 63 plants in the S156G×K8-201 F₂ population were monoecious and gynoecious, respectively, 187 and 69 plants in the K20-95×S060 F₂population were monoecious and gynoecious, respectively, and 226 and 80 plants in the K20-506×K20-507 F₂population were monoecious and gynoecious, respectively (Table 1). Chi-square tests revealed that the recorded segregation in all three F₂ populations fit the expected ratio of 3:1 (Table 1). Our findings suggested that gynoecy is a stable recessive trait controlled by a single gene named Mcgy1.

Table 1

Basic statistics of plant phenotypes based on field investigations

Generation	Monoecy	Gynoecy	Expected ratio	χ²	P
S156G	0	33	-	-	-
K8-201	30	0	-	-	-
S156G×K8-201 F₁	30	0	-	-	-
S156G×K8-201 F₂	228	63	3︰1	1.74	0.19
K20-95	0	35
S060	31	0	-	-	-
K20-95×S060 F₁	30	0	-	-	-
K20-95×S060 F₂	187	69	3︰1	0.52	0.47
K20-506	0	30
K20-507	33	0	-	-	-
K20-506×K20-507 F₁	30	0	-	-	-
K20-506×K20-507 F₂	226	80	3︰1	0.21	0.64

BSA-seq to rapidly identify the Mcgy1 locus

After filtering the raw data generated from the whole-genome sequencing of S156G, K8-201, the monoecy-associated pool, and the gynoecy-related pool data, a total of 52.06 Gb of clean data with a mean sequencing depth of 43.69×, an average Q20 value of 96.44%, and an average coverage of 97.15% (≥1×) were obtained, which indicated that these data were of high enough quality to perform subsequent analyses (Supplementary Table S4). Based on variant calling, a large number of SNPs evenly covering the entire genome were identified (Supplementary Table S5) and then used to execute BSA for identifying the Mcgy1 locus. The results revealed that the Mcgy1 locus was located on the end of pseudochromosome MC01 (hereafter referred to as simply MC01) within a 3.02-Mb region from 18,124,767 to 21,148,382 bp via the ED algorithm (Fig. 2a), and by the use of the SNP index algorithm, the Mcgy1 locus was located to within a 2.55-Mb region from 18,595,926 to 21,148,382 bp on MC01 (Fig. 2b). Overlapping the regions calculated by the ED and SNP index algorithms, we preliminarily delimited the Mcgy1 locus to within a physical distance of 2.55-Mb from 18,595,926 to 21,148,382 bp on MC01 (Fig. 2).

Classic molecular marker linkage analysis validates the Mcgy1 locus

To verify the Mcgy1 locus delimited by BSA-seq, seven polymorphic molecular markers, namely, four InDel and three CAPS markers (Supplementary Table S2), targeting the candidate region were developed and used for genotyping the 291 individuals composing the S156G×K8-201 F₂ population. A partial molecular marker linkage map was constructed on the basis of the genotype and phenotype of each plant (Fig. 3a). The results showed that there were different degrees of linkage between the seven markers and Mcgy1, with genetic distances ranging from 0 cM to 23.6 cM, and four markers, namely, GY_4, GY_5, GY_6, and GY_7, cosegregated with Mcgy1 (Fig. 3a). Notably, the CAPS marker GY_7 (21,144,148 bp) can be used as a boundary for the Mcgy1 locus because, on the basis of the gene annotation of the Dali-11 genome, there is no gene in the 4.23-kb region from GY_7 to the right end of MC01 (Fig. 3b) (Cui et al., 2020). Therefore, the Mcgy1 locus identified by BSA-seq was validated by molecular marker linkage analysis and narrowed to a 476.61-kb region from 20,667,540 to 21,144,148 bp between two flanking markers, GY_3 and GY_7 (Fig. 3b).

Fine-mapping of the Mcgy1 locus

A large S156G×K8-201 F₂ population comprising 5,656 individuals was genotyped using the GY_3 and GY_7 markers to screen recombinant plants. A total of 45 recombinants, namely, 27 gynoecious plants and 18 monoecious plants, were obtained and divided into six haplotypes (Fig. 3c) by the use of GY_4~GY_6 and four new markers, GY_8~GY_11 (Supplementary Table S2). Six markers, namely, GY_9, GY_10, GY_11, GY_5, GY_6, and GY_7, cosegregated with Mcgy1 in all 45 recombinants (Fig. 3c). Combining the genotypes and phenotypes of the 45 recombinants, we ultimately mapped the Mcgy1 locus to a 292.70-kb interval of 20,851,441~21,144,148 bp between the two markers GY_8 and GY_7 on MC01 (Fig. 3c), which reflected a region with a very low recombination rate. Additionally, the results of collinearity analysis showed that there is high similarity in the 292.70-kb region between the Dali-11 genome and the long-read OHB3-1 genome (Supplementary Fig. S1), implying that the assembly of genomic sequences within this 292.70-kb region from the Dali-11 genome is reliable. In other words, the causal gene of Mcgy1 is most likely in this 292.70-kb region. Based on the gene annotation of the Dali-11 reference genome (Cui et al., 2020), there were 29 genes in the 292.70-kb interval (Supplementary Table S6).

RNA-seq of stem apexes provides a preliminary estimate for the candidate genes of Mcgy1

The stem apexes from S156G and S156 at the two-leaf, six-leaf, 10-leaf, and 14-leaf developmental stages displaying the most active flower bud differentiation were subjected to RNA-seq. In total, 162.53 Gb of clean data with a mean Q30 value of 94.47% was obtained, and the mapping percentage of each sample against the Dali-11 reference genome ranged from 93.59% to 96.96% (Supplementary Table S7). Based on the RNA-seq data of stem apexes, we further improved the gene annotation data of the Dali-11 genome generated via transcriptome data from multiple tissues, including roots, stems, leaves, male flowers, ovaries, and fruits (Cui et al., 2020), that optimized 2,442 gene structures (Supplementary Table S8) and found 1,717 new genes (Supplementary Table S9). Four new genes, namely, MC_newGene0100, MC_newGene0101, MC_newGene0102, and MC_newGene0103, were located in the 292.70-kb candidate region delimited by fine-mapping. Hence, adding the 29 genes originally annotated from the Dali-11 genome, there were actually 33 genes in the fine-mapping region (Supplementary Table S6). Moreover, the results of differential alternative splicing analysis showed that there were no alternative splicing events in this region between S156G and S156 at any of the four developmental stages, excluding the possibility that Mcgy1 regulated gynoecy through alternative splicing (Supplementary Tables S10 and S11). In addition, within the fine-mapping interval, we detected only one SNP (20,943,203 bp) locating at 3’ untranslated region (3’ UTR) of MC01g1681, which encodes a cytidine triphosphate synthase (CTPS), and none of the variations in coding sequence of any of the 33 genes (Supplementary Table S12). In S156G and S156, notably, the transcript levels of MC01g1681 showed more obvious difference at two-leaf developmental stage; whereas of the other 32 genes, at all four developmental stages, part of them showed none or only minor difference and the rest showed awfully low or absent (Supplementary Fig. S2).

Meanwhile, BLAST analysis revealed six WIP genes in bitter gourd, namely, MC00g0064, MC01g0142, MC02g0091, MC03g0414, MC05g1109, and MC09g0449 (Supplementary Fig. S3), whose homologues lead to gynoecy in other cucurbit species, such as cucumber, melon, and watermelon (Martin et al., 2009; Hu et al., 2017; Zhang et al., 2020a). Additionally, there are nine ACS genes in bitter gourd, namely, MC01g0962, MC06g0753, MC06g0754, MC07g0325, MC08g0531, MC06g2011, MC10g0187, MC10g0527, and MC10g0702 (Supplementary Fig. S4), whose homologues result in gynoecy in cucumber (Mibus & Tatlioglu, 2004). Taken together, these results showed that the 292.70-kb candidate region does not include any WIP or ACS homologous, indicating that Mcgy1 might be a novel gene responsible for gynoecy in bitter gourd.

Multiple genomic sequence variation analysis to define the candidate genes of Mcgy1

After aligning the 12 inbred lines, namely, the six parents and six additional monoecious lines, to the Dali-11 genome, we obtained 25, 108, 139, 40, 1,293, 244, 38, 3,027, 10, 9, 120, and 3,044 variations in S156G, K20-95, K20-506, K8-201, S060, K20-507, K20-508, THMC170, FOLI112, K7-359, S121, and TR, respectively, within the 292.70-kb interval delimited by fine-mapping (Supplementary Tables S13~S24). Among these variations, only four single-nucleotide variations (SNVs), namely, SNV-1 (20,860,222 bp), SNV-2 (20,943,203 bp), SNV-3 (20,997,811 bp), and SNV-4 (21,144,642 bp), belonged independently to the three gynoecious lines, S156G, K20-95, and K20-506 (Fig. 4a). Of these SNVs, three of them, SNV-1, SNV-3, and SNV-4, are located within intergenic spacers and 2,963 bp, 22,496 bp, and 21,076 bp away from their upstream genes – MC01g_new0504, MC01g1685, and MC01g_new0512, respectively (Fig. 4b). Only two variations, SNV-1 and SNV-3, harbour downstream genes – MC01g1670 and MC01g1686 – and are 2,432 bp and 10,185 bp away, respectively, as SNV-4 is located at the end of MC01, with no downstream gene (Fig. 4b). The other SNV-2 same as the SNP detected by RNA-seq is located within the 3’ UTR of MC01g1681 (Fig. 4b). Therefore, the above six genes should be considered as the candidate genes of Mcgy1, and of which MC01g1681 should be focused on since the distribution of the four SNVs.

Expression analysis suggests that MC01g1681 is an important candidate for Mcgy1

For further define the candidate genes, we evaluated detailly expression levels of the six candidate genes in stem apexes combined with RNA-seq and qRT-PCR analyses. Among the five candidate genes adjacent to the three SNVs locating at intergenic spacers, three, MC01g1670, MC01g1685, and MC01g_new0512, showed hardly any expression levels at all four developmental stages in both S156G and S156, and two, MC01g1686 and MC01g_new0504, had no obvious expression difference at all stages between S156G and S156, except for one minor multiple (only 1.4 times) of significant difference of MC01g_new0504 at six-leaf developmental stage (Supplementary Fig. S5). Furthermore, the results of qRT-PCR also confirmed that there was no obvious expression difference both of MC01g1686 and MC01g_new0504 at all four stages between S156G and S156, especially the expression level of MC01g_new0504 at six-leaf developmental stage was slightly higher (only 1.2 times) in S156 than that in S156G but no significant difference (Supplementary Fig. S6). For the remaining candidate gene, notably, MC01g1681 had the highest FPKM value in S156G at the two-leaf developmental stage and was significantly higher in S156G than that in S156 (Fig. 4c). The transcript levels of MC01g1681 in S156G at the later three developmental stages (six-leaf, 10-leaf, and 14-leaf) were relatively stable and significantly lower than those at the two-leaf developmental stage, whereas the expression levels of MC01g1681 in S156 remained relatively stable at all four developmental stages (Fig. 4c), the results of which were subsequently confirmed by qRT-PCR (Fig. 4d). Taken together, we believed that MC01g1681 is an important candidate for Mcgy1, which might play an important role in sex differentiation in bitter gourd at the early developmental stage.

MC01g1681 can be verified in different bitter gourd germplasms using molecular marker

SNV-2 is a G-A (monoecy-gynoecy) mutation locating at 3’UTR of MC01g1681 obtained by resequencing data of the 12 inbred lines, which was further confirmed using Sanger sequencing (Fig. 4e). To verify the association of MC01g1681 and gynoecy in different bitter gourd germplasms, we designed a dCPAS marker (GY_10) to target SNV-2 by introducing a mismatched base (A) at end of reverse primer to create a Acc Ⅰ restriction enzyme site (Supplementary Table S2). A total of 79 monoecious inbred lines of bitter gourd germplasm were genotyped with the marker GY_10, the results of which showed a complete correlation with type of sexual reproduction (Fig. 4f), which implied that GY_10 developed by SNV-2 can be used to marker-assisted selection for gynoecy in bitter gourd.

Transcriptome analysis reveals the potential molecular formation mechanisms of gynoecy

Based on transcriptome data, we obtained global gene expression pattern of stem apexes of S156G and S156 at the two-leaf, six-leaf, 10-leaf, and 14-leaf developmental stages (Supplementary Table S25). Comparative transcriptome analyses revealed that there were 461, 128, 866, and 149 DEGs between S156G and S156 at the two-leaf, six-leaf, 10-leaf, and 14-leaf developmental stages, respectively (Supplementary Table S26). Among these DEGs, 10 were upregulated at all four developmental stages, and the transcript levels in S156G were significantly higher than those in S156 (Fig. 5a and Supplementary Table S26); and 11 genes were downregulated at all four developmental stages (Fig. 5b and Supplementary Table S26). These 21 core DEGs may play an important role in the formation and maintain of gynoecy in bitter gourd.

KEGG enrichment analysis was then carried out to further explore the metabolic pathway impacting the formation of gynoecy. The results of which showed that plant hormone signal transduction (ko04075) was the most significantly enriched pathway at the two-leaf and six-leaf developmental stages (Fig. 5c and Supplementary Fig. S7), whereas protein processing in the endoplasmic reticulum (ko04141) was the most significantly enriched pathway at the 10-leaf and 14-leaf developmental stages (Supplementary Figs. S8 and S9). 10 DEGs from the plant hormone signal transduction pathway at two-leaf developmental stage were selected for quantification via qRT-PCR to independently verify the RNA-seq data, and the results showed that all of them displayed significant differences in terms of their expression levels between S156G and S156, which was consistent with the results of RNA-seq analysis (Fig. 5d and Supplementary Tables S25 and S27). These findings provided crucial clues into the molecular formation mechanism of gynoecy in bitter gourd.

Gynoecious plants serve as excellent materials for increasing yields and producing hybrid seeds of Cucurbitaceae crop species. Previous studies have established that the C₂H₂ zinc finger transcription factors CsWIP1, CmWIP1, and ClWIP1 are the causal genes of gynoecy in cucumber, melon, and watermelon, respectively (Martin et al., 2009; Hu et al., 2017; Zhang et al., 2020a), and a unique F locus associated with an additional copy of CsACS1, namely, CsACS1G, is responsible for gynoecy in cucumber (Mibus & Tatlioglu, 2004). In addition, no other genes involved in gynoecy have been identified in Cucurbitaceae crops thus far.

In the past two decades, consistent with this study, many studies have shown that gynoecy in bitter gourd is controlled by a single recessive gene (Behera et al., 2006; Ram et al., 2006; Behera et al., 2009), and this inheritance pattern is the same as that in other cucurbits – controlled by the WIP gene. Matsumura et al. (2014) first mapped the gynoecy locus, Mcgy, to a 12.73-cM genetic interval, which could be anchored to a 1.48-Mb physical interval from 19,622,042 bp to 21,103,779 bp on MC01 based on the Dali-11 genome (Cui et al., 2020). Later, Cui et al. (2018) also mapped the gynoecy locus, gy, to a 3.5-Mb physical interval between 17,619,724 bp and 21,144,642 bp on MC01. In this study, we identified the Mcgy1 locus responsible for gynoecy via BSA-seq (Fig. 2), and further delimited it into a 292.70-kb physical interval of 20,851,441 − 21,148,382 bp on MC01 by the use of classic forward genetic strategies, including molecular marker linkage analysis and fine-mapping (Fig. 3). Obviously, these three loci overlap, suggesting that the causal gene responsible for gynoecy is undisputedly located in this region, and the results of our study have greatly narrowed this region. We then performed a collinearity analysis involving the next-generation-sequenced Dali-11 genome and the long-read OHB3-1 genome to ensure that the genome assembly within this 292.70-kb region was correct (Supplementary Fig. S1). Notably, there is no homologous gene of WIP or CsACS1G in this 292.70-kb region, meaning that the causal gene of Mcgy1 in bitter gourd may be a novel gene different from that previously reported in other cucurbits. In addition, on the basis of a high-density genetic map with 2,013 SNP markers, a gynoecy locus, gy-1, was reported to be located in the middle of linkage group 12 (Gangadhara Rao et al., 2018); however, this locus cannot be anchored to the Dali-11 genome for comparison with the Mcgy1 locus due to the lack of sequence information on its flanking markers.

The causal gene responsible for gynoecy in the three gynoecious lines used in this study, namely, S156G, K20-95 and K20-506, is the same, since each was developed from a spontaneous gynoecious mutant, K44 (Cui et al., 2018). Hence, the four SNVs within the 292.70-kb fine-mapping region only belonging to the three gynoecious lines comparing with the nine monoecious inbred lines can be considered as candidate mutations leading to gynoecy (Fig. 4a). Three of these four SNVs, namely, SNV-1, SNV-3, and SNV-4, located at intergenic spacers and at least 2,432 bp away from their adjacent genes (Fig. 4b). Except for a few reports (Liu et al., 2017; Du et al., 2020), mutations occurring intergenic region that are more than 2 kb away from their adjacent genes usually have no effect on gene function. Meanwhile, combining expression analyses (Fig. 4c, d and Supplementary Figs S5 and S6), the most likely candidate mutation, SNV-2, was identified, which was further verified from genotyping 79 monoecious inbred lines (Fig. 4f and Supplementary Table S1). Notably, SNV-2 is located at the 3’ UTR of MC01g1681 based on gene annotation information (Cui et al., 2020). Some studies on the 3’ UTR have shown that this component can regulate mRNA-based processes, such as mRNA localization, mRNA stability, and mRNA translation (Macdonald & Struhl, 1988; An et al., 2008; Mayr, 2019), and the 3’ UTR can also store genetic information that can be transmitted to proteins through the formation of 3’ UTR-mediated protein‒protein interactions (Berkovits & Mayr, 2015). In addition, the 3’ UTRs are often highly conserved (Siepel et al., 2005; Xie et al., 2005). MC01g1681 encodes a CTPS, whose homologous normally play a role in the formation of filamentous structures in various organisms, such as Escherichia coli, Drosophila melanogaster, Saccharomyces cerevisiae and Homo sapiens (Liu, 2010; Barry et al., 2014; Noree et al., 2014; Lynch et al., 2017), and CTPS also plays rate-limiting roles at the final step of de novo synthesis of cytidine triphosphate (CTP), which is essential for DNA, RNA and phospholipid biosynthesis in all organisms (Daumann et al., 2018) and is vital for early endosperm development in Arabidopsis thaliana and Oryza sativa (Hickl et al., 2021; Yoon et al., 2021). These studies supported that SNV-2 can lead to changes in the gene function of MC01g1681, and our findings may be the first report that the 3’ UTR regulates sex differentiation.

Based on RNA-seq analysis, we found that only 21 DEGs showed significantly difference between S156G and S156 at all four developmental stages (Fig. 5a, b and Supplementary Table S26). Among these DEGs, MC02g0607 belongs to the AGL1 MADS-box family, which homologous involves in the regulation of gynoecium and ovule development in Arabidopsis (Flanagan et al., 1996); MC05g0014 is an uclacyanin-3-like gene, which homologous is related to pollen development in rice (Zhang et al., 2020b); and MC04g1310 is annotated as regulating anther development (Supplementary Table S25). The higher expression of MC02g0607 and the lower expression of MC04g1310 and MC05g0014 in S156G compared with S156 (Supplementary Table S26) may be part of an important molecular mechanism for the formation of gynoecy in bitter gourd. Five (MC06g1243, MC06g2033, MC09g1421, MC10g1173, and MC11g0366), one (MC03g0290), two (MC08g1661 and MC10g0031), and two (MC03g0811 and MC03g0812) DEGs, which are annotated as being involved in inorganic ion transport and metabolism, carbohydrate transport and metabolism, the formation of the integral component of membrane, and metal ion binding, respectively (Supplementary Table S25), implied that the molecular regulation for gynoecy may involve transmembrane transport. In addition, the annotations of three DEGs, namely, MC_newGene1691, MC06g1147, and MC10g1236, indicated that the molecular mechanism of the formation of gynoecy involves in DNA replication, RNA transport, and posttranslational modification (Supplementary Table S25). The functions of the remaining five genes are unknown. We believed that these 21 DEGs are vital clues to uncover the molecular mechanism of the formation of gynoecy, but more efforts are needed to elucidate their relationships and concrete mechanisms.

Previous studies have shown that the phytohormone ethylene plays a direct or indirect role in the formation of gynoecy in Cucurbitaceae crop species. For instance, by directly increasing the biosynthesis of ethylene, CsACS1G promotes gynoecy in cucumber (Mibus & Tatlioglu, 2004); this gene also promotes the femaleness in melon (Rudich et al., 1969). Epigenetic changes in CmWIP1 promote the expression of CmACS7 to indirectly enhance the biosynthesis of ethylene, leading to gynoecy in melon (Martin et al., 2009). In this study, we considered the two-leaf developmental stage to be the crucial stage for gynoecy formation on the basis of the expression patterns of MC01g1681 in stem apexes (Fig. 4c, d). Among this pathway, two DEGs related to ethylene signal transduction, MC11g0603, encoding a Constitutive Triple Response 1 (CTR1), and MC04g0109, encoding an Ethylene Insensitive 3 (EIN3), exhibited significantly higher and lower transcript levels in S156G than in S156, respectively (Supplementary Tables S25 and S27). Previous reports have shown that CTR1 plays a positive role in the biosynthesis of ethylene (Leclercq et al., 2002), while EIN has the opposite effect (Wang et al., 2020). These results suggested that ethylene may play an important role in the formation of gynoecy in bitter gourd. There were also three DEGs related to cytokinin signal transduction (Supplementary Tables S25 and S27), whose homologous play an important role in sex differentiation in watermelon (Zhang et al. 2020a). In addition, there are several other genes involved in plant hormone signal transduction, such as that of gibberellin and auxin (Supplementary Tables S25 and S27). Hence, ethylene may not be the only phytohormone influencing the formation of gynoecy in bitter gourd, but the metabolic mechanism by which Mcgy1 regulates sex differentiation needs to be clarified.

The Cucurbitaceae family comprises approximately 97 genera and 940 ~ 980 species (Schaefer & Renner, 2011), while only two pathways regulating gynoecy have been discovered in three species thus far. Here, we proposed a novel candidate gene responsible for gynoecy in bitter gourd. These results contribute to an improved understanding of the molecular formation mechanism of gynoecy in bitter gourd and enrich the knowledge about of genetic pathway of gynoecy in the Cucurbitaceae family.

Acknowledgements

This work was supported by the Key Project of Basic and Applied Research for University in Guangdong Province (2018KZDXM016), the Science and Technology Planning Project of Guangdong Province (2018B020202007), the Guangzhou Science and Technology Plan Projects (202002020086, 202102020800, and 202206010170), and the Science and Technology Plan Projects of Guangdong Province (2019A050520002).

Author contributions

KH and JWC conceived and designed the research. JZ and JJC performed most of the experiments and wrote the manuscript. JL and CFZ performed statistical analysis. FH and JCD provided helpful discussions. All authors read and approved the final manuscript.

Data availability

The data that support the findings of the current study are available from the corresponding author upon reasonable request.

Compliance with ethical standards

Conflict of interest The authors declare that they have no conflicts of interest.

Ethical statement The authors declare that this study complies with the current laws of China.

An JJ, Gharami K, Liao GY, Woo NH, Lau AG, Vanevski F, Torre ER, Jones KR, Feng Y, Lu B, Xu B (2008) Distinct role of long 3' UTR BDNF mRNA in spine morphology and synaptic plasticity in hippocampal neurons. Cell 134:175-187. https://doi.org/10.1016/j.cell.2008.05.045
Barrett SCH (2002) The evolution of plant sexual diversity. Nat Rev Genet 3:274-284. https://doi.org/10.1038/nrg776
Barry RM, Bitbol AF, Lorestani A, Charles EJ, Habrian CH, Hansen JM, Li HJ, Baldwin EP, Wingreen NS, Kollman JM, Gitai Z (2014) Large-scale filament formation inhibits the activity of CTP synthetase. Elife 3:e03638. https://doi.org/10.7554/eLife.03638
Basch E, Gabardi S, Ulbricht C (2003) Bitter melon (Momordica charantia): a review of efficacy and safety. Am J Health Syst Pharm 60:356-359. https://doi.org/10.1093/ajhp/60.4.356
Behera TK, Dey SS, Munshi AD, Gaikwad AB, Pal A, Singh I (2009) Sex inheritance and development of gynoecious hybrids in bitter gourd (Momordica charantia L.). Sci Hortic 120:130-133. https://doi.org/https://doi.org/10.1016/j.scienta.2008.09.006
Behera TK, Dey SS, Sirohi PS (2006) DBGy-201 and DBGy-202: Two gynoecious lines in bitter gourd (Momordica charantia L.) isolated from indigenous source. Indian J Genet Pl Br 66:61-62.
Berkovits BD, Mayr C (2015) Alternative 3' UTRs act as scaffolds to regulate membrane protein localization. Nature 522:363-367. https://doi.org/10.1038/nature14321
Boualem A, Fergany M, Fernandez R, Troadec C, Martin A, Morin H, Sari MA, Collin F, Flowers JM, Pitrat M, Purugganan MD, Dogimont C, Bendahmane A (2008) A conserved mutation in an ethylene biosynthesis enzyme leads to andromonoecy in melons. Science 321:836-838. https://doi.org/10.1126/science.1159023
Boualem A, Troadec C, Camps C, Lemhemdi A, Morin H, Sari MA, Fraenkel-Zagouri R, Kovalski I, Dogimont C, Perl-Treves R, Bendahmane A (2015) A cucurbit androecy gene reveals how unisexual flowers develop and dioecy emerges. Science 350:688-691. https://doi.org/10.1126/science.aac8370
Boualem A, Troadec C, Kovalski I, Sari MA, Perl-Treves R, Bendahmane A (2009) A conserved ethylene biosynthesis enzyme leads to andromonoecy in two cucumis species. PloS one 4:e6144. https://doi.org/10.1371/journal.pone.0006144
Byers RE, Baker LR, Sell HM, Herner RC, Dilley DR (1972) Ethylene: A Natural Regulator of Sex Expression of Cucumis melo L. Proc Natl Acad Sci U S A 69:717-720. https://doi.org/10.1073/pnas.69.3.717
Chen C, Chen H, Zhang Y, Thomas HR, Frank MH, He Y, Xia R (2020) TBtools: An Integrative Toolkit Developed for Interactive Analyses of Big Biological Data. Mol Plant 13:1194-1202. https://doi.org/10.1016/j.molp.2020.06.009
Chen H, Sun J, Li S, Cui Q, Zhang H, Xin F, Wang H, Lin T, Gao D, Wang S, Li X, Wang D, Zhang Z, Xu Z, Huang S (2016) An ACC Oxidase Gene Essential for Cucumber Carpel Development. Mol Plant 9:1315-1327. https://doi.org/10.1016/j.molp.2016.06.018
Chen Y, Chen Y, Shi C, Huang Z, Zhang Y, Li S, Li Y, Ye J, Yu C, Li Z, Zhang X, Wang J, Yang H, Fang L, Chen Q (2018) SOAPnuke: a MapReduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience 7:1-6. https://doi.org/10.1093/gigascience/gix120
Cui J, Luo S, Niu Y, Huang R, Wen Q, Su J, Miao N, He W, Dong Z, Cheng J, Hu K (2018) A RAD-Based Genetic Map for Anchoring Scaffold Sequences and Identifying QTLs in Bitter Gourd (Momordica charantia). Front Plant Sci 9:477. https://doi.org/10.3389/fpls.2018.00477
Cui J, Yang Y, Luo S, Wang L, Huang R, Wen Q, Han X, Miao N, Cheng J, Liu Z, Zhang C, Feng C, Zhu H, Su J, Wan X, Hu F, Niu Y, Zheng X, Yang Y, Shan D, Dong Z, He W, Dhillon NPS, Hu K (2020) Whole-genome sequencing provides insights into the genetic diversity and domestication of bitter gourd (Momordica spp.). Hortic Res 7:85. https://doi.org/10.1038/s41438-020-0305-5
Daumann M, Hickl D, Zimmer D, DeTar RA, Kunz HH, Möhlmann T (2018) Characterization of filament-forming CTP synthases from Arabidopsis thaliana. Plant J 96:316-328. https://doi.org/10.1111/tpj.14032
Delcher AL, Salzberg SL, Phillippy AM (2003) Using MUMmer to identify similar regions in large sequence sets. Curr Protoc Bioinformatics 00:1-18. https://doi.org/10.1002/0471250953.bi1003s00
Dellaporta SL, Calderon-Urrea A (1993) Sex determination in flowering plants. Plant Cell 5:1241-1251. https://doi.org/10.1105/tpc.5.10.1241
Flanagan CA, Hu Y, Ma H (1996) Specific expression of the AGL1 MADS-box gene suggests regulatory functions in Arabidopsis gynoecium and ovule development. Plant J 10:343-353. https://doi.org/10.1046/j.1365-313x.1996.10020343.x
Florea L, Song L, Salzberg SL (2013) Thousands of exon skipping events differentiate among splicing patterns in sixteen human tissues. F1000Res 2:188. https://doi.org/10.12688/f1000research.2-188.v2
Gangadhara Rao P, Behera TK, Gaikwad AB, Munshi AD, Jat GS, Boopalakrishnan G (2018) Mapping and QTL Analysis of Gynoecy and Earliness in Bitter Gourd (Momordica charantia L.) Using Genotyping-by-Sequencing (GBS) Technology. Front Plant Sci 9:1555. https://doi.org/10.3389/fpls.2018.01555
Hamrick JL, Godt MJW (1996) Effects of life history traits on genetic diversity in plant species. PHILOS T R SOC B 351:1291-1298. https://doi.org/doi:10.1098/rstb.1996.0112
Hickl D, Scheuring D, Möhlmann T (2021) CTP Synthase 2 From Arabidopsis thaliana Is Required for Complete Embryo Development. Front Plant Sci 12:652434. https://doi.org/10.3389/fpls.2021.652434
Hill JT, Demarest BL, Bisgrove BW, Gorsi B, Su YC, Yost HJ (2013) MMAPPR: mutation mapping analysis pipeline for pooled RNA-seq. Genome Res 23:687-697. https://doi.org/10.1101/gr.146936.112
Hu B, Li D, Liu X, Qi J, Gao D, Zhao S, Huang S, Sun J, Yang L (2017) Engineering Non-transgenic Gynoecious Cucumber Using an Improved Transformation Protocol and Optimized CRISPR/Cas9 System. Mol Plant 10:1575-1578. https://doi.org/10.1016/j.molp.2017.09.005
Ji G, Zhang J, Zhang H, Sun H, Gong G, Shi J, Tian S, Guo S, Ren Y, Shen H, Gao J, Xu Y (2016) Mutation in the gene encoding 1-aminocyclopropane-1-carboxylate synthase 4 (CitACS4) led to andromonoecy in watermelon. J Integr Plant Biol 58:762-765. https://doi.org/10.1111/jipb.12466
Khan A, Fornes O, Stigliani A, Gheorghe M, Castro-Mondragon JA, van der Lee R, Bessy A, Chèneby J, Kulkarni SR, Tan G, Baranasic D, Arenillas DJ, Sandelin A, Vandepoele K, Lenhard B, Ballester B, Wasserman WW, Parcy F, Mathelier A (2018) JASPAR 2018: update of the open-access database of transcription factor binding profiles and its web framework. Nucleic Acids Res 46:D260-D266. https://doi.org/10.1093/nar/gkx1126
Kim D, Langmead B, Salzberg SL (2015) HISAT: a fast spliced aligner with low memory requirements. Nat Methods 12:357-360. https://doi.org/10.1038/nmeth.3317
Kole C (2020) The Bitter Gourd Genome. Springer, Cham Switzerland. https://doi.org/10.1007/978-3-030-15062-4
Leclercq J, Adams-Phillips LC, Zegzouti H, Jones B, Latché A, Giovannoni JJ, Pech JC, Bouzayen M (2002) LeCTR1, a tomato CTR1-like gene, demonstrates ethylene signaling ability in Arabidopsis and novel expression patterns in tomato. Plant Physiol 130:1132-1142. https://doi.org/10.1104/pp.009415
LEWIS D (1942) The evolution of sex in flowering plants. Biol Rev 17:46-67. https://doi.org/https://doi.org/10.1111/j.1469-185X.1942.tb00431.x
Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25:1754-1760. https://doi.org/10.1093/bioinformatics/btp324
Li Z, Han Y, Niu H, Wang Y, Jiang B, Weng Y (2020) Gynoecy instability in cucumber (Cucumis sativus L.) is due to unequal crossover at the copy number variation-dependent Femaleness (F) locus. Hortic Res 7:32. https://doi.org/10.1038/s41438-020-0251-2
Li Z, Huang S, Liu S, Pan J, Zhang Z, Tao Q, Shi Q, Jia Z, Zhang W, Chen H, Si L, Zhu L, Cai R (2009) Molecular isolation of the M gene suggests that a conserved-residue conversion induces the formation of bisexual flowers in cucumber plants. Genetics 182:1381-1385. https://doi.org/10.1534/genetics.109.104737
Liu JL (2010) Intracellular compartmentation of CTP synthase in Drosophila. J Genet Genomics 37:281-296. https://doi.org/10.1016/s1673-8527(09)60046-1
Livak KJ, Schmittgen TD (2001) Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 25:402-408. https://doi.org/10.1006/meth.2001.1262
Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15:550. https://doi.org/10.1186/s13059-014-0550-8
Lynch EM, Hicks DR, Shepherd M, Endrizzi JA, Maker A, Hansen JM, Barry RM, Gitai Z, Baldwin EP, Kollman JM (2017) Human CTP synthase filament structure reveals the active enzyme conformation. Nat Struct Mol Biol 24:507-514. https://doi.org/10.1038/nsmb.3407
Macdonald PM, Struhl G (1988) Cis-acting sequences responsible for anterior localization of bicoid mRNA in Drosophila embryos. Nature 336:595-598. https://doi.org/10.1038/336595a0
Manzano S, Aguado E, Martínez C, Megías Z, García A, Jamilena M (2016) The Ethylene Biosynthesis Gene CitACS4 Regulates Monoecy/Andromonoecy in Watermelon (Citrullus lanatus). PloS one 11:e0154362. https://doi.org/10.1371/journal.pone.0154362
Manzano S, Martínez C, García JM, Megías Z, Jamilena M (2014) Involvement of ethylene in sex expression and female flower development in watermelon (Citrullus lanatus). Plant Physiol Biochem 85:96-104. https://doi.org/10.1016/j.plaphy.2014.11.004
Manzano S, Martínez C, Megías Z, Gómez P, Garrido D, Jamilena M (2011) The role of ethylene and brassinosteroids in the control of sex expression and flower development in Cucurbita pepo. Plant Growth Regul 65:213-221. https://doi.org/10.1007/s10725-011-9589-7
Mao X, Cai T, Olyarchuk JG, Wei L (2005) Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 21:3787-3793. https://doi.org/10.1093/bioinformatics/bti430
Martin A, Troadec C, Boualem A, Rajab M, Fernandez R, Morin H, Pitrat M, Dogimont C, Bendahmane A (2009) A transposon-induced epigenetic change leads to sex determination in melon. Nature 461:1135-1138. https://doi.org/10.1038/nature08498
Matsumura H, Hsiao MC, Lin YP, Toyoda A, Taniai N, Tarora K, Urasaki N, Anand SS, Dhillon NPS, Schafleitner R, Lee CR (2020) Long-read bitter gourd (Momordica charantia) genome and the genomic architecture of nonclassic domestication. Proc Natl Acad Sci U S A 117:14543-14551. https://doi.org/10.1073/pnas.1921016117
Matsumura H, Miyagi N, Taniai N, Fukushima M, Tarora K, Shudo A, Urasaki N (2014) Mapping of the gynoecy in bitter gourd (Momordica charantia) using RAD-seq analysis. PloS one 9:e87138. https://doi.org/10.1371/journal.pone.0087138
Mayr C (2019) What Are 3' UTRs Doing? Cold Spring Harb Perspect Biol 11:a034728. https://doi.org/10.1101/cshperspect.a034728
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA (2010) The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res 20:1297-1303. https://doi.org/10.1101/gr.107524.110
Mibus H, Tatlioglu T (2004) Molecular characterization and isolation of the F/f gene for femaleness in cucumber ( Cucumis sativus L.). Theor Appl Genet 109:1669-1676. https://doi.org/10.1007/s00122-004-1793-7
Noree C, Monfort E, Shiau AK, Wilhelm JE (2014) Common regulatory control of CTP synthase enzyme activity and filament formation. Mol Biol Cell 25:2282-2290. https://doi.org/10.1091/mbc.E14-04-0912
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL (2015) StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol 33:290-295. https://doi.org/10.1038/nbt.3122
Porebski S, Bailey LG, Baum BR (1997) Modification of a CTAB DNA extraction protocol for plants containing high polysaccharide and polyphenol components. Plant Mol Biol Rep 15:8-15. https://doi.org/10.1007/BF02772108
Ram D, Kumar S, Singh M, Rai M, Kalloo G (2006) Inheritance of gynoecism in bitter gourd (Momordica charantia L.). J Hered 97:294-295. https://doi.org/10.1093/jhered/esj028
Rudich J, Halevy AH, Kedar N (1969) Increase in femaleness of three cucurbits by treatment with Ethrel, an ethylene releasing compound. Planta 86:69-76. https://doi.org/10.1007/bf00385305
Schaefer H, Renner SS (2010) A three-genome phylogeny of Momordica (Cucurbitaceae) suggests seven returns from dioecy to monoecy and recent long-distance dispersal to Asia. Mol Phylogenet Evol 54:553-560. https://doi.org/10.1016/j.ympev.2009.08.006
Schaefer H, Renner SS (2011) Cucurbitaceae. In: Kubitzki K (ed) Flowering Plants Eudicots: Sapindales, Cucurbitales, Myrtaceae. Springer, Berlin Heidelberg, pp 112-174. https://doi.org/10.1007/978-3-642-14397-7_10
Siepel A, Bejerano G, Pedersen JS, Hinrichs AS, Hou M, Rosenbloom K, Clawson H, Spieth J, Hillier LW, Richards S, Weinstock GM, Wilson RK, Gibbs RA, Kent WJ, Miller W, Haussler D (2005) Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomes. Genome Res 15:1034-1050. https://doi.org/10.1101/gr.3715005
Tan J, Tao Q, Niu H, Zhang Z, Li D, Gong Z, Weng Y, Li Z (2015) A novel allele of monoecious (m) locus is responsible for elongated fruit shape and perfect flowers in cucumber (Cucumis sativus L.). Theor Appl Genet 128:2483-2493. https://doi.org/10.1007/s00122-015-2603-0
Van Ooijen JW (2006) Joinmap4, software for calculation of genetic linkage maps in experimental populations. Kyazma BV, Wageningen.
Wang K, Li M, Hakonarson H (2010) ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res 38:e164. https://doi.org/10.1093/nar/gkq603
Wang L, Ko EE, Tran J, Qiao H (2020) TREE1-EIN3-mediated transcriptional repression inhibits shoot growth in response to ethylene. Proc Natl Acad Sci U S A 117:29178-29189. https://doi.org/10.1073/pnas.2018735117
Wang S, Li Z, Yang G, Ho CT, Li S (2017) Momordica charantia: a popular health-promoting vegetable with multifunctionality. Food Funct 8:1749-1762. https://doi.org/10.1039/c6fo01812b
Xie X, Lu J, Kulbokas EJ, Golub TR, Mootha V, Lindblad-Toh K, Lander ES, Kellis M (2005) Systematic discovery of regulatory motifs in human promoters and 3' UTRs by comparison of several mammals. Nature 434:338-345. https://doi.org/10.1038/nature03441
Yoon J, Cho LH, Kim SR, Tun W, Peng X, Pasriga R, Moon S, Hong WJ, Ji H, Jung KH, Jeon JS, An G (2021) CTP synthase is essential for early endosperm development by regulating nuclei spacing. Plant Biotechnol J 19:2177-2191. https://doi.org/10.1111/pbi.13644
Zhang J, Guo S, Ji G, Zhao H, Sun H, Ren Y, Tian S, Li M, Gong G, Zhang H, Xu Y (2020a) A unique chromosome translocation disrupting ClWIP1 leads to gynoecy in watermelon. Plant J 101:265-277. https://doi.org/10.1111/tpj.14537
Zhang J, Shi J, Ji G, Zhang H, Gong G, Guo S, Ren Y, Fan J, Tian S, Xu Y (2017) Modulation of Sex Expression in Four Forms of Watermelon by Gibberellin, Ethephone and Silver Nitrate. Hortic Plant J 3:91-100. https://doi.org/https://doi.org/10.1016/j.hpj.2017.07.010
Zhang YC, He RR, Lian JP, Zhou YF, Zhang F, Li QF, Yu Y, Feng YZ, Yang YW, Lei MQ, He H, Zhang Z, Chen YQ (2020b) OsmiR528 regulates rice-pollen intine formation by targeting an uclacyanin to influence flavonoid metabolism. Proc Natl Acad Sci U S A 117:727-732. https://doi.org/10.1073/pnas.1810968117

Download PDF

Journal Publication

published 23 Mar, 2023

Read the published version in Theoretical and Applied Genetics →

Editorial decision: Major revisions
11 Nov, 2022
Reviewers agreed at journal
09 Oct, 2022
Editor assigned by journal
26 Sep, 2022
First submitted to journal
25 Sep, 2022

You are reading this latest preprint version

Fine-mapping and candidate gene analysis of the Mcgy1 locus responsible for gynoecy in bitter gourd (Momordica spp.)

Status:

Journal Publication

Version 1

Abstract

Figures

Key Message

Introduction

Materials And Methods

Plant materials

Whole-genome resequencing and variant calling

BSA-seq

Molecular marker linkage analysis and fine-mapping

Multiple genomic sequence variation analysis

RNA-seq and comparative transcriptome analysis

Results

Characteristics and inheritance of gynoecy

BSA-seq to rapidly identify the Mcgy1 locus

Classic molecular marker linkage analysis validates the Mcgy1 locus

Fine-mapping of the Mcgy1 locus

RNA-seq of stem apexes provides a preliminary estimate for the candidate genes of Mcgy1

Multiple genomic sequence variation analysis to define the candidate genes of Mcgy1

Expression analysis suggests that MC01g1681 is an important candidate for Mcgy1

MC01g1681 can be verified in different bitter gourd germplasms using molecular marker

Transcriptome analysis reveals the potential molecular formation mechanisms of gynoecy

Discussion

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1