Background & Summary

The wax gourd (Benincasa hispida), also known as ash gourd, white gourd, Chinese watermelon or winter gourd, is the only species in the genus Benincasa. It is an annual cucurbits native to Asia, and has been used as a vegetable and herbal medicine for thousands of years in China and India1,2. In the recent decades, the wax gourd is grown in more and more areas throughout the world for multi-purpose uses. It contains many important nutrients, and some metabolites can be used in treating fever and various disorders3,4. Commonly, it is used as an important vegetable and its young leaves, flower buds, immature and mature fruit is cooked and eaten. For its medicinal properties, it has been recognized in the traditional Chinese medicine and Ayurvedic medicine system over thousands of years, and now there are an increasing number of studies reported its medicinal values5,6. Moreover, the wax gourd is widely used in the food industry for making candied fruit, moon cakes and many kinds of pies as a base filling material. It is of great importance to broaden our knowledge of the wax gourd for promoting to fully exploit its benefits.

Development of a high-quality reference genome would be very useful for molecular genetics, molecular breeding and evolutionary studies of the wax gourd. Previously, we have reported a draft genome sequence of the wax gourd variety B227, and we revealed that the 12 chromosomes of wax gourd represent the most ancestral karyotype of the investigated cucurbits1. The B227 assembly is the only published de novo assembly of the wax gourd genome so far. Its contigs and scaffolds were constructed based on Illumina paired-end reads (approx. 28-fold), mate-pair reads (approx. 12-fold) and PacBio long-reads (approx. 15-fold), and the final pseudo-chromosomes were generated by anchoring scaffolds onto a published genetic map7. Though developed by combing data from multi platforms, it remains incomplete and highly-fragmented (contig N50 of 68.5 Kb, scaffold N50 of 3.4 Mb), and may contain mis-ordered or mis-oriented scaffolds due to technical limitations of the genetic map based pseudo-chromosome construction as reported in some plant species8,9,10. There is still much left to be improved about the wax gourd reference genome, and the availability of genomes of different varieties would provide more resources that can help to understand the genetic variations and evolutionary history of the crop.

Through our continuous efforts, a high-quality and near-complete reference genome assembly has been achieved. Herein, we report a chromosome-scale high quality genome assembly of the wax gourd variety pf3 by combined use of high-coverage PacBio long reads (approx. 86-fold), Illumina short reads (approx. 50-fold) and the Hi-C data. The de novo genome assembly and annotation workflow is as showed in the Fig. 1.

Fig. 1
figure 1

Overview of genome sequencing, assembly and annotation. Data information is shown in rectangles, software and tools are indicated in italic.

Methods

Sample selection, library preparation and sequencing

Previously we divided wax gourd germplasm into four groups according to their genomic variation data1, which were the wild group (W), the landrace (L), the two cultivated groups (C1 and C2). The wax gourd variety pf3 used for sequencing in this study is an inbreed line developed by us at the Vegetable Research Institute, Guangdong Academy of Agricultural Sciences in Guangzhou, China. It derived from a cross between a small fruit landrace (belongs to L group) collected from Yunnan and a giant fruit elite line (belongs to C2 group). It shows moderate fruit size with wax, high-yielding potential and good quality of taste (detailed morphology see Fig. 2). Fresh seedlings of the pf3 were used for high-quality DNA extraction followed by construction of PacBio SMRT Bell library, Illumina short-read library and Hi-C library. The Bell library was sequenced on the PacBio Sequel II platform (CLR mode), and then the output raw subreads bam file was converted to fastq format, generating 86.53 Gb data (Table 1). The Illumina short-read and Hi-C library was sequenced on the Illumina NovaSeq-6000 platform(PE150), generating 50.92 Gb and 99.55 Gb clean data respectively. All the DNA extraction, library construction and sequencing procedures were performed by the Novogene Company (Tianjin, China) according to the manufacturer’s protocols.

