Genomic basis for skin phenotype and cold adaptation in the extinct Steller's sea cow.

Steller's sea cow, an extinct sirenian and one of the largest Quaternary mammals, was described by Georg Steller in 1741 and eradicated by humans within 27 years. Here, we complement Steller's descriptions with paleogenomic data from 12 individuals. We identified convergent evolution between Steller's sea cow and cetaceans but not extant sirenians, suggesting a role of several genes in adaptation to cold aquatic (or marine) environments. Among these are inactivations of lipoxygenase genes, which in humans and mouse models cause ichthyosis, a skin disease characterized by a thick, hyperkeratotic epidermis that recapitulates Steller's sea cows' reportedly bark-like skin. We also found that Steller's sea cows' abundance was continuously declining for tens of thousands of years before their description, implying that environmental changes also contributed to their extinction.

INTRODUCTION

During the Pleistocene, our planet was home to a broad diversity of megafaunal species, including mastodons, mammoths, giant camels, bears, and several large felids. Many of these disappeared suddenly during the Late Pleistocene and Early Holocene. Thus, our knowledge of these species is largely drawn from reconstructions of their ranges and phenotypes based on fossil remains. Several megafauna species, however, went extinct during the later Holocene, often on islands after human colonization. Some of these species survived long enough to be described by European explorers, who provided written details of behaviors and phenotypes that could not have been inferred from their scattered and fragmentary remains. Among the recently extinct megafaunal species is Steller’s sea cow, Hydrodamalis gigas, which was described by Georg Wilhelm Steller in 1741 and became extinct 27 years later, presumably as a result of human-driven habitat change () and overexploitation (). The disappearance of Steller’s sea cow has been regarded as the first historical extinction of a marine mammal as a consequence of human actions () and is an iconic example of the catastrophic impact that human hunting can have on a species (). At the time of Steller’s description, the Steller’s sea cows’ range had already been reduced to shallow waters around the uninhabited Commander Islands, Russia (Fig. 1A), where Steller reported a population of around 1000 individuals ().

Fig. 1.

Molecular basis for Steller’s sea cow’s skin phenotype.

(A) Sirenian distribution according to the International Union for Conservation of Nature Red List (). All sequenced Steller’s sea cow individuals originate from the Commander Islands. (B) Translated multiple sequence alignment of the ALOXE3 and ALOX12B genes showing amino acid sequence conservation corresponding to the human proteins (bold in Steller’s sea cow) and the position of the premature stop codons. (C) Arachidonate lipoxygenases structure, which is composed of the PLAT (Polycystin-1, Lipoxygenase, and Alpha-Toxin) domain and the enzymatic LIPOXYGENASE core domain. Premature stop codons in Steller’s sea cow ALOXE3 and ALOX12B genes are depicted in red. Truncating variants described in human patients and located downstream from the Steller’s sea cow premature stop codons are depicted in black. (D) Left: Steller’s sea cow drawing according to Steller’s description from 1741 (image by R. Ellis). Right: Image of a patient with ichthyosis; detail depicts scaling and hyperkeratosis.

Adult Steller’s sea cows reached a length of approximately 10 m, weighed over 10 metric tons, and stored up to 10 cm of blubber in some areas of the body (). This made them ideal resources for human hunters, who exploited them for meat, fat, and skin. Their skin was thick, hairless, and roughly textured, an attribute that led Steller to describe it as resembling “more the bark of an old oak tree, than the skin of an animal” (). The rough, bark-like skin may have been an adaptation to their shallow water habitat, perhaps preventing abrasion on ice or rocks (). For centuries, naturalists relied solely on phenotypes to describe species, lacking opportunity to explore the genetic underpinnings of those traits. Unlike for many other extinct species, the existence of a thorough phenotypic record for the Steller’s sea cow can facilitate the identification and understanding of the underlying molecular basis of their traits.

RESULTS AND DISCUSSION

