Leaf tissues from multiple pomegranate genotypes ( Table S1 ) obtained from the collection garden of the Genetic Resources Institute of the Azerbaijan National Academy of Sciences were used for DNA extraction. For RNA-Seq, multiple tissue sources, including leaves, stems (new and old), petals, flowers, and roots, were harvested at various time points. The harvested tissue was snap-frozen in liquid nitrogen and stored at −80°C until performing DNA/RNA extractions. All plant materials can be requested from the Genetic Resources Institute, Ministry of Science and Education, Azerbaijan.

DNA extraction was performed using the DNeasy^® Plant Mini Kit (Qiagen, Hilden, Germany) as per the manufacturer’s instructions. Quality control was performed by resolving the DNA sample on a 0.8% (w/v) agarose gel. For RNA extraction, the method described by Zarei et al. (2012) was followed except for resuspending the RNA pellet in 20 µl of RNase-free TE buffer. The extracted RNA samples were resolved on a 1.2% (w/v) denaturing agarose gel to assess the integrity of the nucleic acid. Both the isolated DNA and the total RNA samples were quantified using a spectrophotometer (Thermo-Scientific, Wilmington, DE, USA) at wavelength ratios of A260/280 and A260/230.

Cytogenetic studies were performed on four P. granatum local and two introduced cultivars (Gizili (Azerbaijan), Purpursid (USA), Goynar (Azerbaijan), Valas (Azerbaijan), Achygdona (Uzbekistan), and Fatima (Azerbaijan)) in 2021. Forty to 80 young flower buds (long- (bisexual), middle- (intermediate), and short (male) pistiled) containing meiotic divisions in the PMCs (pollen mother cells) were randomly collected between 9 and 12 a.m. and fixed in glacial acetic acid: ethanol 96% (1:3) for 24–72 h. After thorough washing, the experimental materials were transferred to 70% ethanol and placed in a cool place (+5°C) to be used for subsequent studies. Anthers were stained with 2% aceto-carmine, dissected out, and squashed on microscopic slides in 45% acetic acid medium (Darlington and Lacour, 1969). The cover slip was placed on it and heated on a hot plate repeatedly to avoid boiling the cell suspension. After heating, the slides were pressed firmly to get a good spread of chromosomes. The slides were then sealed and observed under a microscope (Nikon Eclipse) with oil immersion lens.

For the genotype Azerbaijan guloyshasi, two paired-end libraries with 350 and 550 bp inserts were prepared using the Illumina TruSeq^® DNA PCR-Free Library Prep kit as per the manufacturer’s instructions. To further generate a higher diversity of fragments with 2–15 kb inserts, one mate pair library was prepared using Illumina’s Nextera^® Mate Pair Library Prep kit following the Gel-Free protocol as detailed by the manufacturer. The libraries were assessed on the Agilent TapeStation 2200 platform with High Sensitivity D1000 ScreenTape (Agilent Technologies, Santa Clara, CA, USA) following the manufacturer’s protocol. The assessed libraries were further quantified using the KAPA library quantification kit (KAPA Biosystems, Boston, USA) according to the protocol described by the manufacturer prior to subjecting them to sequencing on three lanes of the HiSeq 1500 system (Illumina Inc., San Diego, USA). To minimize sequencing errors for assembly, raw reads from the Azerbaijan guloyshasi genotype were filtered for adaptor contamination and poor quality. Any reads containing more than three consecutive Ns or more than three nucleotides with a PHRED score of ≤20 or a median PHRED score of <20 or a read length of <50 nucleotides for paired end and 25 nucleotides for mate pair data were trimmed. De novo assembly was performed using SOAPdenovo software 1.05_127mer (Li et al., 2008). Reads were first assembled from the short insert size (c. 300–500 bp) libraries into contigs using de Bruijn graph k-mer overlap information. The long insert sizes of mate-paired libraries (≥2 kb) were not initially used, as the chimeric reads common to such libraries can generate incorrect sequence overlaps. To avoid this problem, a hierarchical assembly method was used through step-by-step scaffold construction by data addition from each library separately ranked from smallest to largest according to insert size. The number of reads can be found in Table S1 .

