Background & Summary

Hyenas (also spelled “hyaena” in some parts of the world; Fig. 1) are among the most common large carnivores in Africa, with a widespread distribution occupying most of the habitats of the continent. There are four living species of hyena - spotted hyena (Crocuta crocuta), striped hyena (Hyaena hyaena), brown hyena (Hyaena brunnea), and aardwolf (Proteles cristata). A previous molecular systematics study suggested that hyaenids diverged from their feliform sister group 29.2 MYA, in the Middle Oligocene1. The spotted hyena is the largest member of this family and is known for its laughing call. They are fairly large in build, with body weights up to 64 kg and 55 kg for females and males, respectively2, and have relatively short torsos with lower hindquarters, and sloping backs. They have excellent night-time hearing and vision and can be found in all habitats except central Afrotropical forests, including savannas, grasslands, woodlands, forest edges, subdeserts, and even mountains up to 4,000 meters2.

Fig. 1
figure 1

A photograph of an adult C. Crocuta. The user garywalker” uploaded this image to https://pixabay.com.

The spotted hyena displays a number of unusual features that are unique among mammals. As the most numerous large predators, their prey mostly comes from ungulates, such as wildebeest, zebra, Thomson’s gazelles, cape buffalo, impala, and they also feed on insects and fishes2. Spotted hyenas have an exceptionally robust dentition, and they have the largest premolars compared with any living carnivora species of the same body size3. Adult spotted hyenas can generate powerful bite forces that are associated with their ability to capture prey with body sizes up three times larger than themselves and crush bones using their teeth4. These abilities are related to a unique caudally elongated frontal sinus in spotted hyenas that dissipate bending stresses during bone-cracking5. Unlike other carnivores, spotted hyenas are not only able to splinter the bones of large ungulates, but they are also able to digest them completely, including all organic components2.

However, perhaps the most peculiar feature of spotted hyenas is related to their reproductive biology, which in turn is directly related to their social behavior. Female spotted hyenas are about 10% larger than males and are much more aggressive, resulting in a social system where the masculinized females are dominant to all adult immigrant males6. Furthermore, females have evolved a pseudophallus as a result of a greatly elongated clitoris, the formation of which is independent of androgen hormones but may be related to estrogen signaling7. However, the behavioral aggressiveness of female hyenas and that displayed between cubs soon after parturition to establish dominance may be mediated by unusually high concentrations of androgens8. Therefore, the spotted hyena is a fascinating model species for studying the social behavior, evolution of sexual dimorphism, demography and genetic structure of a gregarious mammalian carnivore. These large predators live in societies that are far larger and more complex than those of any other mammalian carnivore and current studies of spotted hyenas are focused on the social intelligence of hyena societies9. Deciphering the genetic underpinnings of these remarkable traits would be greatly facilitated by the generation of a reference genome for spotted hyenas.

The four extant hyaenid species have a conserved karyotype of 2n = 40, with slight differences in the fundamental number of chromosomal arms. The hyaenid karyotype differs from the ancestral Carnivora karyotype by 4 fusions, 3 fissions, and at least 3 inversions as shown by comparative chromosome painting. As with the majority of autosomes, the X chromosome has a large C-positive centromeric region. G-banding patterns of the spotted hyena are very similar to those of the striped hyena10. To date, only the genomes of two striped hyenas (a female and a sample without sex information) have been sequenced and assembled11 (Genbank accession GCA_003009895.1 and GCA_004023945.1, respectively). The complete mitochondrial genomes have been generated for all four hyena species11,12. Here we present the first draft genome of a male spotted hyena (Crocuta crocuta), which will offer opportunities for unraveling the evolutionary history, population genetics and genetic underpinnings of the unique biological features of this endlessly fascinating species.

Methods

Sample collection, library construction and sequencing

Genomic DNA was obtained from a male specimen of C. crocuta (NCBI taxonomy ID: 9678; Fig. 1) stored in the Frozen Zoo® at the San Diego Zoo Institute for Conservation Research, USA (Frozen Zoo ID: KB4526).