Here, we generate paleogenomic data from Steller’s sea cow bones and explore the genetic underpinnings of their unique phenotype and population history. We extracted ancient DNA from fragmentary remains of 12 Steller’s sea cow individuals recovered from beaches of Bering Island, which we radiocarbon dated to 2205 to 1155 before the present (B.P.) (see section SM1, in the Supplementary Materials). We sequenced the genomes of two of the best preserved of these (SNMB N51667 from the Braunschweig Natural History Museum and SC16.JK045 from the collection of the Kamchatka Branch of the Pacific Geographical Institute, Russian Academy of Sciences Museum) to 15.86 and 15.63× coverage, respectively, and the remaining 10 to an average coverage of 2.78× (range 1.97 to 3.95×) (metrics for libraries and mapping quality are presented in table S1). Our analysis of mitochondrial genomes from the 12 Steller’s sea cows confirms that they are different individuals (section SM1). We mapped the recovered data to a new de novo–assembled draft genome of the dugong (Dugong dugon, originating from the Australian Museum; section SM2), which is the closest living relative of Steller’s sea cow (, ). To explore how genomic changes along the lineage to Steller’s sea cow might explain their unique phenotypes, we generated a comparative dataset consisting of 4877 orthologous genes. To this end, we annotated our new dugong genome and a previously published Florida manatee genome [Trichechus manatus latirostris (); University of California, Santa Cruz (UCSC) GCA_000243295.1] using the MAKER pipeline () and annotations from human, mouse, and elephant genomes (). Because pseudogenization can shape phenotypes (, ), we focused first on identification of inactivated genes and found two arachidonate lipoxygenases (ALOXE3 and ALOX12B) that have been inactivated via premature stop codons along the lineage to Steller’s sea cow (section SM3).

In humans, loss-of-function variants in the lipoxygenase genes ALOX12B and ALOXE3 are the second most common cause of autosomal recessive congenital ichthyosis, which is a disease characterized by hyperkeratotic, dry, thickened, scaling skin (). The premature stop codons in Steller’s sea cow are located at position 104 in ALOX12B and position 394 in ALOXE3 (Fig. 1B) and are present in all our 12 sequenced individuals (section SM3), as well as in an additional individual recovered from the same site (). In humans, inactivating variants with functional relevance in both genes have been described downstream of those present in Steller’s sea cow (–) (Fig. 1C) and are absent in healthy controls supporting the proposed functional relevance of the Steller’s sea cow’s pseudogenizations. Skin appearance as described by Steller in 1741 as resembling “the bark of an old oak tree” is phenotypically similar to skin appearance in ALOXE3- or ALOX12B-deficient humans (Fig. 1D) and the respective gene-deficient mouse models (, ).

On the basis of the presence of inactivating mutations in ALOX12B and ALOXE3 in our 12 sequenced individuals, we calculated the Bayesian credible interval for the population allele frequency to be 0.9 to 1 with a posterior probability of >0.95. This suggests that the mutations were either already fixed or in the process of reaching fixation in the population from which we sampled and may have provided an adaptive advantage to Steller’s sea cows.

To explore this possibility further, we scanned all available mammalian sequences including extant marine mammals. We found that both ALOX12B and ALOXE3 are also inactivated in extant cetaceans through gene loss (ALOX12B) and the accumulation of premature stop codons, deletions, or frameshifts (ALOXE3) but are active in pinnipeds (, ), sea otter, polar bear, and extant sirenians (Fig. 2, A and B, and section SM3). The absence of an ALOX12B ortholog in all extant cetaceans suggests that the gene loss happened in the lineage leading to their most recent common ancestor. Because the inactivation of ALOXE3 occurred via different changes in each cetacean lineage (Fig. 2B), this loss must have happened more recently through several independent events.

Fig. 2.

Genomic locus and function of skin lipooxygenases.

(A) Synteny map of the genomic locus around ALOX12B and ALOXE3 in human, mouse, sea lion, dolphin, blue whale, and dugong. The ALOX12B locus is missing in cetaceans. (B) ALOXE3 inactivation occurred through different changes in cetaceans and the Steller’s sea cow; branches of species with ALOXE3 inactivation are depicted in red. (C) Structure of the epidermis. Acyl-ceramides are produced mainly in the stratum granulosum and partly in the stratum spinosum and stored in lamellar bodies as acyl-glucosylceramides. At the transition to stratum corneum, lamellar bodies release their contents into the extracellular space, where acyl-glucosylceramides are converted to acyl-ceramides. Released acyl-ceramides, fatty acids, and cholesterol form lipid lamellae in the stratum corneum. Some acyl-ceramides are oxidized in the presence of skin lipoxygenases, like enzymes encoded by ALOX12B and ALOXE3. The oxygenation is required to facilitate hydrolysis and subsequent covalent linkage of the ω-hydroxyceramides from the corneocyte-bound lipid envelope to proteins of the cornified envelope, which form the epidermis permeability barrier. ALOXE3’s accumulation of inactivating mutations in the absence of ALOX12B (B) confirms its downstream location in the biochemical pathway (C).

