Genomic insights into the conservation status of the world's last remaining Sumatran rhinoceros populations

Abstract

Small populations are often exposed to high inbreeding and mutational load that can increase the risk of extinction. The Sumatran rhinoceros was widespread in Southeast Asia, but is now restricted to small and isolated populations on Sumatra and Borneo, and most likely extinct on the Malay Peninsula. Here, we analyse 5 historical and 16 modern genomes from these populations to investigate the genomic consequences of the recent decline, such as increased inbreeding and mutational load. We find that the Malay Peninsula population experienced increased inbreeding shortly before extirpation, which possibly was accompanied by purging. The populations on Sumatra and Borneo instead show low inbreeding, but high mutational load. The currently small population sizes may thus in the near future lead to inbreeding depression. Moreover, we find little evidence for differences in local adaptation among populations, suggesting that future inbreeding depression could potentially be mitigated by assisted gene flow among populations.

Fig. 1. Sampling, current range, and phylogeny for Sumatran rhinoceros (Dicerorhinus sumatrensis).

a Geographical origin of the 21 Sumatran rhinoceros genomes analysed in this study. Current distribution is depicted in green. Triangles and circles represent approximate geographical sampling locations of historical and modern samples, respectively. The numbers within the geometric symbols depict the number of samples from the corresponding location. b Unrooted phylogeny (100 bootstrap replicates). Asterisks denote low coverage (< 9X) genomes.

Fig. 2. Population history and timing of population divergence for Sumatran rhinoceros (Dicerorhinus sumatrensis).

Thick coloured lines depict temporal fluctuations in effective population size (N_e) for the three populations with each coloured line representing one individual from each population. Dotted/dashed black curves represent the pseudodiploid sex chromosomes (i.e., X chromosome) used to infer pairwise population divergence times, with the curves going to infinity at the respective time of divergence. The X-axis corresponds to time before present in years on a log scale, assuming a substitution rate (μ) of 2.34 × 10⁻⁸ substitutions/site/generation^, and a generation time (g) of 12 years. The y-axis corresponds to N_e. The grey rectangle depicts the last glaciation. Bootstrap tests were conducted with 100 replicates (Supplementary Fig. 4).

Fig. 3. Inbreeding, relative mutational load and number of loss-of-function (LoF) variants in three modern populations of Sumatran rhinoceros (Dicerorhinus sumatrensis).

a Inbreeding estimated as the average proportion of the genome in runs of homozygosity (F_ROH). Open bars show total proportion of the genome in ROH ≥ 100 kb and solid bars show proportions in ROH of length ≥ 2 Mb. Bars extending from the mean values represent the standard deviation (two-sided pairwise t-test, F_ROH ≥ 100 kb: p_{Borneo-MalayP} = 4.2e−07, p_{Borneo-Sumatra} = 0.0066, p_{MalayP-Sumatra} = 6.4e−07, F_ROH ≥ 2 Mb: p_{Borneo-MalayP} = 2.2e−05, p_{Borneo-Sumatra} = 0.15, p_{MalayP-Sumatra} = 2.2e−05). b Mutational load estimated with GERP scores. Individual relative mutational load was measured as the sum of all derived alleles multiplied by their conservation-score over the total number of derived alleles (see ‘Methods’). Only derived alleles above conservation-score of 1 (i.e., non-neutral) were included (two-sided pairwise t-test, p_{Borneo-MalayP} = 0.00013, p_{Borneo-Sumatra} = 0.001, p_{MalayP-Sumatra} = 0.012. c Number of LoF variants using an annotation of 33,026 gene predictions for a white rhinoceros genome assembly (see ‘Methods’) (two-sided pairwise t-test, p_{Borneo-MalayP} = 2.2e−05, p_{Borneo-Sumatra} = 0.0004, p_{MalayP-Sumatra} = 0.0019). Middle thick lines within boxplots and bounds of boxes represent mean and standard deviation, respectively. Vertical lines represent minima and maxima. n = 14, ***p < 0.001, **p < 0.01, *p < 0.05, ns = non-significant, p-values were not adjusted for multiple comparisons.

Fig. 4. Temporal changes in inbreeding, relative mutational load and number of loss-of-function (LoF) variants in two populations of Sumatran rhinoceros (Dicerorhinus sumatrensis).

a Inbreeding estimated as the average proportion of the genome in runs of homozygosity (F_ROH). Open bars show the total proportion of the genome in ROH ≥100 kb and solid bars show proportions in ROH of length ≥ 2 Mb. Bars extending from the mean values represent the standard deviation (two-sided t-test, F_ROH ≥ 100 kb: p = 0.034, F_ROH ≥ 2 Mb: p = 0.007). b Mutational load estimated with GERP scores. Individual relative mutational load was measured as the sum of all derived alleles multiplied by their conservation-score over the total number of derived alleles (see ‘Methods’). Only derived alleles above conservation-score of 1 (i.e., non-neutral) were included (two-sided t-test, p = 0.047). c Number of LoF variants using an annotation of 33,026 genes for white rhinoceros (see ‘Methods’) (two-sided t-test, p = 0.35). Middle thick lines within boxplots and bounds of boxes represent mean and standard deviation, respectively. Vertical lines represent minima and maxima. n = 6, **p < 0.01, *p < 0.05, ns = non-significant, p-values were not adjusted for multiple comparisons.