Fig. 2
figure 2

Morphology of the sequenced wax gourd cv. pf3. (a) The whole plant in the field. (b) Female flower. (c) Male flower. (d) Mature fruit. (e) Transection of mature fruit. (f) Seeds.

Table 1 Statistics of the sequencing data of the wax gourd variaty pf3.

RNA sequencing

Root, stem, leaf and flower tissue of the pf3 plants were collected for RNA extraction. Total RNA was extracted from each tissue respectively, and then equal amount of them were pooled together. Thereafter, direct-cDNA sequencing and TruSeq RNA-seq library were constructed using the pooled RNA, and the transcriptomes were sequenced on the Nanopore PromethION and Illumina Hiseq4000 platform by Novogene Company (Tianjin, China), respectively. In total, 50.05 Gb full-length RNA-seq data and 8.33 Gb short-read RNA-seq data were obtained (Table 1). These RNA-seq data were used for whole-genome protein-coding gene prediction.

De novo genome assembly

We first converted raw subreads bam file generated by PacBio sequencer into fastq format using the software BAM2fastx. Statistical analysis showed that the average length of the long reads was 17.68 kb, and the N50 length was 21.01 kb. Then we constructed a primary assembly by using MaSuRCA assembler v4.0.911 with default parameters. The assembler built unitigs (high-confidence contigs) through hybrid assembly strategy with high-coverage PacBio long-read and Illumina short-read data, and generated 1897 unitigs with a total size of 974.87 Mb and unitig N50 of 2.43 Mb (Table 2). We then joined unitigs into scaffolds using Hi-C data via Juicer v1.6 and 3D-DNA v180922 pipeline12 with default parameters. We further visualized the raw scaffolds and conducted manual curation using the Juicebox tool package v1.22.0113. After curation, we obtained 1862 high-accuracy scaffolds with a total length of 975.62 Mb and N50 scaffold size of 70.97 Mb. We designated the assembly pf3 v1.1, and 94.92% of the total length (926.05 Mb) is contained in the 12 largest scaffolds that corresponding to the expected 12 chromosomes of the wax gourd (Fig. 3). By comparison with the B227 assembly we reported previously, great improvement on continuity (contig N50 size of 2.43 Mb vs 68.5 Kb), completeness (975.62 Mb vs 912.95 Mb) and chromosome-anchored size (926.05 Mb vs 859.0 Mb) was achieved in the pf3 v1.1 assembly (Table 2).

Table 2 Summary of comparisons of pf3 v1.1 assembly and B227 assembly.
Fig. 3
figure 3

Genome features of the wax gourd cv. pf3. (a) GC content (30–50%) across 12 chromosomes. (b) Repeat percentage (60–100%). (c) Gene density (0–1452). (d) SNP density (0–24764). (e) Syntenic blocks of paralogous genes. a-d are drawn in non-overlapping 1 Mb sliding window.

Repeat annotation

We masked and annotated repetitive sequences and transposable elements (TEs) in the pf3 v1.1 assembly through incorporating de novo and homology-based predictions. We built a de novo repeat sequences library by using RepeatModeler v2.0.114, and performed homology-based predictions by using RepeatMasker v4.1.2-p115 with the Arabidopsis repeat sequences database. The output cat.gz file (contains list and alignment of repeat regions found in the assembly) of de novo and homology-based prediction was merged together and subjected to post-process with RepeatMasker package to produce the final repeat annotation. In total, we identified 770.68 Mb repetitive sequences in the pf3 v1.1 assembly, accounting for 78.99% of its total length (Fig. 3b and Table 3).

Table 3 Repetitive element annotation statistics.

Gene prediction and functional annotation of the genome