Understanding the constraints imposed by this ancient macroevolutionary innovation could provide insight into the mechanism of the human disease. The two encoded proteins are proposed to operate in the metabolic pathway involved in epidermal structural lipid formation (Fig. 2C) (). From an evolutionary perspective, accumulation of different inactivating mutations in ALOXE3, probably as a result of genetic drift in the absence of ALOX12B, supports the proposed involvement of both genes in the same pathway, which becomes nonfunctional in the absence of ALOX12B.

Cetaceans have a high shedding rate of skin’s stratum corneum (), which prevents the buildup of the ichthyotic leathery skin that occurs in the absence of skin lipoxygenases (, ). The exfoliation rate of cetaceans may be due to loss of the desmosome genes DSC1 and DSG4 (), which form the strongest bonds in the outermost skin layer (). We did not identify any inactivating mutations in these desmosome genes in Steller’s sea cow (section SM3), which is likely the reason why Steller’s sea cow, unlike cetaceans, developed the ichthyotic skin appearance, as suggested by Steller’s description.

Gene losses or pseudogenization may occur as an adaptation or because the gene function becomes obsolete and the gene sequence drifts. Previous studies have suggested that the loss of skin lipoxygenases accompanied by hair loss may be an adaptation to a fully aquatic environment (). However, because we find that skin lipoxygenases are convergently inactivated in extant cetaceans and Steller’s sea cow, while the similarly fully aquatic extant sirenians retained active orthologs (Fig. 2, A and B, and section SM3), the marine environment is insufficient to explain this adaptation. A key difference between the environments of Steller’s sea cows and the two extant sirenian species is the temperature of the oceans in which the lineages live. We therefore hypothesize that the loss of lipoxygenase genes may have been adaptive in particular in a cold marine environment and may reflect a different adaptation strategy to retain heat compared to fur-bearing marine mammals.

In further support of convergent evolution in Steller’s sea cow and cetaceans regarding energy metabolism, we identified NPFFR2 as inactivated along both lineages but not in pinnipeds or extant sirenians (section SM3). NPFFR2 encodes for neuropeptide FF receptor-2, a G protein–coupled receptor that has been linked to thermogenesis. Usually, excess caloric intake results in increased fat accumulation and an elevated energy expenditure via diet-induced adaptive thermogenesis. Mice lacking Npfr2, however, develop obesity and fail to activate brown adipose tissue (BAT) thermogenesis [nonshivering thermogenesis (NST)] in response to increased caloric intake (). Selection for large body size is predicted to reduce the importance of BAT for NST by enabling the development of size-dependent heat conservation measures, such as a decrease in relative surface area for heat loss (). In addition to inactivation of NPFFR2, Steller’s sea cows and sequenced cetaceans also lack UCP1 (), which is very important for NST in BAT. In contrast to pinnipeds, which express UCP1 and NPFFR2, UCP1-based NST may therefore not have played an important role in Steller’s sea cow’s adaptation to cold-water habitat, which is similar to cetaceans (). While NPFFR2 absence does not fully explain the increased body size in cetaceans and Steller’s sea cows, it may contribute to their large size by helping to mitigate the cost of heat loss ().

