For the present investigation, young leaf samples from the only persisting mature tree specimen of Acacia pachyceras protected in the Sabah Al Nature Reserve (SANR) of Kuwait ( Figure 1 ) was collected during January 2020. Three to four seed pods from this tree were also gathered, their embryo excised and germinated at KISR laboratories, out of which only one survived. Young leaves were snipped from the tip of this one year old seedling for DNA isolation. Leaf samples from morphologically distinct Acacia genotypes growing in close vicinity of the mother tree (north west, n=5; and north east, n=5) in the SANR were also collected for genetic diversity analysis. Approximately 2-3 months old seedlings of Acacia gerrardii (Suleiman, 2017; Suleiman et al., 2018) growing (n=5) under the green house of Kuwait Institute for Scientific Research (KISR) were also sampled for comparison with field grown accessions. GPS coordinates of all the samples were recorded ( Table S1 ). In total 17 samples were collected and DNA was isolated from the leaf samples as per the CTAB protocol described elsewhere (Habibi et al., 2022b). DNA quality, purity and quantity were estimated through spectrophotometric (Nanodrop UV/Vis spectrophotometer, ThermoFisher) and fluorometric (Qubit, Invitrogen, WA) methods as per the manufacturer’s instructions.

Plant tissues from the standard Panicum miliaceum (2n=4x=36; 2C=2.09 pg DNA) (Kubešová et al., 2010) A. pachyceras, and A. gerrardii were chopped into small pieces ranging between 0.3 to 1.0 mm into a Petri dish filled with a mix of OttoI (0.1 M citric acid + 0.50% Tween20) and OttoII (0.40M of Na₂HPO₄.12H₂O) (1:1). Whole mixture was filtered through 42-μm nylon mesh and centrifuged at 5000 rpm and incubated at 40°C. The pellet was collected and resuspended in OttoI. The suspensions were subjected to RNA degradation using RNase (Sigma Aldrich) for 30 minutes. Nuclear DNA was stained by adding 75 µl of propidium iodide (Sigma Aldrich). Flowcytometry based analysis of stained nuclei was performed on 3 laser flow cytometer MACSquant10 analyzer (Miltenyi Biotech GMBH). A total of 100,000 events were recorded for each sample and florescence was measured in B2 channel and median florescent intensities were measured. Final data was normalized and analyzed by FlowJo™ v10.8 Software (BD Life Sciences). The absolute DNA amount of a sample was calculated based on average primary peak value as per the method defined by (Doležel and Bartoš, 2005).

The single DNA sample from the mature tree of A. pachyceras was lyophilized and shipped to the Beijing Genomics Institute (BGI), Hong kong for paired-end sequencing. Post-quality checks at BGI, the DNA was tagmented and subjected to library preparation. DNA fragments were end repaired, 3’A tailed and ligated with Illumina adapters (Khan et al., 2016). Thereafter the libraries were purified through Agencourt AMPure XP magnetic beads (Beckman Coulter Genomics, Brea, CA, USA) (Habibi et al., 2021). Sequencing was performed on the Illumina HiSeq 2500 (Illumina, San Diego, US) platform using the 2 x 250 cycle chemistry. Quality parameters for the raw data were accessed through FASTQC version 0.119 (Andrews, 2010). The raw reads were trimmed using the Trimmomatic v. 0.17 (Bolger et al., 2014). The filtered high-quality sequences were initially assembled by Platanus-allee 2.0 (Kajitani et al., 2019). JELLYFISH 2.1.4 was used to conduct the K-mer analysis (Marçais and Kingsford, 2011). Genome size estimation was done based on k-mer frequency distribution and the number of 17-mer.

For gene annotation the high-quality paired-end DNA sequencing data was used for de-novo assembly of Acacia pachyceras genome using MaSuRCA-4.0.3 (Zimin and Salzberg, 2022). The primary assembly was filtered to remove scaffolds shorter than 500 bp. A de-novo repeat library for filtered assembly was constructed using RepeatModeler2.0.3 (Flynn et al., 2020) that included RECON (Bao and Eddy, 2002), RepeatScout (Price et al., 2005), and Tandem Repeat Finder (Benson, 1999). The repeat library was then subjected to RepeatMasker 4.1.3 to find and mask the repeats in the assembled genome using rmblast as the default search engine (Smit and Hubley, 2008). The BRAKER2 pipeline was used to perform gene prediction by integrating ab initio gene prediction, based on plant protein sequences, combining both GeneMark-EP and AUGUSTUS 3.2.0 (Brůna et al., 2021). The plant protein sequences from various species used for gene prediction were downloaded from the OrthoDB v10 database. The ProtHint 2.5.0 protein mapping pipeline was used for generating required hints from the plant protein sequences for BRAKER (Hoff et al., 2018). The assembled scaffolds and generated hints from the protein sequences were used for obtaining initial gene structures using the GeneMark-ET tool. The initial gene structures were then used for training by AUGUSTUS to produce the final gene predictions. The predicted genes were filtered to retain those with fully supported hints. GC content of the genome sequence was estimated (Gurevich et al., 2013). The BUSCO v3.0.2 (Benchmarking Universal Single-Copy Orthologs) tool was used to evaluate the completeness of the genome (Simão et al., 2015).

