Transcriptomics reveal useful resources for examining fruit development and variation in fruit size in Coccinia grandis

Introduction The Cucurbitaceae family comprises many agronomically important members, that bear nutritious fruits and vegetables of great economic importance. Coccinia grandis, commonly known as Ivy gourd, belongs to this family and is widely consumed as a vegetable. Members of this family are known to display an impressive range of variation in fruit morphology. Although there have been studies on flower development in Ivy gourd, fruit development remains unexplored in this crop. Methods In this study, comparative transcriptomics of two Ivy gourd cultivars namely “Arka Neelachal Kunkhi” (larger fruit size) and “Arka Neelachal Sabuja” (smaller fruit size) differing in their average fruit size was performed. A de novo transcriptome assembly for Ivy gourd was developed by collecting fruits at different stages of development (5, 10, 15, and 20 days after anthesis i.e. DAA) from these two varieties. The transcriptome was analyzed to identify differentially expressed genes, transcription factors, and molecular markers. Results The transcriptome of Ivy gourd consisted of 155205 unigenes having an average contig size of 1472bp. Unigenes were annotated on publicly available databases to categorize them into different biological functions. Out of these, 7635 unigenes were classified into 38 transcription factor (TF) families, of which Trihelix TFs were most abundant. A total of 11,165 unigenes were found to be differentially expressed in both the varieties and the in silico expression results were validated through real-time PCR. Also, 98768 simple sequence repeats (SSRs) were identified in the transcriptome of Ivy gourd. Discussion This study has identified a number of genes, including transcription factors, that could play a crucial role in the determination of fruit shape and size in Ivy gourd. The presence of polymorphic SSRs indicated a possibility for marker-assisted selection for crop breeding in Ivy gourd. The information obtained can help select candidate genes that may be implicated in regulating fruit development and size in other fruit crops.