The genomic DNA was extracted using phenol-chloroform followed by purification using ethanol precipitation13. The extracted DNA was run and visualized on a 1.5% agarose gel run in 1x TBE buffer to check for the presence of high molecular weight DNA. DNA concentration and purity were quantified on a NanoDrop 2000 spectrophotometer and Qubit 2.0 Fluorometer (Thermo Fisher Scientific, USA) before shipping to BGI-Shenzhen, China. We obtained a total of 372 µg of genomic DNA, with a concentration of 0.418 µg/µL using the Nanodrop 2000 and 0.245–0.399 µg/µL based on four replicate readings using the Qubit 2.0 Fluorometer. The 260/280 ratio of purity was 1.95. We then barcoded the sample using cytochrome b (Cytb) gene. Then, according to the gradient library strategy, we constructed 13 insert-size libraries, with the following insert size lengths: 170 bp, 500 bp, 800 bp, 2 kbp, 5 kbp, 10 kbp, 20 kbp. We used the HiSeq. 2000 sequencer (Illumina, USA) to sequence Paired-End (PE) reads for each library across 14 lanes. A total of about 299 Gb raw data was generated from 13 libraries, achieving a sequencing depth (coverage) of 149.25 (Table 1).

Table 1 Statistics of raw read data, assuming the genome size is 2.0 Gb.

Quality control

To minimize misassembly errors, we filtered raw reads prior to de novo genome assembly according to the following two criteria. First, reads with more than 10 bp aligned to the adapter sequence (allowing <= 3 bp mismatch) were removed. Second, reads with 40% of bases having a quality value less than or equal to 10 were discarded. Finally, we obtained 190.4 G data with a coverage of 95.2 (Table 2).

Table 2 Data statistics following filtering of raw read data.

Estimation of genome size

Three short-insert libraries (two of 170 bp and one of 500 bp) were used to estimate the genome size and genome-wide heterozygosity by k-mer analysis. A total of about 385 M PE reads were submitted to jellyfish14 to calculate k-mer frequency. Then the k-mer distribution was illustrated by Genomescope7 with parameters “k = 17; length = 100; max coverage = 1000”. We obtained an estimated genome size 2,003,681,234 bp, and heterozygosity of 0.325% (Fig. 2).

Fig. 2
figure 2

17-mer estimate of genome size. The x-axis is depth (X), the y-axis is the proportion which represents the frequency at that depth divided by the total frequency of all coverage depths. Without consideration of the sequence error rate, heterozygosity rate, and repeat rate of the genome, the 17-mer distribution should approximate a Poisson distribution.

Genome assembly and assessment

SOAPdenovo (V1.06)15 was employed to assemble the genome de novo, following the filtering of the short insert size data and removing the small peak of the large insert size data. The SOAPdenovo assembly algorithm included three main steps. (1) Contig construction: the short-insert size library data were split into k-mers and constructed using a de Bruijn graph, which was simplified by removing tips, merging bubbles, removing the low coverage of the connection and removing small repeats. We obtained the contig sequence by connecting the k-mer path, resulting in a contig N50 2,104 bp, and total length 2,295,545,898 bp. (2) Scaffold construction: we obtained 80% of all aligned paired-end reads by realigning all usable read on contigs. Then we calculated the amount of shared paired-end relationships between each pair of contigs, weighted the rate of consistent and conflicting paired-ends, and then constructed the scaffolds step by step. As a result, we obtained scaffolds with an N50 7,168,038 bp, and total length 2,355,303,269 bp from short insert-sized paired-ends, to long distant paired-ends. (3) Gap closing: To fill the gaps inside the constructed scaffolds, we used the paired-end information to retrieve the read pairs to do a local assembly again for these collected reads. In summary, we closed 87.7% of the intra-scaffold gaps, or 85.8% of the sum gap length. The contig N50 size increased from 2,104 bp to 21,301 bp (Table 3). The scaffold assembly size was 2,355,303,269 bp, which is close to the assembly-based genome size of 2,374,716,107 bp reported for the striped hyaena, Hyaena hyaena11 (NCBI accession: ASM300989v1). We also retrieved and annotated the mitochondrial genome of the spotted hyena using the MitoZ program16, which has a length of 16,858 bp, similar to the first mitochondrial genomes sequenced for this species12.

Table 3 Statistics of the assembled sequence length.

Assessment of the draft genome was performed by looking at the completeness of single-copy orthologs using BUSCO (version 3.1.0)17, searching against Mammaliaodb9 database which contains 4,104 single-copy ortholog groups. A total of 95.5% of the orthologs were identified as complete, 2.5% as fragmented and 2.0% as missing, indicating an overall high quality of the spotted hyena genome assembly. Given that 99.95% of the short scaffolds (<1k) harbored only 1.2% of the total genome length, we excluded these scaffolds for downstream analysis, including repetitive element and gene feature annotation.

Repetitive element annotation

