Genetic erosion in a tropical tree species demonstrates the need to conserve wide-ranging germplasm amid extreme habitat fragmentation

Palaquium obovatum is a species in Sapotaceae, a diverse family of tropical trees that includes species harvested for timber and latex. Representatives of Palaquium, and in fact most of the Sapotaceae family typically register low abundance in Malayan rainforest (LaFrankie 1996; Turner et al. 1997). Palaquium obovatum is native to the lowland rainforests of Southeast Asia, including Indonesia, Thailand, Malaysia and Singapore, and is sometimes found along coasts. In Singapore, its limited range and low density has rendered it vulnerable, with less than 1000 mature individuals in the wild (Lindsay et al. 2022). This distribution makes the species an excellent candidate for investigating genetic responses to extreme habitat fragmentation in an urban landscape, which are relevant to the preservation of remnant tree populations growing in isolated fragments in many other cities, especially against the backdrop of widespread deforestation in Asia making way for man-made habitats. P. obovatum is a medium to tall tree, with a cylindrical, buttressed bole reaching heights of up to 40 m, and diameter at breast height (dbh) of up to 110 cm. The species has perfect flowers exhibiting characteristics of bat pollination (Phang 2023), which also makes it a pertinent study focus in light of anthropogenic disruptions to bat-plant interactions across the world (Aziz et al. 2021). The fruits are fleshy berries, 2–3 cm in diameter, with 1–2 seeds. Small trees (1–2 cm dbh) grew at a mean rate of 0.04 cm in dbh per year from 1993 to 2004 in a 2-ha primary forest plot within Singapore’s Bukit Timah Nature Reserve (LaFrankie et al. 2005). In Peninsular Malaysia, some congeneric species have been observed to have a similarly slow rate of growth (Soerianegara and Lemmens 1993): trees of Palaquium maingayi and P. rostratum were found to take 100 and 70 years respectively to attain a dbh of 55 cm, and P. gutta a dbh of 50 cm in 50 years; however, a tree of P. rostratum in an arboretum attained a dbh of 57 cm in 40 years. A heritage tree of Palaquium obovatum that was planted in 1897 by Henry Ridley in a part of the Singapore Botanic Gardens used to trial plants of economic potential, has a dbh of 131 cm and height of 31.4 m (National Parks Board 2021); however, growth rates in cultivated areas often differ from those in natural forest (Appanah and Weinland 1993). As a chiroptephilous species, Palaquium obovatum is reliant on bats for both pollination and seed dispersal; when fallen, the ripe fruit are quickly consumed by other fauna such as monkeys, squirrels or birds, which curtails seed survivability (Ng 1972; Soerianegara and Lemmens 1993; Chan et al. 2022; Phang 2023). Seedlings that germinate are often outcompeted due to overshadowing vegetation constraining sufficient light (Soerianegara and Lemmens 1993). The genus, in general, does not regularly flower or fruit but P. obovatum has been observed to flower more often than its congeners (c. every 1–2 years), which increases the accessibility to juvenile material for population studies. The smallest fruit-bearing individual recorded had a dbh of 11 cm (pers. obs).

In Singapore, the species has been assessed as Vulnerable (VU; Lindsay et al. 2022) and has some economic importance as a source of gutta percha—a latex used in a myriad of products including the first undersea telegraphic cables. Nothing is currently known about the species’ genetic diversity, levels of inbreeding and levels of gene flow between populations. Polyploidy in the Sapotaceae family is rare and occurs mostly in the tribe Sideroxyleae, with previous studies finding most representatives, including certain Palaquium species, to be diploid (Johnson 1991; Smedmark and Anderberg 2007).

Individuals of Palaquium obovatum were sampled from populations across their natural range in Singapore. Existing literature and past surveys were consulted (Wong et al. 1994; LaFrankie et al. 2005; Yee et al. 2011, 2019; Neo et al. 2013; Hung et al. 2017; Ngo et al. 2016), and forest reserves were searched to locate the largest extant trees of P. obovatum, which likely survived the widespread deforestation that occurred during British colonisation in the 19th century (Corlett 2013). Sampling included individuals in various life stages, from seedling to reproductively capable adults.

In Singapore, the study area encompassed primary and secondary forest areas as well as presumed cultivated populations (i.e. without historical planting records but found in built-up areas) in the Singapore Botanic Gardens (SBG) and St. John’s Island (SJI). Sampling occurred in the remnant primary forest patches within SBG (which contains a 4 ha area of lowland rain forest that formed part of the original establishment in 1959, known as the “Gardens’ Jungle”; Turner et al. 1996), Bukit Timah Nature Reserve (BTNR) and MacRitchie Reservoir area in the Central Catchment Nature Reserve (CCNR). Secondary forest samples were taken from the Upper Seletar area in CCNR, Bukit Batok Nature Park (BBNP), as well as small islands off mainland Singapore including Sentosa, St. John’s Island and Pulau Ubin. From Peninsular Malaysia, one cultivated population was sampled from Kepong Botanic Gardens at the Forest Research Institute Malaysia, and one herbarium sample from a forest in Langkawi (a district in the state of Kedah) was also included (Fig. 1, coordinates plotted using QGIS: QGIS Development Team 2019). In all, 15 populations were sampled, four from primary forest, eight from secondary forest and three from cultivation (Table 1).

