The study aimed to generate the high-quality reference-level genome of the Indian pomegranate cultivar ‘Bhagawa.’ To achieve this, we opted for short- to long-read sequencing with second- to fourth-generation technologies, such as Illumina, 10X genomics, PacBio, and BioNano Optical mapping. Hence, we could generate enormous genomic data for each technology (Table 1). Using the direct-label and stain (DLS) technology of the Bionano platform, we generated data of 1191.68 Gb covering 340.4 Mb of the genome (∼352.54 Mb of an estimated genome) followed by the use of the HiSeqX platform generating data of 77.08 Gb covering 220 Mb of the genome. The sequence information generated from all these technologies could lead to hybrid scaffolding and the development of a reference-quality genome for ‘Bhagawa.’ The detailed methodology with steps followed to complete the genome is depicted in Figure 1.

From the results of preliminary k-mer genome survey analysis for cv. ‘Bhagawa’ using 10X Chromium datasets, the final genome size was estimated at 352,535,926 bp (352.54 Mb, k = 31). The heterozygosity rate was 0.14% and the repeat ratio was 4.17%, with 65.3% unique sequences and 0.72% error rate, as depicted in Figure 2.

We assembled the genome by using datasets generated from different sequencing platforms and an array of assembling software shown in Table 2. Results of five approaches combining different sequencing technologies and assembly software were compared. Among the five approaches, the Bio-nano genomics (BNG) and BNG-assisted Super Scaffolds coupled with BioNano Access assembly software yielded higher genome coverage (343.27–353.59 Mb) with limited sequences, 78–374, and a higher average N50 of 16.11–16.31 Mb for the scaffolds. Thus, BNG coupled with BioNano Access was the best to obtain high-quality assembly of genomes with higher contiguity than 10X and PacBio contigs coupled with Super Nova, Falcon-unzip, and Scaff10x assembly (Table 2). Furthermore, we polished the genome assembly to obtain the final assembly by using PBJelly2 and Pilon. The final reference-quality genome was 346.08 Mb with only 342 scaffolds and an average N50 of 16.12 Mb and N90 of 1088.62 Kb (Table 2). This assembly covered 98% of the estimated size of the pomegranate genome, 352.54 Mb.

We assessed the genome quality by using the LTR Assembly Index (LAI) based on LTR retrotransposons that account for the largest genome component in most plant genomes. Our genome had an LAI index of 10, indicating a reference-quality genome as per the LAI index scale of ≤ 10 to 20 for reference-quality genomes (Ou et al., 2018). Finally, the qualities of different assemblies were assessed by using BUSCO. The complete genomic landscape of the ‘Bhagawa’ genome is presented in Figure 3.

To judge our genome’s quality, we compared the assembly statistics of the ‘Bhagawa’ genome with the other two draft genomes for pomegranates cv. ‘Taishanhong’ and ‘Dabenzi.’ According to the assembly parameters listed, our genome qualifies for all parameters with high standards as compared with the previous genomes (Table 3). For instance, our final assembly covered 346.08 Mb (98% of estimated genome size) in 342 scaffolds with an N50 of 16.12 Mb and a contig N50 of 6.8 Mb with high GC content (41.01%) (NCBI database Bio project: PRJNA562100, PRJNA505392, PRJNA505397, PRJNA505398, PRJNA505582, and PRJNA445950). A total of 29,435 gene models were discovered, with a mean transcript length of 2,954 bp and a mean coding sequence length 1,090 bp as compared with the previous draft genomes (Supplementary File 1). The LAI index that was calculated for all three genomes was 10 for ‘Bhagawa,’ 8 for ‘Taishanhong,’ and 2.5 for ‘Dabenzi,’ so our genome is of high reference quality as compared with the other draft pomegranate genomes. However, more recently, a chromosome-scale genome assembly of soft-seeded pomegranate cv. ‘Tunisia’ was reported; the assembly involved the combined use of single-molecule sequencing and high-throughput chromosome conformation capture techniques. The genome covers 320.31 Mb, with 39.96 Mb scaffold N50 value and 4.49 Mb contig N50 value, and the genome includes 33,594 protein-coding genes (Luo et al., 2020).