Because evolution does not only proceed by functionally inactivating genes via obvious pseudogenization (nonsense and frame-shifting mutations), we next searched for intact genes evolving at significantly different evolutionary rates in Steller’s sea cow compared to extant sirenians (section SM3). Of the 197 genes evolving significantly faster and 41 evolving significantly slower in Steller’s sea cow, 20 have been shown to modulate body weight and energy metabolism (section SM3). This significant enrichment in genes involved with energy homeostasis and weight regulation (hypergeometric test P value = 0.006) is accompanied by an apparent lack of signal in the gene ontology category related to BAT (GO:0050873; 16 annotated genes in our orthologous set; section SM3), supporting the suggestion that BAT-involving NST indeed did not play a role in thermoregulation in Steller’s sea cows. In addition, several genes with significantly different evolutionary rate on the Steller’s sea cow lineage (ACP6, ACSF3, ACSL5, EHHADH, IVD, PLA2G2A, and PLA2G4F), two of which (ACSL5 and PLA2G4F) are faster evolving in cetaceans, are components of the lipid-degrading metabolic pathways. Some of these are involved in thermogenic futile cycles based on adenosine 5′-triphosphate (ATP)–demanding lipolysis/re-esterification (), which generate heat by ATP hydrolysis and thus contribute to NST (section SM3) (, ).

A shift in the energetic balance toward fat accumulation renders an increased capacity for fasting, which may have been key to both survival and long-distance movements of Steller’s sea cows. Steller’s sea cows lacked dentition (, ), and their diet depended almost exclusively on kelp (). Dragon kelp (Eularia fistulosa), the predominant surface canopy kelp species in the western North Pacific, is present only from about May through September (), such that animals would have entered a fasting state for many months each year. As the Pleistocene fossil record shows that Steller’s sea cows were widespread in shallow Pacific waters from present-day Japan to the Mexican Baja Peninsula (Fig. 1A), fasting may have been critical to dispersing between oceanic islands.

The only published information on historical Steller’s sea cow population size is Steller’s first-hand report from Bering Island during the winter of 1741–1742. That account indicates that Steller’s sea cows numbered in the general range of a thousand animals (). If the pre-Holocene population density of Steller’s sea cows was roughly similar in other kelp forest habitats across the Pacific rim, then this translates to a potential peak abundance of roughly 200,000 animals across North America, Kamchatka, and the Kurile Islands (section SM4).

We used the pairwise sequentially Markovian coalescent (PSMC) () method to estimate effective population size over time. Our results suggest that Steller’s sea cow populations have been declining for at least the last 500,000 years (Fig. 3A and section SM4), perhaps stabilizing during the warm Marine Isotope Stage 5 Interglacial, and then continued along a slow trajectory of decline until the approach loses power to detect demographic change some 10,000 to 30,000 years ago. The demographic trajectories inferred from our two individuals are concordant with that from another undated individual from the Commander Islands (Fig. 3A) (). While there is little doubt that the Steller’s sea cow’s final demise was caused either directly () or indirectly () by humans, this genetic evidence of long-term decline suggests that the species’ demographic trajectory was additionally influenced by environmental factors. We note, however, that all available samples to date come from a single location. Although Steller’s sea cows were probably capable of long-distance inter-island dispersal, undetected population structure, in addition to the small number of available high coverage genomes (N = 3), could confound this reconstructed demographic history.

Fig. 3.

Steller’s sea cow’s population history.

(A) Inferred effective population size trajectories for three Steller’s sea cow individuals. All three trajectories show a consistent decrease in population size from at least 500,000 years ago. (B) Pairwise differences between individuals are estimated as the average number of differences per 10 kbp in nonoverlapping 50-kbp regions in 10,000 randomly chosen blocks. The amount of divergence between individuals declines over time (section SM4). The vertical dashed lines correspond to the median of the number of pairwise differences per 10 kbp. Numbers on the graph represent the 14C uncalibrated date (years B.P.).

To gauge the level of genetic diversity within the Steller’s sea cow population in the nearer term past, we performed pairwise comparisons between our individuals, grouped according to radiocarbon date (Fig. 3B and section SM4). We find that Steller’s sea cow diversity was already declining some 900 to 800 years before their discovery by Steller (Fig. 3B). Our ~15×-coverage individuals, dated to approximately 1250 B.P., show ∼17 differences in 10,000 sites, a level similar to extant brown bears (). This relatively high level of genetic diversity, despite the small long-term effective population size, may have been maintained by inter-island dispersal. These results show that Steller’s sea cow populations had most probably been small for many thousands of years before their observation by Steller on Commander Island.