Other than the single herbarium sample from Langkawi, fresh leaf samples from 103 individuals were collected and stored in silica gel. The sampling consisted of 16 adults and of seedlings or saplings collected mostly within a 5 m radius beneath the canopy of these adults (except for the population in BTNR). A total of 18 samples were classified as saplings (> 30 cm in height, dbh < 10 cm) and 70 as seedlings (< 30 cm in height). No seedlings were found in BTNR, and instead nine representative saplings were collected. The dbh of the adult trees was measured and their height was estimated. The seedlings and saplings were grouped into the above size categories, and combined in downstream analyses as juveniles. GPS coordinates and herbarium vouchers were obtained for each adult tree. For adult trees, leaves were obtained either using a pole pruner for branches less than 3 m high, a Big Shot Catapult to launch a weighted line over high branches, or arborist climbers/cranes for sampling in the Singapore Botanic Gardens.

The Singapore samples were collected from eight distinct localities (BTNR, BBNP, CCNR-MacRitchie, CCNR-Seletar, SBG, Sentosa, SJI, Pulau Ubin) and the Malaysia samples from two localities (Langkawi, Kepong). Three populations were from areas of cultivation: one from SBG, one from SJI and the last from Kepong. Trees around or exceeding 100 cm in dbh, i.e. from the Gardens’ Jungle (primary forest), Pulau Ubin, Sentosa and CCNR-Seletar (secondary forest) are likely to be close to or more than 100 years old. Assuming 20–30 years per generation, the oldest trees from wild areas other than BBNP and SJI would have been alive for at least three overlapping generations. The smaller trees from the wild areas of SJI and BBNP, as well as the parent from Kepong Botanic Garden, are unlikely to have produced more than one generation; the larger cultivated instances in SBG and SJI, not more than two generations. There were no seedlings found in the entire 2 ha primary forest plot within BTNR; the sole reproductively capable adult sampled was a young tree only 14 cm in dbh, which is unlikely to have produced the saplings of heights between 1 and 8 m found within the same plot.

Genomic DNA was extracted from fresh leaf material using the cetyl-tri-methylammonium bromide (CTAB) method, with 24:1 chloroform: isoamyl alcohol and overnight precipitation in isopropanol at -20 °C (Doyle and Doyle 1987; Cullings 1992). Samples were purified with SeraMag magnetic carboxylate-modified microparticles (Cytiva, USA) diluted to 5%, and the concentration of DNA quantified with a Nanodrop 2000c Spectrophotometer (Thermo Scientific, Singapore). Double-digest Restricted Associated DNA (ddRAD) libraries of 16–99 individuals per library (pools depended on collection and lab access timings) were built using protocols from Peterson et al. (2012) and adapted by Niissalo et al. (2018; 2020), with DNA digested using PstI (rare-cutting) and ApeKI (frequent-cutting) restriction enzymes, and standard Illumina indexing barcodes of 5 bp, adapters and polymerase chain reaction (PCR) primers (Integrated DNA Technologies, Singapore). DNA concentration of 400 ng per sample was used, and SeraMag particles applied to size select 450–550 bp fragments, with quality and size confirmed with a Tapestation 4200 System using High Sensitivity D1000 ScreenTapes (Agilent Technologies, CA, USA). To enrich digested fragments, 12 PCR cycles were conducted for each library, with annealing temperature of 64 °C. Libraries were analysed with a Qubit analyser (Life Technologies Corporation, Carlsbad, CA), pooled to equimolar amounts and sequenced by NovogeneAIT Genomics (Singapore) using an Illumina Novaseq 6000 instrument (San Diego, CA, USA).