Introduction
Fruit is a characteristic feature of angiosperms and can be a chief source of vitamins, fibres, carbohydrates and other nutrients necessary for human diet (Zhao et al., 2021).Fruits in flowering plants usually develop from matured ovary or in some cases, accessory floral tissues contribute in fruit formation (Seymour et al., 2011;Zhao et al., 2021).In addition to being the edible part of most fruit-bearing plants, they also help in seed dispersal after maturation (Azzi et al., 2015), thereby allowing the successful adaptation of angiosperms on Earth.Fruit has diverse morphology in order to cope with changing environments (Wang et al., 2022) and usually undergoes three distinct stages of development; fruit set, fruit growth, and fruit maturation (Azzi et al., 2015).In the first stage of development, the ovary's growth is followed by successful fertilization, which leads to the decision for fruit setting.In the second stage of fruit development, fruit grows due to cell division.The formation of seed and embryo development takes place at this stage.Hereafter, the fruit cells don't divide, rather the fruit increases its volume by cell expansion.This stage is the longest and decides the size of a fruit.The third stage is known as fruit maturation where deposition of different types of storage products in fruits occurs and at this stage, fruit size increases very slowly (Gillaspy et al., 1993;Mu et al., 2009;Zhao et al., 2021).
Fruit size and shape constitute important agronomic traits and also contribute to plant fitness in the course of evolution (Hussain et al., 2020).Fruits that are larger in size and weight are usually preferred by consumers, lending them greater economic value (Huang et al., 2023).Many external and internal factors influence fruit size including genetics, environment, cultivation practices etc. and have been studied in various reports (Huang et al., 2023;Jahed and Hirst, 2023).Apart from these, molecular mechanisms such as cell division, cell expansion, and endoreduplication, play a vital role in determining fruit size and shape (Chevalier et al., 2014;Okello et al., 2015;Mauxion et al., 2021).More detailed analyses have revealed the role of phytohormones, transcription factors, quantitative trail loci (QTLs), miRNAs, and signalling pathways in fruit development (Guo et al., 2010;Libault et al., 2010;Xu et al., 2013;Hussain et al., 2020;Huang et al., 2023).Nowadays, developing plants with gene knock-out and knock-in is feasible due to advanced technologies for genome editing and plant transformation protocols.Hence the only challenge is to identify genes and regulatory elements controlling fruit sizes that can be later implemented for crop improvement through genome editing and through marker-assisted selection (Hussain et al., 2020).
Ivy gourd (Coccinia grandis) belongs to the Cucurbitaceae family and is a nutrient-rich vegetable found in sub-Saharan Africa, tropical Asia, neo-tropics, and India (Mohanty et al., 2017).The cucurbitaceae family consists of many members that produce a variety of fruits and vegetables for everyday consumption.Great variations have been observed for fruit sizes in this family which comprises of many important vegetables such as pumpkin and bottle gourd (Xanthopoulou et al., 2017;Zhang et al., 2020).Such morphological differences in fruits of the same lineages can be used to understand the genetic entity that controls the fruit size.We have selected two varieties of Ivy gourd namely "Arka Neelachal Kunkhi" and "Arka Neelachal Sabuja", which have considerable differences in fruit size and shape, with "Arka Neelachal Kunkhi" being elongated in shape and "Arka Neelachal Sabuja" being ovoid shaped."Arka Neelachal Kunkhi" is a high yielding variety that produces around 800 fruits per season with yield potential of 15-20 tonnes/ha.It is a striped Ivy gourd consumed both in raw (in salads) as well as in cooked form.On the other hand, "Arka Neelachal Sabuja" is a vigorously growing plant that yields 70-80 harvests per season with a yield potential of 20-25t/ha (https://www.iihr.res.in/central-horticultural-experiment-stationches-bhubaneswar).The fruits are consumed widely in India and are a good source of several vital minerals and vitamins.Therefore, analyzing the factors associated with fruit morphology at a molecular level will be helpful in future gene-editing efforts.
Various genes regulating phytohormones such as auxin, gibberellins, cytokines, abscisic acid and ethylene are known to regulate fruit development in plants.For example, over expression of GmYUCCA5 gene resulted in short siliques in Arabidopsis (Wang et al., 2017).Various transcription factor families such as B3 family, YABBY family, Zinc finger family are known to regulate fruit sizes.In addition to that, microtubule play an active role in cell division and hence play a role in fruit size determination (Hussain et al., 2020).With the progress in high-throughput sequencing techniques with computational, statistical methods and development of multiomics, comparative transcriptomics has become a norm for identification of key genes underlying important biological phenomena such as control of fruit growth and size (Peng et al., 2022;Jahed and Hirst, 2023).Transcriptomic analysis can reveal genes regulating important agronomic traits in plants and can also reveal differences in gene expression in response to environmental factors (Zhang et al., 2020).The method has been used in numerous studies where the researchers have identified genes based on their differential expression in various tissues and developmental stages of fruits in the cucurbits (Xanthopoulou et al., 2017;Abbas et al., 2020;Zhang et al., 2020;Luo et al., 2021;Xanthopoulou et al., 2021).The present study aims to add to the limited collection of genes and regulatory elements that govern fruit development and morphology in cucurbits.Here we have generated the transcriptome of two cultivars of Ivy gourd ("Arka Neelachal Kunkhi" and "Arka Neelachal Sabuja") with different fruit sizes and to identify genes, regulators, and key pathways involved in fruit development and size regulation.Information obtained will be very useful for rearchitecting the shape and size of fruit using gene-based technology.DAA and 20 DAA, for both the varieties of C. grandis.Leaf tissue was also collected from both plants to use for the generation of transcriptome to generate better coverage for the transcriptome which will allow better assembly.All samples were collected in triplicates (biological), snap-frozen in liquid N2, and stored at -80°C till further processing.Fresh fruits were used for morphological studies including measurements of fruit length and width.Anatomical studies were also performed through the observations of thin cross sections from different stages of fruit development in both the varieties under light microscope (Model Stemi508, Make Zeiss) and Scanning electron microscope (SEM -FESEM-EDX, Model JSM-IT800, Make JEOL).The SEM images (Scale 10 µm) were exported to the ImageJ software for cell number count and area measurement in µm 2 .These steps were performed with three biological replicates.

RNA isolation and RNA-Seq
Total RNA was extracted from 100mg of frozen plant tissues using a Nucleospin Plant and Fungi RNA extraction kit (REF 740120.50,Düren,Germany)following the manufacturer's protocol.The quantity and quality of the RNA were analyzed on 1% agarose gel, and Qubit 4.0 fluorometer (Thermo Fisher Scientific, USA) and an Agilent TapeStation system (Agilent, USA).A total of 1mg of total RNA with RNA integrity number (RIN) > 7 was used for library preparation.The libraries were prepared taking three biological replicates of each sample (fruits of 5 DAA, 10 DAA, 15 DAA, 20 DAA, young leaves) and using the TruSeq Stranded mRNA Library Prep kit (Illumina, San Diego, USA) following the manufacturer's instructions.Library quality was evaluated using (Thermo Fisher Scientific, USA) and Agilent TapeStation system (Agilent, USA).After quality assessment, a total of thirty libraries representing three biological replicates each of five samples of "Arka Neelachal Kunkhi" (hereafter depicted as "ANKu") and "Arka Neelachal Sabuja" (hereafter depicted as "ANSa") were pooled and paired-end sequenced on the NextSeq550 (2x150) platform (Illumina, USA) at Institute of Life Sciences, Bhubaneswar, Odisha.

Processing of reads and transcriptome assembly
A total of 405 million reads in "ANKu" and 566 million reads in "ANSa" were obtained from all the libraries after demultiplexing with bcl2fastq software (Illumina, USA).The raw reads have been submitted to NCBI SRA database with the BioProject ID PRJNA851184.The reads were screened for quality and trimmed with Trimmomatic v. 0.39 with the following parameters: (ILLUMINACLIP: TruSeq3-PE2.fa:2:30:10LEADING:5 TRAILING:5 SLIDINGWINDOW:5:10 MINLEN:50).The highquality reads after removal of adapter sequences were used for de novo assembly using the Trinity pipeline (version 2.11) with the default parameters (Haas et al., 2013).CD-HIT-EST (Li and Godzik, 2006) was performed at 98% similarity to remove redundant contigs.

Assembly quality assessment and annotation
The quality of the de novo transcriptome assembly was assessed in two ways (i) by calculating the percentage of reads mapped onto the assembled unigenes and (ii) by using the curated datasets at BUSCO (Simão et al., 2015).For the first method, the filtered reads from each sample were mapped onto the assembled unigenes using Bowtie2 (Langmead and Salzberg, 2012) and percentage of mapped reads were plotted in a graph on MS Excel.The BUSCO analysis was performed using the eudicot_odb10 data set (Simão et al., 2015).

Transcription factor family identification and SSR mining
Transcription factors were identified in the Ivy gourd transcriptome in two stages.In the first pass, the HMM profiles of all available gene families were downloaded from Pfam (http:// ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz) and used to search for the conserved Hidden Markov Model (HMM) profiles in the peptide sequences of Ivy gourd.However, the HMM profiles for all the transcription factor families were not present in the Pfam database.So, in a second pass of TF prediction, peptide sequences for individual TF families of Arabidopsis were downloaded from TAIR and aligned to generate HMM profiles using HMMER v. 3.3 (http://hmmer.org/).These profiles were then used to search for TF families in the peptide sequences of Ivy gourd.
Simple sequence repeats (SSRs) were identified in the unigenes using MISA software (Beier et al., 2017) with the following parameters: Mono-nucleotides 10 repeats, di-nucleotide 8, triand tetra-nucleotide 4 repeats, penta-nucleotide and hexanucleotide 3 repeats.A total of 20 SSR markers associated with differentially expressed genes were randomly selected having a length of 15bp to 24bp.Primers for genic SSRs sequences were designed by using MISA web tool (https://webblast.ipkgatersleben.de/misa/index.php?action=1) based on the following criteria; a G/C content between 40% and 70%, an annealing temperature between 54°C to 63°C, minimum product length of 100bp and primer length 10-24 nucleotides.All the candidate SSR primer pairs were synthesized by Eurofins Scientific India Pvt Ltd (Bengaluru, India).The PCR reactions were conducted with 25ul reaction volume, and following conditions: initial denaturation at 95°C for 3min followed by 30 cycles of denaturation at 95°C for 30sec, annealing at 54°C for 45sec and extension at 72°C for 45 sec.The final extension was done at 72°C for 5min and hold at 4°C.PCR amplifications were performed using genomic DNA isolated from 15 DAA stage of fruit (using NucleoSpin ™ Plant II DNA extraction kit, Macherey-Nagel) as template and the PCR products were verified with 3% metaphor agarose (Lonza, Basel, Switzerland) and 100bp DNA ladder (GeNei ™ , Bangalore, India).

Differential gene expression analysis, GO term enrichment and biological pathway analysis
In silico gene expression analysis was carried out as described by (Suranjika et al., 2022).Briefly, we used Bowtie2 (Langmead and Salzberg, 2012) to map the short reads from all sample libraries (including replicates), onto the assembled transcriptome and abundance was estimated by RSEM i.e RNA-Seq by Expectation-Maximization (Parrish et al., 2011).Empirical Analysis of Digital Gene Expression (EdgeR) was used to obtain differentially expressed unigenes in Ivy gourd as well as to normalize the expected counts for relative expression (Robinson et al., 2010) with FDR correction ≤ 0.05, P value ≤ 0.001, fold change ≥ 2. The normalization factor of effective library size was performed by Trimmed Mean of M-values (TMM).Finally, the unigenes with threshold FDR ≤ 0.05 and log fold change |log 2 FC| (≥ 1) were considered to be differentially expressed and selected for further analysis.

Quantitative real time PCR analysis
In order to check the expression of few genes to relate their roles in fruit size regulation, qRT-PCR analysis was performed by using C. grandis 18S rRNA as endogenous control to normalize the gene expression.Primers were designed by using the Primer ™ Quest tool of Integrated DNA technology (IDT) software (Supplementary Table S1).Total RNA of all the stages of developing fruits (5 DAA, 10 DAA, 15 DAA, and 20 DAA) and with three biological replicates was isolated.First strand cDNA synthesis was accomplished by using 1mg of freshly isolated RNA using First strand cDNA synthesis kit (ThermoFisher Scientific, Waltham, MA, USA) following the specifications of the manufacturer.The cDNA products were diluted five times (1:5) and stored in -20°C.For Real Time PCR, each sample reaction was set up in a reaction mixture of 10ml containing 5ml of 2X GoTaq qPCR Master mix (Promega Corporation, Madison, WI, USA), 1 ml of gene-specific primer mix, 0.5 ml of 1:5 diluted cDNA and nuclease-free water for volume adjustments (Promega Corporation, Madison, WI, USA).The qRT-PCR assay was performed using the QuantStudio5 system (Applied Biosystems, ThermoFisher Scientific, Waltham, MA, USA), with the following reaction conditions; 95°C for 10min (initial denaturation) followed by 35 cycles of 95°C for 15s and 55°C for 1min.A melt curve program was performed by initiating at 95°C for 15s and reducing the temperature by 1.6°C/s until it reaches 60°C for 1 min followed by an increase in 0.15°C/s until it reaches 95°C.The relative fold changes were calculated by using the comparative 2 -DDCt method (Livak and Schmittgen, 2001).The statistical analysis was done by calculating the standard deviation and the level of significance was determined by the Student's t-test (P ≤ 0.05).

Morphological and anatomical characterization of Ivy gourd var "Arka
Neelachal Kunkhi" and "Arka Neelachal Sabuja" Phenotypic observations of two varieties "Arka Neelachal Kunkhi" (ANKu) and "Arka Neelachal Sabuja" (ANSa) showed remarkable size and shape differences in all the developmental stages taken in the present study (Figure 1).Two fruit variables such as length and width from all the collected samples (5 DAA, 10 DAA, 15 DAA, 20 DAA) were measured.At maturity, the mean length of ANKu (20 DAA) was found to be 9.02 cm and that of ANSa was 6.7 cm.Likewise, the diameter or mean width of ANKu and ANSa were measured as 7 cm and 8.92 cm respectively (Supplementary Table S2).Anatomical studies of fruits at each stage showed that the locule of ANSa expands gradually in every growing stage and the number of seeds in mature ANSa fruits is higher in comparison with ANKu (Figure 2).More detailed studies using SEM also showed differences in cell number in different stages of fruit development in both varieties (Figure 3).More number of cells were observed in ANKu when compared to ANsa.The cell numbers were observed to increase rapidly during the early stages of fruit development (5 DAA, 10 DAA) in both the varieties of C. grandis.However, this rapid increase stagnated as the fruits reached maturity and no rapid increase in number was observed (Supplementary Figure S1A).The number of cells were found to be higher in ANKu than in ANSa.To examine the cell expansion in both varieties, the cell area of different developmental stages was measured.Interestingly, we observed, a larger cell area in ANSa than in ANKu (Supplementary Figure S1B).Moreover, maximum expansion of cells was observed in 10 DAA to 20 DAA fruits of ANSa.

Assembly and annotation of the Ivy gourd transcriptome
A total of 971 million reads were obtained after short read sequencing of the two varieties of Ivy gourd.These reads were filtered to remove low quality reads and adapter contamination to obtain 948 million reads, giving us 96% retention of high-quality reads on average (Supplementary Table S3).A total of 155205 unigenes remained after processing, with an N50 value of 1472 bp (Table 1).During quality assessment, we observed that on average, 97.6% of the RNA-seq reads could be mapped onto the assembly (Supplementary Figure S2A).Analysis using the conserved orthologs in BUSCO database showed that 95.4% of complete BUSCOs from eudicots_odb10 database were represented in the Ivy gourd transcriptome (Supplementary Figure S2B).
We annotated the transcriptome of Ivy gourd against various databases and found that only 55,377 unigenes could be annotated, which means that only ~36% of the unigenes could find a match in any of the databases (Figure 4A).Using the results of Blastx against the Uniprot-Swissprot database, we assigned GO terms to the unigenes and found categories like "metabolic process" and "cellular process" were over-represented under the term "Biological Process" while categories like "catalytic activity" and "binding" accounted for the majority of unigenes under the "Molecular Function term" (Figures 4B-D).Transdecoder tool (https://github.com/TransDecoder/TransDecoder)predicted 47,296 complete CDS in the transcriptome of Ivy gourd which accounts for 30% of the total unigenes.

Differential gene expression analysis
We generated expression profiles for the entire transcriptome after comparing the fruits from the two varieties at individual developmental stages i.e at 5, 10, 15 and 20 DAA (Figures 5A-D; Supplementary Figure S3).We found that at 5 DAA, compared with ANKu, fruits of ANSa showed up-regulation of genes encoding U- regulated in ANSa with most of them being upregulated while genes like Cytochrome P450, MYB TFs were downregulated (Supplementary Table S5).A comparison of the DEUs at each stage of fruit development revealed 433 unigenes that were common to all stages (Figure 6).These comprised of genes coding for MYB TF, heat stress transcription factor, trihelix TF, UDPglycosyltransferase etc.
In addition to analyzing the unigenes obtained in this study, we also investigated the expression patterns for genes previously reported in fruit development and morphology in cucurbits (Dou et al., 2018;Pan et al., 2020;Boualem et al., 2022) We extracted the sequences for these genes using the information provided in their respective publications and identified their homologs in the C. grandis transcriptome using BLASTN and subjected them to insilico gene expression analyses (Supplementary Table S4).We observed that most of the genes related to ethylene biosynthesis such as homologs of CpETR1, CpERS1, CpCTR2 and CpEIN3.2 were downregulated in ANKu as compared with ANSa.It was interesting to note that homologs of the genes Cla011257 was downregulated in all stages of development in AnKu, except at 15 DAA while another homolog of CpCTR1 was upregulated at 10, 15 and 20 DAA in ANKu as compared with ANSa (Figure 7).It was also curious to note that few of the crucial genes like CpACS7 did not find a homolog in C. grandis transcriptome.

Gene enrichment at different stage of fruit development
Taking ANKu as a reference, we analyzed the DEUs that were enriched in each stage of fruit development of the varieties of Ivy gourd to look for molecular functions that could be involved in regulating fruit morphology in ANKu and ANSa.The results showed that various binding activities like "nucleotide binding", "ATP binding" and catalytic activities like "hydrolase" and "transferase" activity were differentially regulated at 5 DAA.These activities were also found to be enriched at 10 DAA with the addition of "transporter" and "pyrophosphatase" activity.The later stages of 15 and 20 DAA showed a drastic change in the gene enrichment pattern with a clear reduction in the representation of catalytic activity like "transporter" (Supplementary Figures S4A-D).

Transcription factors regulating fruit development and morphology in C. grandis
We identified a total of 7,635 unigenes coding for some of the most widely studied transcription factors in the transcriptome of C.   8).The expression profiles for these three families of TFs in the two varieties of C. grandis showed that a number of Trihelix TFs were differentially expressed in ANSa as compared to ANKu.Some of these unigenes such as TRINITY_DN8473_c0_g1_i1, TRINITY_DN4341_c0_g1_i9, TRINITY_DN14419_c1_g2_i5, TRINITY_DN94151_c0_g1_i1, TRINITY_DN8473_c0_g1_i1 and TRINITY_DN6940_c0_g1_i2 were upregulated in the ANSa variety (Figure 9A).On the other hand, the unigenes TRINITY_DN6318_c0_g1_i6, TRINITY_DN7952_c0_g1_i7, TRINITY_DN8318_c1_g2_i1 and TRINITY_DN5103_c0_g1_i7 were downregulated in ANSa (Figure 9A).Such clear preferential expression was not observed for the MYB and zf-C3HC4 TFs (Figures 9B, C).However, we observed that the unigenes TRINITY_DN 4808_c0_g1_i3 and TRINITY_DN8591_c0_g1_i3, coding for MYB TFs were downregulated in the ANSa variety while TRINITY_DN7452_c0_g2_i13 was upregulated in ANSa as compared with ANKu (Figure 9B).

SSRs identified in C. grandis transcriptome
We searched for simple sequence repeats (SSRs) in the unigenes of Ivy gourd and identified 98768 SSRs in 155205 unigenes.Majority of these SSRs were mononucleotide repeats (64914), followed by trinucleotide repeats (18702) (Table 2).Since we have identified these SSRs in the transcriptome of Ivy gourd, these would be regarded as genic SSRs and have a high possibility of being associated with protein coding genes.Hence, we also analyzed the SSRs in the differentially expressed subset of unigenes and found that a number of these SSRs were present in the differentially expressed unigenes which were annotated to code for various transcription factors, genes involved in biosynthesis of metabolites, ubiquitin mediated proteolysis, structural elements and signalling pathways (Supplementary Table S6).
We designed primers for 20 SSRs associated with DEUs and were able to successfully amplify 16 through PCR in three biological replicates of 15DAA fruit tissue (Supplementary Figure S5).These SSRs were associated with unigenes coding for Leucine-rich repeat extensin-like protein, Trihelix transcription factor and 3hydroxybutyryl-CoA dehydrogenase (Supplementary Table S7).The initial observations suggest that few of the SSRs, such as CgSSRs 10, 11 and 14, could be polymorphic (Supplementary Figure S5).

qRT-PCR based gene expression analysis
We performed quantitative real time PCR assay to assess the correlation between the in silico expression data and qRT-PCR based expression analysis.The analysis suggests a good correlation between the two sets of data (Supplementary Table S8) where we used leaf tissue samples as the reference.We then performed qRT-PCR analysis for 8 genes annotated as MYB59 (TRINITY_DN8591_ c0_g1_i3.p1),P2C18 (TRINITY_DN12531_c0_g1_i5.p2), SR43C  In-silico expression of genes reported to play crucial roles in fruit development in cucurbits.The scale represents Log 2 FC between the fruits of Kunkhi and Sabuja varieties generated using RSEM and edgeR.
( T R I N I T Y _ D N 8 6 2 _ c 0 _ g 1 _ i 4 .p 1 ) , F P A ( T R I N I T Y _ DN7200_c0_g1_i4.p1),BRK1 (TRINITY_DN3675_c0_g1_i1.p1), ASIL2 (TRINITY_DN4808_c0_g1_i3.p1), WRKY11 (TRINITY_ DN7218_c0_g1_i4), ENL 2 (TRINITY_DN3745_c0_g1_i1) to observe their expression in the two varieties of Ivy gourd taking ANKu as reference and compared with those of ANSa to make a comparative expression analysis between the developing fruits of the two varieties (Figure 10).The qRT PCR analysis showed significant upregulation of BRK-1 gene in all the developmental stages of ANSa as compared with ANKu, with maximum upregulation of BRK-1 in 5 DAA stage of ANSa.The transcription factor MYB 59 showed slight upregulation in all developmental stages of ANSa.Other genes such as P2C18, ASIL 2 were found to be downregulated in all the stages of development in ANSa.SR43C showed upregulation in only in the 5DAA tissue sample in ANSa while ENL2 was upregulated in 20DAA tissue.Two genes FPA and WRKY11 were generally down regulated or showed no significant change in expression in ANSa.

Discussion
Cucurbitaceae is the second largest family bearing diverse groups of plant with broad range of fruit characteristic (Shaina and Suhara Beevy, 2012;Ma et al., 2020).Studies related to fruit development and variations in fruit sizes are also reported in other plants like Arabidopsis, tomato, pear, pumpkin, cucumber, bottle gourd, apple, Apricot, sweet cherry among many others (Jiang et al., 2015;Xanthopoulou et al., 2017;Di Marzo et al., 2020;Yuste-Lisbona et al., 2020;Zhang et al., 2020;Kuhn et al., 2021;Liu et al., 2023;Wang et al., 2022;Jahed and Hirst, 2023).RNA-Seq analyses have become the most popular method to identify the key genes regulating fruit morphology in several cucurbit vegetables (Zhang et al., 2016;Ma et al., 2020;Chen et al., 2021).Ivy gourd is a wonderful vegetable for its nutritive as well as medicinal value.However, research in this crop is very limited, and genomic information is scarce.Here, we attempt to address the patterns of gene expression that influence fruit size and factors that are heavily influenced by the process of cell division and cell expansion.

Anatomical characterization reveals the role of cell numbers and size in fruit size variation between two varieties of C. grandis
Generally, during development, a fruit undergoes significant growth in size due to cell divisions and cell enlargement that determine fruit size.Interestingly, in the present study, we found Expression profile of the most abundant transcription factors in Ivy gourd transcriptome; (A) Trihelix (B) MYB and (C) Zinc finger C3HC4 type.The heat map was generated using MeV software after averaging the normalized expression values obtained after edgeR.
a similar pattern of growth in cell numbers during fruit developmental stages of both the ANKu and ANSa varieties.The number of cells in fruits of each variety start increasing from 5 DAA to 15 DAA and then stagnate at maturity.This pattern agrees with the normal growth pattern of fruit where more cell division takes place at the early stages of fruit development that determines its fruit size and at maturity cell division slows down or stops to check further growth (Yang et al., 2013).However, the number of cells was found to be higher in ANKu fruits than in ANSa.Reduced cell number are known to result in smaller fruit size in cucumber (Yang et al., 2018) and some previous histological studies also suggest the role of cell number in fruit size determination of melon (Higashi et al., 1999), sweet cherry (Olmstead et al., 2007), peach (Scorzal et al., 1991), pear (Zhang et al., 2007).Therefore, the present study suggests that the longer fruits of ANKu are a result of a higher cell number.Apart from cell number, cell area or expansion also plays an important role in fruit size or shape determination (Liu et al., 2020).In the present study, greater cell area was observed in ANSa compared to ANKu.Maximum expansion was found in 15 DAA and 20 DAA fruit of ANSa and ANKu as compared to earlier stages of development, suggesting that cells expand during mature stages of fruit development before ripening.Such observations lead us to believe that greater cell area and lower cell number in ANSa makes the fruit short and oval or conical while the fruits of ANKu may be long and slender due to larger number of cells and less expansion.Previous studies on apple and pineapple also report fruit size determination through cell number and expansion (Harada et al., 2005;Li et al., 2010).These observations act as a guide to look for the molecular bases of these specific processes for a better understanding of fruit development in Ivy gourd.

Fruit shape and size in C. grandis may be a result of varied gene regulatory networks
So far, we understood that, variation in cell number leads to fruit size variations and that cell number is influenced by cell division.Various genetic factors, hormones can influence the rate of cell division and knowing those factors could help us to identify fruit size regulators.We have incorporated gene expression analysis to identify DEUs that may regulate fruit size in Ivy gourd.Cell division and cell expansion is often associated with reorganization of structural elements such as microtubules and actins (Gibiezǎ and Petrikaite, 2021).It has been reported that BRICK 1 or BRK 1 gene plays an important role in promoting actin polymerization (Djakovic et al., 2006).In our study, BRK1 was up regulated in 5 DAA fruits and 20 DAA fruits of ANSa as compared to the respective stages of ANKu, emphasizing the role of this gene in fruit development.Higher expression of BRK1 in early stage of fruit development in ANSa might have role in extending its width and down-regulation in ANKu probably helps in reducing the width and focusing on increase Expression analysis of selected DEUs (A-H) from fruit tissue of both the varieties ANKu and ANSa by qRT-PCR using 18S rRNA as reference gene.
Values with asterisk (*) are significantly different from other variety at p ≤0.05, as analyzed by unpaired t test.
in length.The higher expression of BRK1 in earlier stages of fruit development also aligns with studies that suggest that the fruit size is determined in early stages of its development (Pan et al., 2020).
Transcription factors like MYB are known to be involved in numerous roles in fruit development, quality, ripening etc (Wang et al., 2023;Zhang et al., 2024).MYB is involved in auxin biosynthesis and regulates fruit size in pear and citrus fruits (Wang et al., 2022;Zhang et al., 2023).MYB 59 belongs to the MYB transcription factor family that promotes cell growth in vegetative parts of plants including roots and shoots (Fasani et al., 2019).In this study, MYB59 showed negligible change in expression between the ANKu and ANSa varieties suggesting that the gene may be involved in processes of fruit development that do not relate to fruit size in Ivy gourd.Protein phosphatases are a group of proteins that have not been extensively studied in plants, however, fruit development and ripening can be regulated by protein phosphatase 2C by some unknown mechanism (Liang et al., 2021).Previous reports also suggest that protein phosphatase negatively controls ethylene biosynthesis (protein phosphatase deficient plants synthesize more ethylene) by regulating the enzyme 1-aminocyclopropane 1carboxylate synthase (ACS) (Park et al., 2021) and ethylene is a negative regulator of cell expansion in most of the tissues (Skottke et al., 2011).In the present study, we have examined the expression of P2C18 gene, which encodes probable protein phosphatase 2C, to be lower in ANSa which suggests that, P2C18 downregulation may promote high ethylene production in fruit tissue of ANSa that lowers the vertical expansion of fruit axis, making the fruit shorter in length.
In other words, upregulation of P2C18 in ANKu may cause low ethylene production thereby supporting higher cell expansion of fruit axis that makes the fruit elongated.ASIL2 is a transcription factors belonging to the trihelix DNA binding protein family.Role of ASIL in fruit development is not yet reported, however it has been reported to have role in seed maturation and vegetative growth (Hussain et al., 2020;Ruiz et al., 2021).In the present study, ASIL 2 is downregulated in all the developmental stages of ANSa, and therefore, can be considered to have higher expression in ANKu, which suggests its role in elongation of fruit by some unknown mechanism.In addition to this, we also observed expression of gene SR43C which codes for chloroplast signal recognition particle and has role in chlorophyll synthesis and fruit quality (Huo et al., 2023) et al., 2010).In the present study, FPA has similar expression in ANSa as compared with ANKu which suggests that the fruit morphology is unaffected by the FPA gene.The WRKY gene family is one of the important families of transcription factors, known to have a role in plant stress and various plant physiological processes, however, no reports have been found about the involvement of WRKY in fruit size regulation (Ali et al., 2018).However, the members of the WRKY TF family are reported to be involved in fruit ripening and enlargement (Liu et al., 2023).In the present study, the WRKY 11 was found to be generally down-regulated in ANSa which suggests that the homolog of WRKY11 in Ivy gourd could have a role in fruit elongation.The expressions of ENL2 gene that encodes early nodulin like protein 2, was also examined in fruit tissues of both the verities.The early nodulin like proteins play major role in transport of nutrients, solutes, amino acids and hormones in plants (Denance et al., 2014).A previous study of watermelon, reported that, downregulation of nodulin like protein 2 inhibits fruit expansion (Gao et al., 2023).In the present study, ENL2, was found to be downregulated in all the stages of fruit development of ANSa except the matured stage (20 DAA).This result can conclude that, fruit expansion is inhibited or slow during early stages of fruit development that makes the fruit short and sudden up-regulation in matured stage may contribute to its lateral expansion making the fruit conical or oval.
There are reports in cucurbits that indicate the process of ethylene biosynthesis is inversely co-related to fruit elongation (Martıńez et al., 2013).In this study, we found the homologs of some of the members of ethylene biosynthesis pathway to be downregulated in the elongated ANKu fruits as compared with the shorter ANSa.It was also interesting to note that another well characterized gene studied in watermelon, Cla011257 (Pan et al., 2020), was downregulated in all stages of ANKu except in 15 DAA.This could be a discrepancy related to the completeness of the transcriptome and warrants further investigation.The next step would be to generate a more comprehensive genomic map of C.grandis to include noncoding elements and complete gene sequences.Such studies have been undertaken in C.grandis where the authors have studied the evolution of Y chromosome (Janousek et al., 2022).Similar efforts need to be made to add to the genomic resources of C.grandis to facilitate future efforts in plant improvement.However, all the above hypotheses are based on the data provided through RNA-seq and molecular characterization using quantitative PCR analyses.These genes need to be characterized for their functional roles in plant system using tools for gene manipulation which have been developed recently for cucurbits (Devani et al., 2020) and can be taken up in future to report a more definitive role of these genes in fruit development in an economically important crop like C. grandis.However, the transformation of cucurbits with target genes is still in its infancy and will require considerable standardizations for Ivy gourd.Therefore, a more feasible approach would be to identify molecular markers associated with the traits of interest and use them for selective breeding.

Transcription factors may be implicated in fruit development and morphology regulation
Transcription factors (TFs) are a varied group of DNA binding protein that interact with cis-acting elements in promoter sequence of target genes thereby regulating gene expression (Cheng et al., 2019).In this study, we have observed that the Trihelix transcription factor (TTF) gene family, also referred to as the GT transcription factor family, is differentially expressed in ANKu and ANSa.The members of this family have been reported to play important roles in overall plant developmental process (Cheng et al., 2019).The ASIL1 gene, a member of TTF DNA binding protein family, regulates silique size in Arabidopsis (Gao et al., 2009).Another TF family that was conspicuous in this study is the MYB transcription factor family, which primarily function in cell cycle regulation, particularly in the G2/M phase (Kranz et al., 2000;Stracke et al., 2001).C3HC4-type RING finger proteins represent one of the largest transcription factors known for various plant processes like regulation of growth and development, signalling network and abiotic stress and fruit development (Wu et al., 2014;Agarwal and Khurana, 2020).Although genes encoding TTF, MYB and C3HC4-type RING finger containing proteins were most abundant in the Ivy gourd transcriptome, given that the unigenes encoding Trihelix TFs showed better differential expression compared with MYB and C3HC4-type RING finger proteins, they could be good candidates for further characterization and it would be interesting to identify molecular markers associated with these genes.
4.4 SSR markers can be developed as a valuable resource for marker assisted selection DNA markers play crucial role in plant breeding programmes to develop new superior varieties with desirable traits and simple sequence repeats (SSRs) or microsatellites markers are one of the most widely studied groups given that maximum mutation rates are observed in SSRs making them a crucial component of genome evolution (Kashi et al., 1997).SSRs are desirable molecular markers with several applications such as determination of functional genetic variations such as paternity determination, genetic diversity assessment, and population genetics studies and for the development of a genetic map (Xanthopoulou et al., 2017).Nowadays transcriptome sequencing is highly evolved that provides a rapid, cost effective and reliable source of expression datasets in non-model species and also facilitates the development of SSRs using bioinformatics tools (Marioni et al., 2008;Zhang et al., 2012).Transcriptomes are a good source of genic SSRs which are linked with the loci that mostly focus on morphology (Lal et al., 2017).In comparison with genomic SSR markers, the genic SSR markers may facilitate the identification of candidate functional genes and could expand the efficiency of marker assisted selection (Gupta and Rustgi, 2004;Zhang et al., 2012).In Ivy gourd, presence of SSRs in the differentially expressed unigenes coding for transcription factors, biosynthesis of metabolites, ubiquitin mediated proteolysis, structural element and signalling pathways which are crucial in fleshy fruit development and ripening (Jia et al., 2023;Wang et al., 2023), indicates their potential for developing functional molecular markers.Our initial assessment also indicates that few of the SSRs are polymorphic in nature, and associated with genes encoding Leucinerich repeat extensin-like protein, Trihelix transcription factor and 3hydroxybutyryl-CoA dehydrogenase, implying that analysis and subsequent sequencing efforts at larger scales will help identify numerous polymorphic SSRs in C.grandis that can be used for marker assisted selection for crop breeding as well as for diversity analysis in plants (Kapoor et al., 2023;Li et al., 2023).

Conclusion
This study reports the transcriptome of Ivy gourd (Coccinia grandis), an important and popular food crop of India, and attempts to explore the differences in patterns of gene expression in two cultivated varieties of the plant with considerably different fruit morphologies.It showcases one of the most widely utilized functions of the modern-day "omics" platforms that allow understanding the molecular basis of important biological processes of an organism.This study involves various stages of fruit development and has led to the identification of genes that may have a crucial role in fruit development as well as in determining fruit size and shape.This study has also generated a variety of genomic resources for Ivy gourd in form of genic SSRs which will provide new avenues for research into fruit morphology in Cucurbitaceae family.These resources will not only enrich the existing repertoire of molecules available for cucurbits, but also shed light on a crop like Ivy gourd, which has great economic potential.Polymorphic molecular markers will prove to be an especially valuable genomic resource for the crop and facilitate future breeding attempts for improving crop quality.Overall, this study provides a foundation for further exploration of molecular basis of fruit size regulation in C. grandis varieties.
FIGURE 4 Functional annotation; (A) Venn diagram representing the annotation of protein coding sequences of C. grandis with various publicly available databases (Uniref, Pfam, COG, KEGG, Swissprot), GO terms classified into three major classes namely Biological Process, Molecular function and Cellular Components (B-D).

FIGURE 6 A
FIGURE 6A comparative analysis of differentially expressed unigenes in both the varieties at each stage of fruit development.S refers to ANSa and K refers to ANKu in the Venn diagram which was created using the webtool at https://bioinformatics.psb.ugent.be/webtools/Venn/.

FIGURE 8
FIGURE 8 Distribution of different transcription factor families in the transcriptome of Ivy gourd.A bar graph representing Ivy gourd unigenes encoding transcription factors belonging to the most well represented families.

FIGURE 7
FIGURE 7 Two varieties of C. grandis with contrasting fruit shape and size were collected from Central Horticultural Experiment Station Bhubaneswar (20.2446°N, 85.7812°E).Seedlings of C. grandis var."Arka Neelachal Kunkhi" (Long, slender fruit) and C. grandis var."Arka Neelachal Sabuja" (Short, oval fruit) were grown in the fields at ILS, Bhubaneswar and flowers were tagged on the day they opened completely.Fruits were collected at four different stages of development: 5 Days after anthesis (DAA), 10 DAA, 15 Abbreviations: DEGs, Differentially expressed genes; qRT-PCR, Quantitative real-time PCR; COG, Clusters of Orthologous Groups; KEGG, Kyoto Encyclopaedia of Genes and Genomes; GO, Gene Ontology; Blast, Basic Local Alignment Search Tool.

TABLE 1
Assembly statistics for the transcriptome of Ivy gourd.

TABLE 2
Details of the SSRs identified in the transcriptome of Ivy gourd.
. Higher expression of SR43C in early stage of fruit development (5 DAA) in ANSa may imply a potential role of SR43C gene in determining fruit size through some unknown mechanisms.Early stage of fruit development is dependent on flowering time, which is regulated by flower time control protein or FPA.It is a nuclear protein that regulates flowering by silencing the mRNA encoding the repressor protein FLC (Hornyik