Chromosome-level reference genome for the medically important Arabian horned viper (Cerastes gasperettii)

Abstract

Venoms have traditionally been studied from a proteomic and/or transcriptomic perspective, often overlooking the true genetic complexity underlying venom production. The recent surge in genome-based venom research (sometimes called "venomics") has proven to be instrumental in deepening our understanding of venom evolution at the molecular level, particularly through the identification and mapping of toxin-coding loci across the broader chromosomal architecture. Although venomous snakes are a model system in venom research, the number of high-quality reference genomes in the group remains limited. In this study, we present a chromosome-resolution reference genome for the Arabian horned viper Cerastes gasperettii (NCBI: txid110202), a venomous snake native to the Arabian Peninsula. Our highly contiguous genome (genome size: 1.63 Gbp; contig N50: 45.6 Mbp; BUSCO: 92.8%) allowed us to explore macrochromosomal rearrangements within the Viperidae family, as well as across squamates. We identified the main highly expressed toxin genes within the venom glands comprising the venom's core, in line with our proteomic results. We also compared microsyntenic changes in the main toxin gene clusters with those of other venomous snake species, highlighting the pivotal role of gene duplication and loss in the emergence and diversification of snake venom metalloproteinases and snake venom serine proteases for C. gasperettii. Using Illumina short-read sequencing data, we reconstructed the demographic history and genome-wide heterozigosity of the species, revealing how historical aridity likely drove population expansions. Finally, this study highlights the importance of using long-read sequencing as well as chromosome-level reference genomes to disentangle the origin and diversification of toxin gene families in venomous snake species.

Keywords: gene synteny; genomics; toxin evolution; transcriptomics; venom.

Figure 1:

(A) Reference genome for C. gasperettii, including BUSCO score, GC content, coverage level, and the main toxins found within the genome. Macrochromosomes are shown in light orange; microchromosomes are shown in bright orange. Sex chromosomes are shown in gray. DIS: disintegrins; HYAL: hyaluronidases; LAAO: L-amino acid oxidase; CRISP: cysteine-rich secreted proteins; CTL: C-type lectins. (B) HiC contact map for the macrochromosomes (above), including the sex chromosomes (Z and W), and microchromosomes (below).

Figure 2:

Chromosome-level analyses for one Elapidae (N. naja), one Crotalinae (C. adamanteus), and one Viperinae (C. gasperettii) species, with A. sagrei as the outgroup. The 4 smallest scaffolds (14, 15, 16, and 17) of A. sagrei were removed because no orthologous groups were found with other species. Borders of regions showing evidence for chromosomal rearrangements are shown in black. Estimates for branch times obtained from TimeTree.org based on divergence times between Iguania and Serpentes, Elapidae and Viperidae, and Crotalinae and Viperinae, respectively.

Figure 3:

Main toxins found in both the transcriptome and proteome of C. gasperettii. (A) Transcriptomic results with genes upregulated and exclusively found in the venom gland for both individuals. Each column represents a different tissue type per sample. Rows show the different genes, and colors correspond to different expression levels. VG: venom gland; EY: eye; OV: ovary; LI: liver; LU: lung; PA: pancreas; SP: spleen; GB: gallbladder; HE: heart; KI: kidney; LI: liver; BR: brain; TO: tongue; TE: testis. (B) Proteomic results of venom composition for a pool of 2 individuals of C. gasperettii. The pie chart displays the relative abundances of the toxin families found in the proteome of the C. gasperettii venom. PDE, phosphodiesterases.

Figure 4:

(A) Local synteny analyses for the SVMP toxin family in N. naja, C. adamanteus, and C. gasperettii. Different colors indicate orthologous genes unique to C. gasperettii, crotalids, true vipers, or elapids. ADAM28 (right) as well as flanking genes (left) are also indicated. (B) Phylogeny of SVMPs. Bold type indicates groups that contained SVMPs from C. gasperettii; purple indicates the gene is exclusively found in C. gasperettii. (C) Local synteny analyses for PLA₂ in N. naja, A. feae, C. adamanteus, and C, gasperettii. Nontoxic PLA₂ and flanking genes are also shown. (D) Phylogeny of the PLA₂ gene family, with 2 non-toxic PLA₂s as outgroups. Some samples that did not fit in any category have been removed. For a complete phylogeny see Supplementary Fig. S9. Note that PLA₂-gK is present in the phylogeny but not in the local synteny analyses, as any of the studied species contains it. (E) Local synteny analyses for SVSPs for C. adamanteus and C. gasperettii. Flanking genes are also shown. (F) Phylogeny for SVSPs with a nontoxic outgroup. For the 3 different phylogenies the groups that contained toxins from C. gasperettii are highlighted in bold.

Figure 5:

(A) Genome-wide diversity for 6 different venomous snakes: B. jararaca, C. gasperettii, C. viridis, N. kaouthia, N. naja, and S. tergeminus. (B) PSMC analysis recovering the ancient demographic history of C. gasperettii. Generation time was set to 3 years and the substitution rate to 2.4 × 10⁻⁹ per site per year. Shaded lines represent 10 bootstrap estimates. Two last glacial periods are shown with gray lines.