Our findings exemplify the power of paleogenomics to reveal unique aspects of species’ historical ecology and evolution, here addressing the phenotypic traits and population decline of Steller’s sea cow, an iconic member of the now-extinct Pleistocene megafauna. We find that Steller’s sea cows’ population decline began during the Middle Pleistocene and continued until their extinction. We also explored the genetic basis of Steller’s sea cow’s skin phenotype, exceptionally large size, and reportedly thick blubber. Here, we found evidence for convergent evolution with cetaceans in regard to inactivations of genes that are active in the extant sirenians. The thorough phenotypic record provided by Steller studied within a paleogenomics framework revealed the molecular basis of Steller’s sea cows’ peculiar skin appearance almost three centuries after their description.

MATERIALS AND METHODS

For detailed methods, please refer to the Supplementary Materials.

Sampling, DNA extraction, library preparation, and sequencing

We sampled bone fragments from 25 H. gigas individuals that were collected along the beaches of Commander Island, Russia and one sample (SNMB N51667) that is part of the collection at the Braunschweig Natural History Museum (section SM1). Of these, 12 samples contained sufficient endogenous DNA to allow sequencing of the nuclear genome (section SM1). For each specimen excluding SNMB N51667, we removed a thin layer of surface material and then powdered the remaining sample using a Mixer Mill MM400 (Restch). We divided the resulting powder and used ~30 mg for ancient DNA extraction and the remainder for radiocarbon dating, which was performed at the W. M. Keck Accelerator Mass Spectrometry Center at University of California, Irvine using ultrafiltration and following Shammas et al. () (section SM1).

With the help of a computed tomography (CT) scan (IZW Berlin; section SM1), we identified the densest regions of the petrous bone that was then drilled to obtain bone powder from SNMB N51667. This sample was radiocarbon dated at the Klaus-Tschira-AMS facility in the Curt-Engelhorn-Centre of Archaeometry in Mannheim (laboratory number GMP391).

Laboratory work for SNMB N51667 was performed in dedicated ancient DNA facilities at the University of Potsdam and included the parallel processing of blanks to detect possible contamination of laboratory equipment and reagents. DNA was extracted from 49.5 mg of bone powder following the protocol of Dabney and colleagues (). The DNA extract was converted into a single-stranded, double-indexed Illumina library and treated with uracil-DNA glycosylase and endonuclease VIII before library preparation, effectively reducing DNA damage (). The optimal cycle number for indexing was established by quantitative polymerase chain reaction [qPCR; Thermo Fisher Scientific PikoReal Real-Time PCR System ()]. The library was sequenced at the SciLifeLab facilities on an Illumina HiSeq 2500 using 151-cycle paired-end runs.

The remaining 11 samples were processed at the Paleogenomics Laboratory at the University of California, Santa Cruz in a clean laboratory facility dedicated to ancient DNA research. These samples were not treated with uracil-DNA glycosylase and endonuclease VIII before library preparation. We pretreated each powdered sample with a 0.5% sodium hypochlorite solution to remove surface contaminants (), extracted DNA as described in (), and eluted in 50 μl of 10 mM tris-HCl, 1 mM EDTA, and 0.05% Tween 20. We prepared dual-indexed, single-stranded DNA libraries (), which we pooled and sequenced at the SciLifeLab facilities on an Illumina HiSeq 2500 using 151-cycle paired-end runs adapted to double-indexed libraries. SC16.JK045 was sequenced at the University of California, San Francisco Center for Advanced Technology on an Illumina NovaSeq 6000 S4 2x100 lane.

To generate a de novo assembly of D. dugon, high–molecular weight DNA was extracted by staff at the Australian Museum from a fresh kidney sample (specimen M39374.001) using the Bioline Isolate ll Genomic DNA Kit (Bioline). DNA was shipped to the Swedish Museum of Natural History under the Convention on International Trade in Endangered Species (CITES) regulation of shipments between registered scientific institutions entitled to the exemption (Article VII, paragraph 6, of the Convention). The DNA was subsequently converted into a linked-read library using the Chromium Genome v2 kit (10x Genomics) and sequenced at the National Genomics Infrastructure (SciLifeLab) in Stockholm using two lanes on an Illumina HiSeqX instrument [2 × 150 base pairs (bp)] to reach approximately 60× genomic coverage.