To assess the syntenic relationships between the newly assembled ‘Bhagawa’ genome and the two previously reported draft genomes, we performed global whole-genome alignments using Minimap2 with parameter “asm10”; the synteny relationships were visualized as a Jupiter Plot (Jones et al., 2014). A Circos-based genome assembly consistency plot was used to view large-scale translocations and other large structural variations (Krzywinski et al., 2009; Jones et al., 2014). The connecting bands within the circle represent regions of synteny, whereas the blocks on the arc of the circle represent the largest scaffolds in the assembly. The lack of diagonal lines extending from the middle of the scaffold block suggests no definite breaks in synteny between the two assemblies at 10-Kb resolution. We found perfect macro synteny between the ‘Dabenzi’ vs. ‘Bhagawa’ genome compared to ‘Taishanhong’ vs. ‘Bhagawa’ (Figure 4). These results indicated higher syntenic relationships between ‘Bhagawa’ and ‘Dabenzi’ followed by ‘Taishanhiong.’ Several scaffolds had multiple hits in the ‘Taishanhong’ vs. ‘Bhagawa’ plot, thus indicating low levels of genome rearrangements between these two genomes. Similar results were observed with the ‘Dabenzi’ genome plotted against ‘Taishanhiong.’

An SV represents a major source of genetic diversity, we investigated overall SVs in the pomegranate genome. We aligned the ‘Bhagawa’ genome with the ‘Taishanhong’ and ‘Dabenzi’ genomes separately by using Assemblytics (Supplementary Figure 2). Here, we tried to identify SVs for six different classes: insertion–deletion, repeat expansion and contraction, and tandem expansion and contraction, focusing on 1–10 Kb between our ‘Bhagawa’ genome and the other two draft pomegranate genomes. We identified 5.7 Mb and 9.2 Mb SVs for ‘Taishanhong’ and ‘Dabenzi’, respectively, with ‘Bhagawa’ as a reference (Supplementary Tables 1A,B and Supplementary Figure 2). Although, all the three cultivars belong to the same species, P. granatum, they have accumulated much higher SVs at the nucleotide level in the genome. The SV graph indicates that except for tandem contraction, much higher SVs were observed between ‘Bhagawa’ and ‘Dabenzi’ than between ‘Bhagawa’ and ‘Taishanhong’ for all the variant classes. Also, the SVs frequency decreased with the increasing size of the nucleotide sequences analyzed for each class. Among the SVs, insertions and deletions contributed much greater to total variations in the genomes, followed by repeat expansion and contraction.

Most plant genomes have highly repetitive regions with transposable elements. We have analyzed the pomegranate genome for these aspects to reveal the basic structural features of the genome. De novo analysis of repetitive elements using the Repeat Modeler software revealed a large proportion of repetitive DNA, comparable to the other higher eukaryotic genomes (Table 4). Repeat Modeler generated 1,765 different families of repeats known as the de novo repeat reference library. A total of 496,570 repetitive elements were in the genome, covering 179,092,850 bp of sequence. Most of the repetitive elements were the interspersed type, consisting of Class I (retro transposons) and Class II (DNA transposons) elements followed by unclassified elements, small RNA, and satellites. Simple direct repeats and low complexity repeats represented only 1.03% and 0.18% of the total repetitive elements. Classification of the observed transposable elements into known classes revealed that most repetitive sequences were retrotransposons (18.27%), whereas 1.10% were DNA transposons (Table 4). The most abundant repeats were long-terminal repeat elements (17.33%), of which 13.57% were Gypsy-type elements and 2.95% Copia-type elements (Table 4).

The complete genome sequence of ‘Bhagawa’ (∼346 Mb) was analyzed for gene prediction ab initio and by using a homology-based approach with the BRAKER/MAKER pipeline. A total of 29,435 genes were predicted with an average size of 2,954 bp, with average exon and intron sizes of 286 and 368 bp, respectively (Supplementary File 2). All predicted genes were functionally annotated by following a consensus approach of known homologous or predictive sequence signatures by Cluster of Orthologous Groups, Gene Ontology, InterProScan, Kyoto Encyclopedia of Genes and Genomes, Uniprot, and EggNOG (Figures 5, 6). In total, 96.04% of genes had sufficient similarity entries in databases to tentatively assign gene functions. Only 3.96% of genes remain unannotated. The overall GC content of the genome was 41.01% (Table 3). Additionally, 617 transfer RNA genes were predicted in the genome.