Raw reads were quality-checked with FastQC v.0.11.9 (https://​www.​bioinformatics.​babraham.​ac.​uk/​projects/​fastqc). Sequences were then filtered to eliminate poor quality reads, demultiplexed, and reads assigned to individuals by specifying the parameters -q -c -r in the process_radtags program within STACKS 2.53 (Catchen et al. 2013). KMC v. 3.1.0 (Kokot et al. 2017) was used to calculate k-mers, and heterozygous k-mer pairs to assess ploidy levels of demultiplexed reads were estimated using Smudgeplot v. 0.2.3dev (Ranallo-Benavidez et al. 2020). A graph of the percentage of the minor and major common variants was plotted against the kmer count; a peak near 0.5 is expected if the individual is diploid, and 0.167 if hexaploid.

A reference genome for Palaquium obovatum, assembled from short reads, was obtained from the ongoing project to assemble whole genomes of all native species found in BTNR (Niissalo et al., in prep; N50 = 2,381 bp; BUSCO = 64.0%; assembly size = 1.4 Gb: the sizes of other Palaquium species sequenced in the same project were also referenced), and indexed with bwa v.0.17.17 (Burrows-Wheeler Aligner; Li and Durbin 2009; https://​github.​com/​lh3/​bwa), Reads for each individual were mapped to the reference genome with the bwa-mem setting in bwa v.0.7.17, then sorted with the sort function in Samtools v.1.14 (Li et al. 2009) to obtain BAM files per sample, which were evaluated by using the samtools stats and depth settings to generate mapping statistics and coverage.

To identify genetic variation, GATK v.4.2.6.1 (https://​github.​com/​broadinstitute/​gatk) was used, along with Picard v.2.19 (https://​github.​com/​broadinstitute/​picard) and Samtools v.1.14 (https://​github.​com/​samtools/​) to sort, mark duplicates, and index the data, followed by HaplotypeCaller with ploidy set at the estimated level (according to Smudgeplot results) to call SNPs on individual samples to create GVCF files with alignment to the most likely haplotypes. These GVCF files were consolidated using GenomicDBimport and jointly genotyped with GenotypeGVCF. Indels and SNPs of low base quality and depth of coverage were removed using the following hard-filtering parameters in GATK: “QD < 2.0 || QUAL < 30.0 || SOR > 3.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0”. To minimise the influence of singletons, further filtering was then done with bcftools v.1.14 (https://​github.​com/​samtools/​bcftools) to remove multiallelic SNPs, those below a set frequency (maf 0.01) and all with missing data.

Genetic deterioration through the accumulation of deleterious genes is normally detected in measures of excess homozygosity and inbreeding, which eventually manifests in reduced fitness of organisms (Bosse and van Loon 2022). Heterozygosity was assessed by generating related metrics (i.e., distribution and frequency of heterozygous and homozygous genotypes in the dataset) with the stats -s parameter in bcftools. A custom script was then applied to only retain common alleles present as one or two copies out of a total of six in the whole sample set, effectively converting the vcf into one with diploid calls allowing direct comparison of the same genetic parameters across all samples (Supplementary File 1). To address linkage disequilibrium between SNP pairs, the dataset for unlinked SNPs was generated using the pruning parameters --indep-pairwise 50 5 0.5 (i.e. window size of 50 variants, window shift of 5 variants and r² of 0.5) in plink v.1.90 (Chang et al. 2015). This pruned vcf file was examined for genetic variation by generating heterozygotic and inbreeding sample-specific statistics using the --het parameter in plink v.1.90 (Chang et al. 2015), as well as nucleotide diversity (π) and Tajima’s D averaged across loci using parameters --window-pi and --TajimaD in vcftools v.0.1.16, with reference to Tajima (1989) for critical values indicating deviation from mutation-shift balance. Paired t-tests were run to determine whether differences in heterozygotic and inbreeding measures were significant between cohorts of age classes and type of sampling locality (i.e. between wild and cultivated populations, as well as between primary and secondary forest populations).

To assess genetic similarity, a pairwise comparison of heterozygous single nucleotide polymorphisms (SNPs) was performed using the script from https://​github.​com/​mattiniissalo/​Genepop_​clonality_​script but adapted by manually converting the STRUCTURE file to genepop format (typically used by the Stacks pipeline). Specifically, each heterozygous SNP present in a sample was compared to its corresponding locus in another sample. A heatmap was also generated for these pairwise comparisons using the R package ggplot2 v.3.4.0 (Wickham 2011) in Rstudio v.4.1.3. After converting the pruned vcf file to STRUCTURE format with PGDSpider2 v1.1.5 (Lischer and Excoffier 2012), STRUCTURE (Pritchard et al. 2000) and Structure_threader v1.3.10 (Pina-Martins et al. 2017) were used to infer population structure, with K values (i.e. number of genetic clusters) ranging from one to 15, using 10 replicates for each value of K, with a burn-in of 20,000 and Markov Chain Monte Carlo (MCMC) repetitions at 100,000. The optimal value of K was estimated with the software Kfinder2 v.1.0.0.0 (Wang, 2019), which generates estimators of the best K using the ΔK (Evanno et al. 2005) and Pr[X|K] (Pritchard et al. 2000) methods, which take into consideration the rate of change in log probability between consecutive values of K (ΔK). The R package pophelper v.2.3.1 (Francis 2017) was applied to produce barplots for each value of K.

Plink was used to convert the consolidated VCF file to eigenvector/eigenvalue format for principal components analysis (PCA), with the first eight dimensions visualised using the tidyverse v.1.3.2 package (Wickham et al. 2019) in Rstudio v.4.1.3. The association between genetic and geographical distance was assessed by a Mantel test with 9999 permutations performed with the ade4 v.1.7–20 package (Dray and Dufour 2007) in Rstudio v.4.1.3. To further assess population differentiation, weighted F_ST (Bhatia et al. 2013) was calculated using vcftools v 0.1.14 (Danecek et al. 2011). Gene flow amongst populations was investigated using the ABBA-BABA test in D-suite (Malinsky et al. 2021); the sample from Langkawi was fixed as the outgroup (it had emerged as one of the earliest diverging samples in the Palaquium obovatum clade from a previous phylogenetic study, Phang et al. 2023), and the dataset was split into 20 jackknife blocks.