A total of six RNA-Seq libraries were prepared with the Illumina TruSeq^® Stranded RNA LT Kit according to the protocol described by the manufacturer. The libraries were assessed on the Agilent TapeStation 2200 platform with D1000 ScreenTape (Agilent Technologies, Santa Clara, CA, USA) following the manufacturer’s protocol. Each RNA-Seq library was prepared with a unique indexing primer, and all six libraries were multiplexed at an equimolar concentration obtained from the TapeStation to generate a single pool. The multiplexed pooled sample was quantified using the KAPA library quantification kit (KAPA Biosystems, Boston, USA) according to the protocol detailed by the manufacturer. The quantified sample was subjected to pair-end sequencing using the HiSeq 3000 system (Illumina Inc., San Diego, USA).

The raw sequence reads were filtered using a custom perl script and Cutadapt v1.9 (Marcel, 2011). Low-quality reads (reads with >10% bases with Q ≤20) and adaptor sequences were removed from the sequenced reads. Data trimming involved the removal of reads that had three or more consecutive unassigned Ns with a PHRED score of ≤20. Prior to the de novo transcriptome assembly step, sequence reads that were less than 50 bp were also discarded. The filtered data was primarily assembled using the transcriptome assembler, SOAPdenovo-TRANS (Xie et al., 2014), with a k-mer size of 101. Further, bubble, fork, and complex loci obtained from SOAPdenovo-TRANS assembly were further combined using the CAP3 assembler (Huang and Madan, 1999) with a minimum overlap of 50 bp and 95% identity to generate longer sequences. Transcripts having a total length of less than 240 bp were removed, as these were considered shorter than the length of a single pair of the sequence. The assembled transcripts were aligned with BLASTX against the NCBI non-redundant (Nr) protein database and the UniRef100 database version 1.0 with an E-value of <10⁻¹⁰.

Pairwise whole genome alignment between different de novo assembled reference genomes was conducted using MUMmer V3.23 (Delcher et al., 2002) with a maximum gap length of 500 bp and a minimum cluster length of 100 bp. The resulted coordinates were visualized using the R package “circlize” (Gu et al., 2014) to present the synteny between chromosomes/scaffolds for different reference genomes.