The complete genome represented 93.68% of the 1440 ortholog genes, with 64 missing and 27 fragmented genes, according to step-wise BUSCO assessment of the Embryophyta lineage. We also successfully validated the four transcriptome data samples (NCBI database: SRR5187757, SRR5187758, SRR5187763, and SRR5187764) of bacterial blight-challenged pomegranate tissues by mapping them onto the reference genome. As a result, 85% of 95% of reads showed significant matches to the genome (Supplementary Figure 3) (Kozlov et al., 2019). These results indicated an extremely low proportion of misassemblies in the gene-rich regions. The transcriptome assembly does not necessarily represent all pomegranate genes; it is restricted to the tissue types used, and genes expressed at low levels are likely under-represented.

To reconfirm the phylogenetic position of our newly sequenced pomegranate, cv. ‘Bhagawa,’ and to identify expansion, contraction, and rapidly evolving gene families in the Punica clade in relation to four other sequenced plant genomes (Supplementary Figure 1), we performed whole-genome phylogenetic analysis. For this, we used CAFE software after pruning the gene families containing a single species. From these analyzes, pomegranate and E. grandis diverged 64 (60–70) million years ago (MYA), after the paleotetraploidy event (109 MYA) was identified in the E. grandis genome (Figure 7). Therefore, this whole-genome duplication event is shared by pomegranate and E. grandis (Myburg et al., 2014; Yuan J. et al., 2018). We also found a higher divergence period between the grape and pomegranate genomes, which suggests that pomegranate and grape did not undergo a recent genome duplication as per results from syntenic block analysis between pomegranate, Eucalyptus, and grape genomes (Myburg et al., 2014; Yuan J. et al., 2018).

To study the gene families that had expanded or contracted in pomegranate, we have estimated and compared gene gain and loss rates (λ) using the two-lambda model for gene families from pomegranate with those of four other representative ancestral species. In total, 965 gene families had expanded in number; 5,406 gene families had contracted in number; and 115 families had rapidly evolved in the pomegranate genome, with λ value 0.0059 as compared with its most recent common ancestor, Eucalyptus (Supplementary Table 2 and Figure 7). However, the grape genome had the most divergent species with gene family evolution (λ 0.0038) as compared with the other genomes.

Resistance genes in the ‘Bhagawa’ genome were identified by using the protein sequence of predicted genes and the Plant Stress Protein Database (PSPDB). A total of 1,573 protein-coding genes were identified as resistance genes, classified into 32 domains. Among these, Kinase, TM (577); Kinase, LRR, TM (203); LRR, TM (194); and kinase (132) domains increased in number (Table 5). Transcription factors were predicted by using PlantTFDB, with 1,533 proteins identified as transcription factors in ‘Bhagawa.’ These transcription factors are divided into 57 transcription factor families. Genes related to MYB were highest in number, followed by bHLH, ERF, NAC, and C2H2 families (Table 6).

The masked genome was used to identify simple sequence repeats (SSRs) by using two different tools (MISA and PERF). Since PERF does not identify complex/compound SSRs, we retained outcomes from MISA only. A total length of 2.52 Mb or approximately 0.7% was identified as SSRs in the genome. Then SSRs were divided into three categories according to their length as follows: (1) ≤ 12 bp, (2) 12–19 bp, and (3) ≥ 20 bp. SSRs ≤ 12 bp were higher in number than the other SSR categories, and mononucleotide SSRs were in very high abundance (Table 7).

miRNAs are important regulators of several biological processes, such as plant growth and development. These are 20–24 nt in length. A total of 314 copies of miRNAs belonging to 26 different families were identified in the ‘Bhagawa’ genome (Figure 8).