Both tandem repeats and transposable elements (TE) were searched for and identified across the C. crocuta genome. Tandem repeats were identified using Tandem Repeats Finder (TRF, v4.07)18 and transposable elements (TEs) were identified by a combination of homology-based and de novo approaches. For the homology-based prediction, we used RepeatMasker version 4.0.619 with the settings “-nolow -no_is -norna -engine ncbi” and RepeatProteinMask (a program within RepeatMasker package) with the settings “-engine ncbi -noLowSimple -pvalue 0.0001” to search TEs at the nucleotide and amino acid level based on known repeats (Fig. 3). RepeatMasker was applied for DNA-level identification using a custom library which combined the Repbase21.10 dataset20. At the protein level, RepeatProteinMask was used to perform RMBlast against the TE protein database. For ab initio prediction, RepeatModeler (v1.0.8)21 and LTR_FINDING (v1.06)22 were applied to construct the de novo repeat library. Contamination and multi-copy sequences in the library were removed and the remaining sequences were classified according to the BLAST result following alignment to the SwissProt database. Based on this library, we used RepeatMasker to mask the homologous TEs and classified them (Fig. 4). Overall, a total of 826 Mb of repetitive elements were identified in the spotted hyena, comprising 35.29% of the whole genome (Table 4).

Fig. 3
figure 3

Distribution of divergence rate of each type of transposable element (TE) in the Crocuta crocuta genome assembly based on homology-based prediction. The divergence rate was calculated between the identified TEs in the genome using a homology-based method and the consensus sequence in the Repbase database20.

Fig. 4
figure 4

Distribution of divergence rate of each type of TE in the Crocuta crocuta genome assembly based on ab initio prediction. The divergence rate was calculated between the identified TEs in the genome by ab initio prediction and the consensus sequence in the predicted TE library (see Methods).

Table 4 Transposable element content of the Crotuta crotuta genome assembly.

Protein-coding gene annotation

We used ab initio prediction and homolog-based approaches to annotate protein-coding genes as well splicing sites and alternative splicing isoforms. Ab initio prediction was performed on the repeat-masked genome using gene models from human, domestic dog, and domestic cat using AUGUSTUS (version 2.5.5)23, GENSCAN24, GlimmerHMM (version 3.0.4)25, and SNAP (version 2006-07-28)26, respectively. A total of 22,789 genes were identified by this method. Homologous proteins of, Homo sapiens, Felis catus and Canis familiaris (from the Ensembl 96 release) were mapped to the spotted hyena genome using tblastn (Blastall 2.2.26)27 with parameters “-e 1e-5”. The aligned sequences as well as their query proteins were then submitted to GeneWise (version 2.4.1)28 for searching an accurate spliced alignment. The final gene set (22,747) was collected by merging ab initio and homolog-based results using a customized pipeline (Table 5).

Table 5 General statistics of the number of protein-coding genes based on ab initio (de novo) and homology-based prediction methods.

Gene function annotation

Gene functions were assigned according to the best match obtained by aligning translated gene coding sequences using BLASTP with parameters “-e 1e-5” to the SwissProt and TrEMBL databases (Uniprot release 2017-09). The motifs and domains of genes were determined by InterProScan (v5)29 against protein databases including ProDom30, PRINTS31, Pfam32, SMART33, PANTHER34 and PROSITE35. Gene Ontology IDs for each gene were obtained from the corresponding SwissProt and TrEMBL entries. All genes were aligned against KEGG proteins, and the pathway in which the gene might be involved was derived from the matched genes in the KEGG database36. In summary, 22,166 (97.45%) of the predicted protein-coding genes were successfully annotated by at least one of the six databases (Table 6).

Table 6 Number of genes with predicted homology or functional classification according to alignment to different protein databases.

Gene family construction and phylogeny reconstruction

To gain insight into the phylogenetic history and evolution of gene families of Crocuta crocuta, we clustered gene sequences of seven species (Felis catus, Canis familiaris, Ailuropoda melanoleuca, Crocuta crocuta, Panthera pardus, Panthera leo, Panthera tigris altaica) and Homo sapiens as the outgroup (Ensembl release-96, Panthera leo from unpublished data) into gene families using orthoMCL (v2.0.9)37. The protein-coding genes for the eight species were retrieved by selecting the longest transcript isoform for each gene for downstream pairwise assignment (graph building). We performed an all-against-all BLASTP search on the protein sequences of all the reference species, with an E-value cut-off of 1e-5. Gene family construction employed the MCL algorithm38 with the inflation parameter of ‘1.5’. A total of 16,271 gene families of C. crocuta, H. sapiens, F. catus, A. melanoleuca were clustered. There were 11,671 gene families shared by these four species, while 292 gene families containing 1,446 genes were specific to C. Crocuta (Fig. 5). Noticeably, the gene families C. crocuta and F. catus shared were less than C. crocuta and H. sapiens shared, which could result from that H. sapiens had a more complete genome and annotation.