We conducted protein-coding gene prediction with integration of evidence data from ab initio training, transcript and homologous proteins alignment through BRAKER pipeline v2.1.616. We prepared transcript evidence by mapping the full-length RNA-seq and short-read RNA-seq data to pf3 v1.1 assembly by using minimap2 v2.23-r111117 and STAR v2.7.9a18 respectively, followed by sorting alignment bam files with samtools v1.719. In running the pipeline, the soft-masked pf3 v1.1 assembly was used as input genome (–genome option), the sorted RNA-seq alignment files were used as RNA-seq evidence (–bam option), and peptides of the wax gourd B227 assembly and other three cucurbit species (Cucumis sativus ChineseLong v3, Lagenaria siceraria v1 and Cucurbita moschata v1) were used as homologous protein data (–prot_seq option) and aligned to the genome using GenomeThreader v1.7.020 (–prg gth). Briefly, the pipeline started with generating seed genes by GeneMark-ES v4.69_lic with supported by RNA-Seq evidence. Subsequently it performed ab initio training of AUGUSTUS v3.4.021 with the seed genes, RNA-Seq and protein alignment information. Finally, it conducted gene prediction with AUGUSTUS through integrating of the training output, RNA-Seq and homologous protein alignment information. After filtered out genes encoding protein sequence shorter than 50 aa (amino acids), and genes containing internal stop codon, illegal start or stop codon, a total of 37,092 genes were annotated in the pf3 v1.1 assembly (Fig. 3c and Table 4).

Table 4 Summary of gene annotation.

Furthermore, we performed functional annotation of the predicted genes through searching proteins against InterPro database (v88.0). We carried out this by submitted protein sequences to InterProScan 5 webservice22 via a perl script (iprscan5_lwp-nodie.pl) with default parameters. We further used PANNZER223 and eggNOG-mapper v2.1.724 to annotate proteins by Gene Ontology (GO) terms. In total, 31,562 genes were functionally annotated, and 22,707 genes were assigned to specific GO term (Table 4).

Comparison of the assembly of B227 and pf3

To infer the synteny and colinearity the pf3 v1.1 assembly and B227 assembly, we ran the Quick Genome Dot Plot plugin (parameters: Blast e-value 1e-3, Num of BlastHits 5) embedded in TBtools v1.09872625. We found that the two assemblies were quite syntenic (Fig. 4a), but there were large-scale inversions at the ends of some chromosomes. We further analysed the genomic rearrangements and local sequence differences between the two assemblies by using the SyRI v1.626 software. We identified a lot of Mb-sized structural rearrangements including inversions, translocations and duplications, and these rearrangements were mainly located at the end of chromosomes (Fig. 4b).

Fig. 4
figure 4

Whole-genome comparison of the pf3 v1.1 with B227 assembly. (a) Dot plot for the syntenic blocks. (b) Chromosome-level local sequence differences.

Considering the close genetic background between the two cultivars, and that the pseudo-chromosome of B227 assembly was developed based on a genetic map, it could be deduced that most of the rearrangements may be not real existing but rather mis-oriented or mis-placed in the B227 assembly. To check this, we further split the B227 assembly and re-constructed scaffolds (pseudo-chromosomes) by using Hi-C data of a genetically very close cultivar B418 that also belongs to C2 group. We detected much higher collinearity between the Hi-C-based B227 assembly and the pf3 v1.1, and most of large rearrangements discovered previously disappeared (Supplementary Fig. S1). We also examined the differences between the Hi-C based and genetic map based assemblies of the B227, and detected several large inversions at end of chromosomes (Supplementary Fig. S2). These evidences suggest that most of rearrangements showed in the Fig. 4d are indeed errors in the genetic map based B227 assembly.

Discovery of genomic variations