Steller's sea cow sample collection and experiments are in accordance with ethics committee regulations from the University of Postdam and the University of California, Santa Cruz. The dugong tissue sample was transferred from the Australian Museum to the Swedish Museum of Natural History in accordance with the CITES Convention exemption provided by Article VII, paragraph 6.

De novo genome assembly, genome annotation, and orthology assignment

We assembled the dugong genome initially using Supernova (version 2.1.1) (), followed by another round of scaffolding using the ARKS+LINKS pipeline (versions 1.0.2 and 1.8.6) () with the 10X linked reads. Parameters used for ARKS+LINKS were as follows: -c 5, -k 30, -j 0.5, -z 500, -e 30000, -m 50-10000, -l 5, and -a 0.3.

We evaluated assembly quality using BUSCO (version 3.0.2) () with both eukaryotic and mammalian ortholog datasets. We used QUAST (version 4.5.4) () for contiguity statistics.

We annotated the dugong (D. dugon) and the manatee [T. manatus latirostris (); UCSC GCA_000243295.1] genomes using the MAKER pipeline () and RNA and protein sequences of three high-quality genomes (human, GRCh38.p13; mouse, GRCm38.p6; and elephant, Loxafr3.0) (), as previously described (). Triplet orthologs between human, mouse, and elephant were downloaded from Ensembl 98 (). Manatee and dugong genes were considered orthologs to a triplet if the ortholog assignment from MAKER agreed with the orthologous gene assigned in each triplet orthologous pair.

Data processing, selection analysis, and demographic history reconstruction

We used the Illumina Bustard software for base calling, deML for demultiplexing (), and leeHom () for adapter trimming and merging of overlapping paired-end reads. We merged sequences from a given library using SAMtools () and aligned reads from our individuals, as well as the publicly available data from Sharko et al. () to the dugong genome using Burrows-Wheeler Aligner (BWA) () with parameters adjusted for ancient DNA (). We removed PCR duplicates using bam-rmdup (https://doi.org/10.5281/zenodo.4616582) and filtered for DNA fragments of length ≥ 32 bases with a mapping quality of 30 or higher.

We used the ancient DNA damage-aware genotyper snpAD () on alignments to the dugong genome to call genotypes on scaffolds >100 kbp. We calculated pairwise distances of all samples based on mitochondrial positions with mapping quality (MQ) ≥ 30, genotype quality (GQ) ≥ 20, and coverage ≥10 identified using SAMtools mpileup () to ensure that the 12 sequenced specimens originated from different individuals.

Using the snpAD variant calling on the Steller’s sea cow SNMB N51667, which was treated with uracil-DNA glycosylase and endonuclease VIII and had the lowest ancient DNA damage (section SM1), and the dugong gene predictions, we lifted the corresponding coordinates and inferred the Steller’s sea cow sequence by replacing heterozygous/homozygous alternative positions with the alternative allele. For positions where Steller’s sea cow displayed two alternative alleles, we randomly picked one to reconstruct the sequence. To this end, we considered positions with coverage ≥5× and GQ ≥ 20, while all others were deemed as unresolved and replaced by “N”s. Codon-based alignments were obtained using MACSE version 2 (), followed by exclusion of regions containing gaps in any of the species.

We used CODEML () to scan for differently evolving genes under a branch model (model = 2 versus model = 0 compared via likelihood ratio test) using the set of orthologs as defined above in the five species, and the inferred Steller’s sea cow sequence. Significance of differently evolving genes was established using a chi-square test with one degree of freedom.

We estimated the allele frequency of stop codons in ALOXE3 and ALOX12B using the noninformative Jeffreys’ prior for a binomial distribution incorporated in the jeffreysci function from the ratesci package in R ().

For pairwise differences between individuals, we used Consensify () to generate a consensus pseudohaploid genome sequence. PSMC () was run using the parameters -t 40 and -p 4+25*2+4+6. We estimated the mutation rate using polymorphic sites between Steller’s sea cow and dugong or manatee and the divergence time for each pair of species, respectively (). Prehistorical population estimates for regions beyond Bering Island were obtained by extrapolating the number of sea cows per square kilometer of habitat at Bering Island to the square kilometer of habitat (x) in any region of interest, i.e., as (x km2 / 429 km2) × 1000.