Fig. 5
figure 5

Venn diagram showing comparison of shared and unique protein-coding genes among spotted hyena, human, domestic cat and domestic dog based on orthology analysis.

We identified 6,601 single-copy orthologous genes to reconstruct the phylogenetic tree of the eight species. Multiple sequence alignments of amino acid sequences for each gene were generated using MUSCLE (version 3.8.31)39, and trimmed using Gblocks (0.91b)40, achieving well-aligned regions with the parameters “-t = p -b3 = 8 -b4 = 10 -b5 = n -e = -st”. We performed phylogenetic analysis using the maximum-likelihood method as implemented in PhyML (v3.0)41, using the JTT + G + I model for amino acid substitution (Fig. 6). The root of the tree was determined by minimizing the height of the whole tree via Treebest (v1.9.2; http://treesoft.sourceforge.net/treebest.shtml). Finally, we estimated the divergence time among the eight lineages using MCMCTree from the PAML version 4.4 software package42. Two priors based on the fossil record were used to calibrate the substitution rate, including Boreoeutheria (91-102 MYA) and Carnivora (52-57 MYA)43. Consistent with previous studies, the spotted hyena groups with the four species included from the Felidae in a clade defining the suborder Feliformia, which diverged from the Caniformia (represented by the domestic dog and giant panda) 53.9 Mya44.

Fig. 6
figure 6

Phylogenetic tree of C. crocuta and seven other species constructed by the maximum likelihood method based on 6,601 single-copy orthologues. The divergence time was estimated using the two calibration priors derived from the Time Tree database (http://www.timetree.org), which are marked by a red rhombus. All estimated divergence times are shown with 95% confidence intervals in brackets.

Data Records

Raw reads from Illumina sequencing are deposited in the NCBI Sequence Read Archive (SRA) database with accession numbers SRP21580045, and Bioproject accession PRJNA554753 and are also deposited in CNGB Nucleotide Sequence Archive (CNSA) database with accession number CNR0105011-CNR010502346,47,48,49,50,51,52,53,54,55,56,57,58 and Bioproject accession CNP0000511. The genome assembly of C. crocuta generated in this study was deposited in NCBI Assembly under the accession number GCA_008692635.159 and in CNSA with the accession number CNA000352060. Copies of all of annotation outputs including genes, functional assignments, and copies of the gene families and its statistics, and final tree in newick format for 8 species are deposited in figshare database61.

Technical Validation

DNA quality control

Quantification of the DNA sample using both NanoDrop and a DNA fluorometer were performed before library construction (see method). DNA sample was also identified by DNA barcode of Cytb gene to avoid a mislabeling.

Comparison for genome assembly of striped hyena

The previous released genome assemblies of striped hyena (Hyaena hyaena) were relatively fragmented and had scaffold N50 of 66,490 bp and 2,001,327 bp, respectively, which are significantly shorter than the presented spotted hyena genome assembly has scaffold N50 of 7,168,038 bp (Table 7). We chose the better assembled genome of striped hyena (GCA_003009895.1) to conduct a synteny analysis with spotted hyena. The synteny analysis revealed that the two genome assemblies had overall a high ratio of syntenic region, with 96% for spotted hyena and 90.5% for striped hyena can be mapped to each other. However, the spotted hyena genome assembly has less breakpoints and has more contiguous than striped hyena genome (Fig. 7). On average, each spotted hyena scaffold can be mapped to 1.07 scaffolds of striped hyena. Overall, the spotted genome assembly has much higher quality and can be valuable for future comparative genomics study for hyena and other mammals.

Table 7 Overall statistic of syntenic analysis between spotted hyena and striped hyena.
Fig. 7
figure 7

A case of syntenic relationship of spotted hyena genome with striped hyena genome assembly (GCA_003009895.1). The right dark green scaffold belongs to spotted hyena, whereas the left light green scaffolds belong to striped hyena. Links connect the location of homeologs blocks between two genome assemblies, based on the comparison of sequence information using LASTZ (version 1.02.00, http://www.bx.psu.edu/~rsharris/lastz/) under the default settings. The scaffolds of length >1 Mb and links (syntenic block) of length >100 kb are showed on figure.