The most commercially grown and popular Indian pomegranate cultivar ‘Bhagawa’ was chosen for genome sequencing. The variety is currently being maintained at field Gene Bank of ICAR-NRCP, Solapur, India. This variety has a medium plant height (1.5 mt), bearing deep red fruits with large size (>250 g), thick rind having dark red-colored bold arils with the soft seed type. The soft red arils of this variety are mainly preferred for edible and processing purposes because of their higher juice percentage (45%), TSS (15.9°Brix), vitamin C (14.60 mg/100 g fresh arils), anthocyanin (360 mg/100 g), iron (0.32 mg/100 g), and zinc (0.50 mg/100 g) content. This variety is a late-maturing type (180 days from flowering to harvesting) with a high yield potential (20 tonnes/ha) widely grown across the country. The newly emerged leaves were taken for DNA isolation and subsequent downstream operations.

High-molecular-weight genomic DNA (>50 Kb) extraction and purification were performed using the Genomic-tip 100/G genomic DNA isolation kit (Qiagen). DNA concentration was measured by using NanoDrop (Thermo Fisher Scientific) and Qubit 2.0 (Invitrogen) instruments. The integrity and quality of genomic DNA were confirmed by using the Bio-Rad^® CHEF Mapper^® XA Pulsed Field Electrophoresis system. High-quality sequencing data were generated using four different sequencing and mapping technologies (i.e., Illumina, 10X Chromium, PacBio Sequel, and DLS BioNanoSaphyr optical mapping). Illumina sequencing libraries were constructed by using the NEBNext UltraTM DNA Library Prep Kit (Illumina). The 10X Chromium genomic libraries were prepared by using the Chromium Genome HT Library Kit and Gel Bead Kit v2. The SMRTbell library was prepared using SMRTbell Express Template Preparation Kit. For the construction of optical genome maps, standard BioNano protocols, nicking, labeling, repair, and staining processes were implemented. Specifically, DNA was digested by the single-stranded nicking endonuclease Nt.BspQI. Optical maps were assembled with BioNano Irys View analysis software; only single molecules with a minimum length of 100 kb and six labels per molecule were used.

We performed a genome survey using the 10X Chromium datasets before the genome assembly. Initially, the data were de-barcoded by using scaff_reads from Scaff10x v4. Kmer-genie was used to generate histograms at multiple k-mers. The histogram from k = 31 was used to run GenomeScope v1.0 at a maximum k-mer coverage of 10,000 for estimating the genome size, heterozygosity, and repetitive and unique content of the genome (Vurture et al., 2017). After this, we proceeded with genome assembly step by step by using various datasets obtained from the sequencing. The genome was assembled by using the Falcon and Phased falcon-unzip module (Chin et al., 2017). The unzipped genome was polished by using Arrow as a part of the unzipping pipeline. Then the genome was scaffolded iteratively twice by using Scaff10X v4³ with the de-barcoded 10X Chromium reads. BioNano Tools 1v.3.8041.8044 and Solve 3.3_10252018 were used to build a consensus genome map by using molecules > 150 Kb long and hosting 8 labels (Stankova et al., 2016).

The de novo assembled Optical Maps were then used to super-scaffold the Scaff10X scaffolded genome by using the Hybrid Scaffolds tool from the same BioNano solve release (Stankova et al., 2016). Gaps introduced by both Scaff10X and Hybrid Scaffolds toolkits were filled by using PBJelly2 once and later by polishing using Pilon, and BUSCO v3 used with Embryophyta lineage constituting 1440 orthologs was performed at each step to determine the completeness of the genome assembly (English et al., 2012; Walker et al., 2014; Simao et al., 2015). Additionally, on the Pilon-polished genome, we mapped transcriptome data (SRR5187757, SRR5187758, SRR5187763, and SRR5187764) of four samples by using HiSAT2 to check for gene-model completion and to aid in the annotation process (Kim et al., 2014; Geib et al., 2018). Deep-Sequenced Illumina short-read datasets were mapped to the genome by using BWA, and subsequent alignments were used to polish the genome by using Pilon. Assemblies from the ‘Dabenzi’ and ‘Taishanhong’ cultivars were aligned to ‘Bhagawa’ by using the “asm5” parameter from minimap2. Various other genomes, Arabidopsis thaliana, Vitis vinifera, E. grandis, and Theobroma cacoa, were aligned with ‘Bhagawa’ by using the “asm10” parameter. Variants were called from the alignments by using Assemblytics (Nattestad and Schatz, 2016), commit df5361f from GitHub from sizes 1–10,000, with a unique anchor length of 10 kb.

