The genome of the zebra mussel, Dreissena polymorpha: a resource for comparative genomics, invasion genetics, and biocontrol

Abstract

The zebra mussel, Dreissena polymorpha, continues to spread from its native range in Eurasia to Europe and North America, causing billions of dollars in damage and dramatically altering invaded aquatic ecosystems. Despite these impacts, there are few genomic resources for Dreissena or related bivalves. Although the D. polymorpha genome is highly repetitive, we have used a combination of long-read sequencing and Hi-C-based scaffolding to generate a high-quality chromosome-scale genome assembly. Through comparative analysis and transcriptomics experiments, we have gained insights into processes that likely control the invasive success of zebra mussels, including shell formation, synthesis of byssal threads, and thermal tolerance. We identified multiple intact steamer-like elements, a retrotransposon that has been linked to transmissible cancer in marine clams. We also found that D. polymorpha have an unusual 67 kb mitochondrial genome containing numerous tandem repeats, making it the largest observed in Eumetazoa. Together these findings create a rich resource for invasive species research and control efforts.

Keywords: Dreissena polymorpha; RNA-Seq; genome; shell formation; stress response; thermal tolerance; zebra mussel.

Figure 1

Zebra mussel biogeography and genome sequencing strategy. (A) Photo of D. polymorpha (by Naomi Blinick). (B) Phylogenetic tree showing the evolutionary divergence between D. polymorpha and other sequenced bivalve genomes. For context, the evolutionary divergence of humans, mice, zebrafish, manta rays, nematodes, and fruit flies are shown. Grey text indicates that a genome sequence for that organism is not publicly available. Divergence times and tree construction based on Kumar et al. (2017). (C–E) Maps depicting the spread of D. polymorpha in the United States of America from 1988 through 2018. Data from US Geological Survey, Non-indigenous Aquatic Species database (USGS 2019a). (F) Map showing the extent of zebra mussel infestation in Minnesota lakes as of 2018 and depicting the location where the specimens for genome sequencing and scaffolding were collected (left). Summary of the sequencing and annotation strategy (right).

Figure 2

D. polymorpha genome and mitogenome structure and content. (A) Plots depicting the gene content, repeat and transposon density, and GC content of the 16 D. polymorpha chromosomal scaffolds. (B) Proposed circular mitochondrial genome structure. GC content plots (blue) based on 40 bp sliding window. Annotations based on sequence similarity to previously published partial mitochondrial genome (Soroka et al. 2018). Coding regions are in green and red, and the three large repeat blocks are colored turquoise, blue, and purple. (C) Plot of long (>25 kb) Oxford Nanopore (red) and PacBio (grey) reads supporting the proposed 67 kb circular mitogenome structure. Orientation of mitochondrial genome (blue) is the same as in (B).

Figure 3

SLEs in the D. polymorpha genome. (A) Schematics depicting the eight SLE copies, each with two LTRs flanking the longest ORFs among all similar elements in the D. polymorpha genome. (B) Maximum likelihood phylogenetic tree of nucleotide sequences from the RNaseH-integrase domain of Gag-Pol in D. polymorpha and other bivalve SLEs. The selected model (Anisimova et al. 2011) of DNA sequence evolution was the GTR + G (rates Γ-distributed, α = 1.190) + I (estimated proportion of invariant sites = 0.011). The tree was rooted on the Polititapes aureus 2/3/Mercenaria mercenaria branch (bottom) and bootstrap support values > 70 are shown. Colored boxes A, B, and C contain taxa involved in all HTT events within bivalves that were identified previously (Metzger et al. 2018). Arrows label HTT events 1 and 2, identified previously (Metzger et al. 2018) and HTT 3, which we identified based on the same criteria. Together these account for two independent insertions of SLEs into zebra mussels. Clade D contains SLE sequences from the zebra mussel genome; “D. polymorpha C” = chromosomal location of the SLE, with letters to order multiple insertion sites. Taxon labels include NCBI Accession number, taxon, followed by isolate number or code. ∗ = Sequence is from full length ORF encoding Gag-Pol, † = pseudogene sequence (one or more stop codons), § = sequence derived from neoplastic hemocytes (Metzger et al. 2016).

Figure 4

Tissue-specific gene expression patterns: mantle gene expression analysis. (A) D. polymorpha: lateral view of the left valve with the right valve and the covering mantle fold removed to reveal the organs dissected for transcriptomes. In purple is the margin of the mantle tissue within the left valve. In D. polymorpha, the mantle tissue is fused to form the siphons. Inhalent and exhalant siphon openings are pictured, as is the gill (ctenidium). Modified from Yonge and Campbell (2012). (B) Heatmap depicting Z-scores for tissue-specific gene expression in the foot, gill, and mantle. (C) List of the most highly expressed mantle-specific genes (tau > 0.95). (D) Gene ontology term enrichment analysis for the mantle-specific genes.

Figure 5

Byssal genes. (A) SEM image of byssus, consisting of threads and plaques. (B) Chromosomal location of the 37 loci predicted to encode 38 byssal protein variants. Chromosomal contigs (blue shaded ovals) are numbered (italics) in order of decreasing size. Byssal genes labeled above the chromosomes are on (+) strands; below are on (−) strands. Byssal protein Dpfp7 has three (α, β, γ) and Dpfp 11 has two (α, β) classes of divergent variants. Chromosome lengths and gene coordinates are in megabases (Mb). To the right of panel B are chromosomes 8 and 9, on which byssal genes are abundant. Modified from McCartney (2021).

Figure 6

Response of D. polymorpha to thermal stress. (A) Overview of experimental set-up. Animals were subjected to low (24°C), moderate (27°C), and high (30°C) thermal stress (n = 4 animals per condition). (B) Top 20 genes upregulated during moderate thermal stress by log2 fold-change. (C) Top 20 genes upregulated during high thermal stress by log2 fold-change. (D) Top 20 genes downregulated during moderate thermal stress by log2 fold-change. (E) Top 20 genes downregulated during high thermal stress by log2 fold-change. Genes highlighted in red encode chaperone proteins.