Whole genome sequence alignment for the five de novo assembled cultivars Azerbaijan guloyshasi, Bhagawa, Taishanhong, Tunisia, and Dabenzi was conducted using the software MUGSY V1.2.3 (Angiuoli and Salzberg, 2011) with its default parameters. Scaffolds with a length of <10 kb were excluded from this analysis. Only aligned segments that existed in all five genomes were used to avoid bias due to incompleteness of genomes (e.g., the assembly size for Taishanhong was only 274 Mb while the estimated genome size was 322.7 Mb). The phylogenetic tree was constructed using IQ-TREE software V1.6.5 (Nguyen et al., 2015) with the model finder option “mtree” which finds the best model from 286 DNA models. The best model was “GTR + F + I” which represents the General time reversible model (GTR; Tavaré, 1986) with unequal rates (+I: estimate the proportion of invariable sites to consider the heterogeneity rate between lineages and across sequence positions) and unequal base frequencies (+F: Empirical state frequency observed from the data). Phylogenetic trees were visualized with FigTree v1.4.4 software (https://github.com/rambaut/figtree/).

Six cultivars were selected for pan-genome analysis and re-sequenced with coverage ranging between 10× and 80×. The dataset included four cultivars from Azerbaijan, one from Uzbekistan, and one from the USA, plus the de novo assembled genomes (Azerbaijan guloyshasi, Taishanhong, and Dabenzi). The reads for each cultivar were mapped to the Tunisia genome using the software BOWTIE2 V2.3 (Langmead and Salzberg, 2012). Unmapped reads were pooled and de novo assembled for each cultivar independently using MaSuRCA V4.0.2 software (Zimin et al., 2013) to produce new sequences that are not located in the Tunisia genome, and only contigs with size >1,000 bp were kept. Assembled contigs were compared to all available bacterial, plastid, and mitochondrial genomes on NCBI (downloaded on 21 August 2021). Only hits with E-values <10⁻⁸ and length >100 bp were kept. Contigs that had aggregated matches >25% of their lengths with >70% identity were filtered out and considered contamination. Remaining contigs were annotated using MAKER2 software (Holt and Yandell, 2011), while Augustus (Stanke et al., 2006) and SNAP (Schmid et al., 2003) software were used for de novo gene prediction using publicly available RNA and protein data for pomegranate (downloaded on 23, 8, 2021) as well as the RNAseq data produced in the present study. The cleaned and non-reference contigs were pooled, and redundant sequences were eliminated using the software CD-HIT (Li and Godzik, 2006). The remaining contigs were re-aligned to the Tunisia reference to ensure that they do not exist in the reference. These contigs and their annotations were pooled and added to the reference Tunisia to define the gene presence and absence variation among the cultivars and to build the graph reference. The graph reference was generated using vg V1.37 software (Garrison et al., 2018).

Whole genome data for the eleven pomegranate cultivars (six re-sequenced cultivars as well as five de novo assembled genomes) were mapped to the Tunisia reference genome after adding the pooled non-reference contigs produced in the previous step using BBMap V35.85 software (Bushnell, 2014). The resulted BAM files were analyzed with the software SGSGeneLoss V0.1 (Golicz et al., 2015) to determine the gene presence and absence variation of the pooled genes.

Genotype likelihoods at each genomic position from the BAM files generated when aligning the reads of each cultivar to the Tunisia reference were calculated using the “mpileup” command implemented in bcftools (Li, 2011). The output was processed with the bcftools command “call” to generate a VCF file for each genome with the -m option (alternative model for multiallelic and rare-variant calling). Finally, the VCF files for all cultivars were merged with the “merge” command implemented in bcftools.

The re-sequenced cultivars in the present study were added to the 26 cultivars reported in Luo et al. (2020). Linkage disequilibrium decays for soft-seeded and hard-seeded subpopulations were estimated as the r2 value (Hill and Robertson, 1968) calculated between each pair of loci located on the same chromosome with the R package “snpstats” (Clayton, 2015). The r2 values were plotted against the physical distance between each pair of loci, and the second-degree Loess smoothing line was fitted and plotted using R. Principal component analysis (PCA) was calculated with PLINK V1.9 software (Purcell et al., 2007), and the first two PCs were plotted with R. Population structure was inferred using the software ADMIXTURE (Alexander et al., 2009). A hundred replicates of the analysis were conducted, with the number of underlying subpopulations (K) ranging from 2 to 10 with 20 cross-validations. The optimal K value was defined as the lowest average cross-validation value over the 100 replicates, which was K = 3 ( Figure S2 ). However, different replicates showed the lowest cross-validation value at K = 4, so we plotted the results for K = 2, 3, and 4. The phylogenetic tree was constructed using IQ-TREE software V1.6.5 (Nguyen et al., 2015) with the model GTR + ASC, which is designed to handle ascertainment bias in SNP data. The analysis was run with 1,000 bootstraps and 1,000 replicates of the SH-like approximate likelihood ratio test (Guindon et al., 2010).

Genomic loci that were differentiated between hard-seeded and soft-seeded subpopulations were selected to have the reference allele fixed (or have a maximum of one heterozygote cultivar) in one subpopulation while the other subpopulation had the alternative allele fixed (or have a maximum of one heterozygote cultivar). The advantage of repeating this analysis in our study over that reported in Luo et al. (2020) is that we included the soft-seeded cultivar “Purpursid,” which was genetically clustered with the hard-seeded cultivars. Therefore, it is expected to have fewer and shorter differentiated genomic regions between soft- and hard-seeded subpopulations compared to the results presented in Luo et al. (2020).

Cytological studies on the studied pomegranate cultivars showed that they all possessed 2n = 16 chromosomes, and at diakinesis and M1, the association of 8II was noted ( Figure S1 ). Table S1 shows a summary of sequencing reads for each cultivar. Based on empirical testing for genome assembly, an optimal k-mer size of 91 was used for SOAPdenovo assembly. Initially, reads from short insert sizes were assembled into 908,448 contigs, with an N50 of 1,568 bp. After adding the long insert size mate-paired libraries, the final assembly consisted of 12,363 transcripts (>500 bp), representing a cumulative length of 371.6 Mbp with an N50 of 220,197 bp ( Table S2 ). For RNA-Seq data, the optimal k-mer size of 101 was identified after empirical testing. Initial assembly resulted in 78,359 transcripts, representing a cumulative length of 67.7 Mbp with an N50 of 1,657 bp ( Table S2 ). Following that, the CAP3 program was used to analyze and assemble a set of scaffolds that were identified as specific loci and contained multiple sequence entries described as forks, bubbles, or complexes. The final assembly consisted of 58,064 transcripts after filtering transcripts <241 bp, with a total assembly length of 56,497,814 bp and an N50 length of 1,561 bp ( Table S2 ). All transcripts from this assembly were BLASTX against the non-redundant (Nr) and UniRef100 databases, allowing the identification of a total of 29,216 transcripts (50%) with significant matches to protein databases (28,972 and 29,019, respectively).

The genome of the AG2017 cultivar showed a very high synteny with the cultivar Tunisia, which was assembled to a chromosome level. All scaffolds were completely and continuously mapped to a single chromosome except for scaffold4026, which was split into the centromeric region of chromosome 4 and the telomeric region of chromosome 7. Figure 1A shows the relation between scaffolds >250 kbp (400 scaffolds) of the cultivar AG2017 and Tunisia chromosomes. These large scaffolds seem to be in the telomeric regions, while short scaffolds are concentrated more in the centromeric regions. These two cultivars have eight diploid chromosomes. Interestingly, they both showed very high synteny with the Chinese cultivar Dabanzi, which has nine pairs of chromosomes, as well as the Indian cultivar Bhagawa, which has eight pairs of chromosomes. However, when compared to the other Chinese reference with nine diploid chromosomes, “Taishanhong,” 11 scaffolds were mapped to two different chromosomes, while one scaffold, “MTJX01000087.1,” was mapped to three chromosomes (1, 4, and 8). Figure 1B presents the synteny between these scaffolds with the Dabanzi and Tunisia sequences. It is clear from the figure that the breakpoints of these scaffolds have been mapped to different locations on both genomes, and the synteny of these locations is matched on the Dabanzi and Tunisia genomes.

A whole genome alignment of the five cultivars (only scaffolds >10 Kbp) was conducted, and only genomic regions that were found in the five genomes were used to infer the phylogenetic tree. The total length of the matched sequences was 267.7 Mbp, which is equal to 98.9% of the total size of Taishanhong scaffolds with length >10 Kbp (270.8 Mbp) that were used in the analysis. Dabenzi and AG2017 have the highest relatedness among all cultivars ( Figures 1C, D ), followed by Tunisia and Bhagawa. This is a very interesting finding, giving that both cultivars have different numbers of chromosomes. Tunisia had an equal distance from both cultivars, while Taishanhong showed the highest distance from all cultivars. However, Taishanhong was most related to the Dabenzi cultivar, and they both have nine diploid chromosomes ( Figure 1D ).

The pan-genome dataset (six re-sequenced and four de novo assembled cultivars) added up to a total of 1,941 contigs to the Tunisia reference genome with an aggregated sequence length of 4.1 Mbp ( Table 1 ). The number of contigs ranged from 25 for Achygdona to 1,001 for Taishanhong, with total lengths ranging between 34.6 Kbp and 1.8 Mbp, respectively ( Table 1 ). These results indicated very limited insertion/deletion variation among the studied cultivars despite the wide geographical and botanical variations in the population. The variant summary between the re-sequenced panel and the Tunisia genome can be found in Table S3 . The total number of genes in the pan reference was 31,016 of which 25,210 (81.3%) existed in all eleven cultivars. A total of 552 (1.8%) were missing from only one cultivar, of which 324 (58.7%) were missing from Taishanhong and Tunisia ( Figure 2C ). The previous two categories were considered core genes (83.1%) that existed in almost the whole population ( Figure 2A ). On the other hand, 3,916 (12.6%) genes were cultivar-specific genes that existed in a single cultivar only, of which the majority (3,805, 97.2%) belong to the cultivars Taishanhong and Tunisia ( Figure 2C ). These genes were assumed to be cultivar-specific genes. The remaining 1,338 genes (4.3%) existed in different numbers of cultivars, ranging between two and nine cultivars, which were considered the dispensable set of genes. Interestingly, the gene content of the cultivars Taishanhong and Tunisia represents 98.9% (30,667) of the whole pan-genome genes, while the cultivar Purpursid added another 198 (0.64%) new genes to the set. The remaining eight cultivars combined added only 151 (0.49%) genes to the gene content of the reference pan-genome ( Figure 2B ).

Our re-sequenced panel was added to the 26 cultivars that were previously used by Luo et al. (2020), which included eight hard-seeded cultivars originated from China, to have a deeper understanding of the domestication of pomegranates and the diversion between soft- and hard-seeded subpopulations. Like Luo et al. (2020), hard-seeded cultivars seem to diverge from semi-soft and soft-seeded cultivars. However, the soft-seeded cultivar “Purpursid” from our re-sequenced panel was clustered with the re-sequenced hard-seeded panel (PAN) in the PCA, phylogeny, and ADMIXTURE analyses ( Figures 3A–C ). The PAN panel in the present study (all hard-seeded except Purpursid) formed an independent cluster that is most related to the hard-seeded cultivars of Luo et al. (2020). The hard-seeded cultivars from the PAN panel showed the highest heterozygosity over all other groups ( Figure 3D ) and adding them to the hard-seeded subpopulation had a substantial influence on reducing LD decay, especially within a distance of <1 Mbp ( Figure 3E ). The Azerbaijani cultivar Fatima and the Uzbekistani cultivar Achygdona showed some admixture with the Chinese hard-seeded cultivars, while the Chinese cultivar Pom23 showed some admixture with the re-sequenced group. Two soft and semi-soft Italian cultivars (Pom14 and Pom15) and one soft cultivar from the USA (Pom6) that belong to one of the soft-seeded subpopulations (colored blue in Figure 3C ) showed some admixture with the Chinese hard-seeded cultivars.

Having a soft-seeded cultivar with a genetic background related to the hard-seeded cultivars in the present dataset, “Purpursid,” increased the power of detecting genes that control seed hardness and reduced the confounding effect of population structure. Our analysis detected 10 genomic regions that are highly differentiated between the soft- and hard-seeded subpopulations, of which four were common with the 24 regions detected in Luo et al. (2020). The description of these genomic regions and the annotations of the closest genes to them are presented in Table 2 . Of these, chromosome 5 seems to have the largest region between 11.36 Mbp and 11.50 Mbp, which involves a cluster of four candidate genes. All SNPs within this region had R² values >0.95. Many of the detected genes along the genome have dehydrogenase and CoA synthase functions, as well as other seed development and cell wall-related activities.