Circos Jupiter graphs were plotted to check the assembly continuity and any misassembles between the ‘Dabenzi’ and ‘Taishanhong’ genomes with respect to the ‘Bhagawa’ genome. We used the largest scaffolds representing 75% of the genome length, consisting of scaffolds > 100 Kb long and with no breakages for at least 50 Kb. Scaffolds were broken into fragments if there was a continuous stretch of 100 Kb, then “N”’s had to be removed (i.e., false alignments). The alignments were computed by using Minimap2 with the parameter “asm10”⁴.

To identify and classify different repeats in the pomegranate genome sequences, we used the RepeatModeler-open-1.0.10 pipeline to construct a de novo repeat library (Smit et al., 2008). Repeat sequences related to the Punica genus were obtained from the RepeatMasker library to classify the repeats by using RepeatClassifier. De novo identified repeat sequences, and plant-related repeat sequences were merged to create a custom library. Masking of the genome involved using RepeatMasker v 4.0.9 with the custom library (Tarailo-Graovac and Nansheng, 2009).

We created an orthologs protein set of 4 plant species, including V. vinifera, P. granatum (‘Dabenzi’), T. cacao, and E. grandis by using OrthoVenn 2 (Wang et al., 2015). Ortholog protein sets shared by at least two species were submitted to BRAKER2.0 using –prg = gth and –prot_seq options. GTFs from the mapped transcriptome data were merged by using StringTie v1.3.6 to create a catalog of putative transcripts (Pertea et al., 2015). BRAKER2.0 was used to train Augustus v3.3 with the merged GTF and protein alignments (Stanke et al., 2006). Gene models were further improved via MAKER2 by providing the gff file generated by using Augustus v3.3, the protein alignment file generated through genome threader (gth), and a CDS file generated by using StringTie v1.3.6; MAKER2 was run iteratively twice to refine the gene predictions (Cantarel et al., 2008).

The predicted gene set was annotated by using EggNOG mapper and blasted against the UniProt⁵ Viridi plantae dataset and InterProScan to filter out false-positives generated by MAKER2 (Van Nguyen and Dominique, 2009; Jones et al., 2014; Huerta-Cepas et al., 2017). The curated gene sets were re-annotated by using the same databases. Finally, genome Annotation Generator (GAG) v2.0. was used to add start and stop codons and identify various metrics based on the annotated GFF (Geib et al., 2018). Genes flagged as overlapping or contained by GAG were removed by using BedTools Intersect.

High-throughput SSRs were identified by using MISA (Beier et al., 2017) and PERF v0.2.5 with the same parameters (Avvaru et al., 2017). The parameters used were minimum SSR motif length of 10 bp and repeat length mono-10, di-6, tri-5, tetra-5, penta-5, and hexa-5; the maximum size of interruption allowed between two different SSRs in a compound sequence was 100 bp. Concordant SSRs between the two approaches were chosen by using Bedtools Intersect (Quinlan, 2010) with a 90% reciprocal overlap.

We used Orthofinder 2.3.3 to define gene families across five genomes: P. granatum cv. ‘Bhagawa,’ A. thaliana, E. grandis, T. cacao, and V. vinifera downloaded from NCBI(Emms David and Steven, 2015). We used the options Diamond to blast, MSA using MAFFT, and tree inference using raxml-ng (Katoh et al., 2005; Buchfink et al., 2015). CAFE 4.2.1 was used to identify expansion, contraction, and rapidly evolving gene families in the Punica clade (De Bie et al., 2006; UniProt Consortium, 2014; Kozlov et al., 2019; Singh et al., 2021). We estimated gene gain and loss rates with a two-lambda model.