To explore genetic variation pattern of the wax gourd germplasm referenced to the newly developed pf3 genome, 31 representative wax gourd accessions (Supplementary Table S1) that we sequenced previously were selected and subjected to mapping and variation discovery procedure. Sequencing data of each accession was mapped onto the pf3 v1.1 assembly by using bwa v0.7.17-r118827 with default parameters, and the alignment files were sorted and indexed with samtools. We called variations of all the 31 accessions together by using bcftools v1.828 with mpileup and call commands. We applied pre-call filtering with -q 30 and -Q 20 to skip poor mapped reads and low-quality bases when run mpileup command, resulting in an initial total of 36,401,973 variant sites. We evaluated summary metrics of the raw variants, and then filtered them by using VCFtools v0.1.1629 and bcftools based on quality score, depth, average mapping quality and other criteria as described in the Supplementary Note 2. We finally got more than 12 million high-quality variations, including 12,366,466 single-nucleotide polymorphisms (SNPs) and 286,201 small insertions and deletions (InDels). The distribution of SNPs across the pseudo-chromosomes was as shown in the Fig. 3d. Furthermore, we investigated numbers of SNPs in subset of samples, we found that the C2 group contained the minimum SNPs among the four groups (Table 5), and there was only 593,107 SNPs between B227 and pf3.

Table 5 Number of SNPs in subset of samples.

Data Records

The sequencing data, genome assembly and annotation data reported in this paper have been deposited in the Genome Warehouse in National Genomics Data Center (NGDC), Beijing Institute of Genomics, Chinese Academy of Sciences/China National Center for Bioinformation30,31, under the BioProject accession number PRJCA010475 that is publicly accessible at https://ngdc.cncb.ac.cn/gwh. All the clean genome sequencing data including PacBio long-read, Illumina short-read and Hi-C data, as well as RNA sequencing data including Nanopore full-length RNA-seq and Illumina short-read RNA-seq data, were deposited in the Genome Sequence Archive (GSA) of NGDC under the accession number CRA007486. The pf3 v1.1 assembly and annotation data have been deposited in the Genome Assembly Sequences and Annotations (GWH) of NGDC under accession number GWHBJVO00000000. The DNA and RNA sequencing data were also submitted to the National Center for Biotechnology Information (NCBI) SRA database with accession number SRR23081782, SRR23081783, SRR23081784, SRR23081781 and SRR23096591 under BioProject PRJNA89881932,33,34,35,36. The genome assembly has also been deposited at DDBJ/ENA/GenBank under the accession GCA_027475165.137. Sequencing data of the 30 wax gourd cultivars used for genetic variation discovery are available in the GSA under project accession number PRJCA001140, and sequencing information for these cultivars is summarized in Supplementary Table S1.

Technical Validation

Assessment of the genome assembly

To evaluate the completeness of the wax gourd pf3 v1.1 assembly, we first mapped Illumina short-read and PacBio long-reads data back to the assembly, and analysed the alignment file with Qualimap v.2.2.238. The mapping rate of both libraries was above 98% (Table 2), and more than 96.5% of the assembly have at least 20× coverage of Illumina short-read and PacBio long-reads respectively. We then performed the Benchmarking Universal Single-Copy Orthologs (BUSCO) v5.2.239 with eudicots dataset (n = 2,326) to assess the completeness of the assembly. We identified 2,280 complete BUSCOs (98.02%) out of the 2,326 BUSCO groups, including 2,183 complete and single-copy BUSCOs and 97 complete and duplicated BUSCOs (Table 6). The number of fragmented BUSCOs and missing BUSCOs was 9 (0.4%) and 37 (1.5%), respectively.

Table 6 BUSCO assessment results.

Moreover, we evaluated the result of Hi-C based pseudo-chromosomes construction. We mapped the Hi-C data to the 12 pseudo-chromosomes, and then analysed and visualized with Hicexplorer v3.740. As the heat-map of Hi-C contact displays in Fig. 5, the signal intensities of interaction between the two bins were clearly divided into 12 distinct groups, indicating the high-quality of the pseudo-chromosomes assembly.

Fig. 5
figure 5

Hi-C contact map of the chromosome-level assembly of pf3. The intensity of interactions was calculated using a bin size of 10 K.