The GMATA 2.0 tool was used to filter microsatellite regions in the genome (Wang and Wang, 2016). The recognition criteria for the number of di-, tri-, tetra-, penta- and hexa-nucleotides repeats were set as 5 and above. The in-built Primer 3 (v 3.0) in GMATA was used to design primers against the filtered SSR motifs (Untergasser et al., 2012). The parameter setting of primer design was 18 ~ 23 bp primer size, annealing temperature at 57–62 °C, GC content at 30–70% and 100 ~ 400 bp final product length. A total of 110 primers were randomly chosen and synthesized for PCR amplification (Thermo Fisher Scientific, Waltham, MA, USA). All the primers were dissolved in nuclease free water and adjusted to 1 μM working stock. PCR was conducted on 17 genotypes of Acacia. A reaction volume of 20μL was assembled in a 200 μL PCR plate comprising 10 ng of DNA, 1× PCR Master Mix (Solis BioDyne, Tartu, Estonia), and 0.3 μM of each forward and reverse primers. Thermal cycling was performed in a Veriti Thermal Cycler (Applied Biosystems, Grand Island, NY) by initially inactivating the DNA polymerase at 95°C for 12 mins, followed by 30 cycles of denaturation at 95°C for 30 sec, annealing at 55°C for 30 s and extension at 72°C for 50 sec (Mustafa et al., 2017). The PCR products were visualized on 2.0% agarose gel at 10 V/cm for 60 mins. Gel was visualized on a gel documentation system (Chemidoc MP, BioRad, USA) under UV light. All the samples amplified and produced clear, reproducible bands generating 186 loci, which were scored as present (1) or absent (0). A binary matrix (1/0) was generated and analyzed through the GenAlEx 6.5 software (Peakall and Smouse, 2006). Population structure was determined by the neighbor-net split decomposition network generated by SPLITS Tree 5.0 and bootstrapping runs of 1000 replicates (Kloepper and Huson, 2008) and principal coordinate analysis (PCoA).

Two 300-400 bp (insert size) paired-end libraries yielded 108 GB, raw reads. Per base quality score ranged between 34.0-37.0 (mean=36.0). Approximately 10% of low-quality data was filtered and the remaining 96 GB (ca. 89%) was used for downstream analysis. Among these about, 645 Mb qualified for further analysis. K-mer analysis provides an estimation of genome size based on substrings of length k contained within a biological sequence. Besides, this it also indicates low quality or contamination in the sequences. The 17-mer frequency analysis identified the main peak at 120 X depth based on the k-mer number of 86,462,013,672 ( Figure 2 ). Upon dividing the k-mer number by 120 the genome size was predicted as ~720 Mb (720,516,781bp). The other peak at ½ of the main peak (50X) is most likely due to the high level of heterozygosity (6.84%), whereas an additional peak at a similar depth (25X) is because of repeat rates at the position of multiple integers (56.89%).

A total of 640 million paired-end reads were used to construct a raw genome assembly of A. pachyceras. The assembled genome consisted of 169,210 scaffolds (1.19 Gb). Removal of scaffolds < 500 Mb yielded a filtered assembly of ~1.043 Gb containing 155,442 scaffolds ( Figure 3A ). The number of sequences excluding Ns was 1,042,649,449 bp. The final assembly contained very low N content (~0.07%). The largest scaffold was 120,115 bp with an N50 of 12,073, N90 of 3,238, L50 of 24,820 and L90 of 88,523. The genome contained more adenosine bases (340,306,430) followed by thymine (338,837,198) > cytosine (181,783,987) and guanine (181,721,834) ( Figure 3A ). The GC content was predicted to be 34.84% employing a sliding window algorithm at every 10kb sequence. ( Figure 3B ). BUSCO assessment of completeness of genome revealed 93.2% of assembly to be complete (C: 2,168). There were 1,127 (48.5%) complete and duplicate BUSCOs (D), 1,041 (44.8%) complete and single copy BUSCOs (S), and a mere of 97 (4.2%) of fragmented BUSCOs (F). Total 2,368 BUSCOs were searched among which only 2.5% (67) were missed (M) ( Figure 3C ). The higher number of C as compared to F and M BUSCOs is suggestive of a high-quality genome assembly.

Generic repetitive DNA often assist in expression of typical coding sequences and important in organizing functions required for genome replication and progression to the following generation of cells. Therefore, we characterized the repetitive sequences in the assembled genome and found its total length as ~488.3 Mb, accounting for ~46.8% of the draft assembly. Among these, ~34% of the repeats were unclassified. DNA transposons corresponded to ~2%, whereas, retroelements corresponded to ~9% of the genome. A complete list of repeats along with the content in the draft genome has been shown in Table 1 . Total interspersed repeats were 45.42%.

Gene prediction using BRAKER2 pipeline based on ab-initio method, and ortholog protein sequence- resulted in a total of 123,610 genes corresponding to 127,367 transcripts. However, these include both supported and unsupported transcripts. The transcripts that are fully supported by hints are the most confident predictions. Hence, we used a script, ‘selectSupportedSubsets.py’ provided by the BRAKER2 pipeline to extract the most confident predictions from the complete list of transcripts. This resulted in a total of 34,374 transcripts corresponding to 33,650 genes. These set of transcripts were considered for further analysis. The statistics of most confident transcript and protein sequences are shown in Table 2 . The mean length of the coding sequence was 1,027 nts (Max 11,262 nts; Min-201nts), whereas that of protein sequence was 342 aa (Max-3,754 aa; Min-67aa).

Mining the assembled genome for microsatellite motifs in GMATA yielded 901,755 SSR motifs. Among these, dinucleotides were the largest in number (796,441; 88%), followed by the trinucleotides (90,769; 10.06%), tetranucleotides (12,435; 1.3%), pentanucleotides (1,651; 0.18%), and, hexanucleotides (418; 0.05%) ( Figure 4A ). In total 11,596 primer pairs were designed against these SSR motifs and 11,181 among them were unique ( Figure 4B ). Approximately 98.71% of SSR motifs were positioned in the genome that did not meet the standard primer designing parameters hence no primers could be designed for them. The GBS based SNP markers are recommended to be investigated to overcome this limitation.

A subset of 110 primer pairs representing each repeat class were randomly chosen for PCR validation ( Table S2 ). Amplification yielded 186 genetic loci for downstream analysis. Among, these 141 (75.65%) were polymorphic (PIC=0.27). These loci were imported in the GenAlEx v6.5 software for estimating the pairwise genetic distances, analysis of molecular variance (AMOVA) and principal coordinate analysis. The PCoA analysis revealed the Acaciaaccessions to be distributed into three clusters. Variations at 1^st 2^nd and 3^rd axis were 24.26%, 38.54% and 49.06%, respectively. The PCoA distributed the Acacia genotypes into a small and a large cluster ( Figure 5A ). The smaller cluster had only two genotypes AP01(single specimen of A. pachyceras growing in SANR) and AP12 (seedling of the mother tree). Rest all the genotypes formed a part of the larger cluster. The UNK-NE, UNK-NW and the AG-S existed as three overlapping sub-clusters within the larger cluster. Genotype AP-17 admixed with UNK-NW accessions whereas AP3, AP4, AP11 and AP16 formed a separate group. Based on the grouping of UNK-NW and UNK-NE accessions with AG-S in the larger cluster we hypothesize these morphologically distinct genotypes growing near the single specimen of A. pachyceras in SANR area to be A. gerrardii. Pairwise genetic distance (Nie’s D) between the genotypes ranged from 6 (AP09-AP10) to 77 (AP07-AP12) ( Table S3 ). We constructed a split decomposition tree on this genetic distance and observed AP01 and AP12 to be closely related. The grouping of other genotypes was in agreement with the PCoA analysis ( Figure 5B ).

Analysis of molecular variance (AMOVA) partitioned the genetic diversity within populations to be 75% whereas, among population diversity was 25% with a PhiPT (r) equalling to 0.249 (p<0.001). Total 999 permutation events were performed on 17 samples and 189 loci. Pairwise genetic distance within populations were 47 (AP-MT), 34.2 (UNK-NW), 25.4 (UNK-NE) and 30 (AG-S).

The genome size of A. pachyceras was also estimated through flow cytometry with reference to the standard P. milaceum ( Table 3 ). The 2C content of A. pachyceras was calculated as 2.46 pg DNA with a 1C value of 1.23 pg DNA. The holoploid genome size was estimated as 1204.9 Mb by applying a factor of 1 pg of DNA equivalent to 978 megabases. As compared to the standards’ 2n peak, A. pachyceras showed three distinct peaks suggestive of polyploidy in the collected specimen. The first peak or G1 peak is the diploid peak, whereas the second and the third peaks are most likely due to tetraploidy and hexaploidy. The median florescent intensity (MFI) of primary peak (G1 peak) had slightly higher value of 7.38 than the reference peak (6.26) whereas the secondary and tertiary peaks for A. pachyceras exhibited MFI equivalent to 13.20 (1.7 times higher) and 24 (3.2 times higher) respectively ( Figures 6A, B ). The difference in peak heights is attributed to different cells exhibiting varied ploidy levels in the same plant tissue. Based on the ploidy level analysis the 1Cx pg of DNA (monoploid DNA content) was estimated as 0.41 pg suggesting a monoploid genome size of A. pachyceras as 400.98 Mbp. Ploidy and genome size of A. gerrardii were also analyzed through flow cytometry. The 2C, 1C and 1Cx values were in parallel with A. pachyceras ( Figures 6C, D ). This species also possesses a hexaploid genome alike A. pachyceras.