# Environmental transcriptomes of invasive dreissena, a model species in ecotoxicology and invasion biology

Romain Péden, Pascal Poupin, Bénédicte Sohm, Justine Flayac, Laure Giambérini, Christophe C. Klopp, Fanny Louis, Sandrine Pain-Devin, Marine Potet, Rémy-Félix Serre, et al.

► To cite this version:

Romain Péden, Pascal Poupin, Bénédicte Sohm, Justine Flayac, Laure Giambérini, et al. Environ- mental transcriptomes of invasive dreissena, a model species in ecotoxicology and invasion biology. Scientific Data, 2019, 6 (1), pp.234. 10.1038/s41597-019-0252-x. hal-02337309

HAL Id: hal-02337309

https://hal.science/hal-02337309

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN

DATA DESCRIPTOR

# Environmental transcriptomes of invasive dreissena, a model species in ecotoxicology and invasion biology

Romain Péden<fnref n="1" />,²,⁵, Pascal Poupin<fnref n="1" />,⁵, Bénédicte Sohm<fnref n="1" />,⁵, Justine Flayac<fnref n="1" />, Laure Giambérini<fnref n="1" />, Christophe Klopp<fnref n="1" />,³, Fanny Louis<fnref n="1" />, Sandrine Pain-Devin<fnref n="1" />, Marine Potet<fnref n="1" />, Rémy-Félix Serre⁴ & Simon Devin<fnref n="1" />*

Dreissenids are established model species for ecological and ecotoxicological studies, since they are sessile and filter feeder organisms and reflect *in situ* freshwater quality. Despite this strong interest for hydrosystem biomonitoring, omics data are still scarce. In the present study, we achieved full *de novo* assembly transcriptomes of digestive glands to gain insight into *Dreissena polymorpha* and *D. rostriformis bugensis* molecular knowledge. Transcriptomes were obtained by Illumina RNA sequencing of seventy-nine organisms issued from fifteen populations inhabiting sites that exhibits multiple freshwater contamination levels and different hydrosystem topographies (open or closed systems). Based on a recent *de novo* assembly algorithm, we carried out a complete, quality-checked and annotated transcriptomes. The power of the present study lies in the completeness of transcriptomes gathering multipopulational organisms sequencing and its full availability through an open access interface that gives a friendly and ready-to-use access to data. The use of such data for proteogenomic and targeted biological pathway investigations purpose is promising as they are first full transcriptomes for this two *Dreissena* species.

## Background & Summary

*Dreissena polymorpha* and *D. rostriformis bugensis*, also known as the zebra and the quagga mussels, are model species in ecology and ecotoxicology since the early eighties. These invasive species are now commonly found in freshwaters of the northern hemisphere. *D. polymorpha* belongs to the 100 of the World's Worst Invasive Alien Species identified by the IUCN (http://www.iucngisd.org/gisd/100_worst.php). It spread from Ponto-Caspian basin to northern and western Europe in the 19th century, and to North America and Middle East (Turkey) in the late 20th century, while *D. r. bugensis* colonized both North America and Western Europe between the late 20th and at the beginning of the 21st century<fnref n="1" />,². They can reach high densities in invaded ecosystems³, inducing important ecological and economic damages⁴,⁵. Their distribution and invasion dynamic is well documented, as well as ecological features of colonized ecosystems. Dreissenids thus become model species in ecology of biological invasion to identify pathway of invasion, and genetic mechanisms associated to colonization at various spatial and temporal scales.

Their ability to tolerate a wide range of environmental contaminants and their presence in almost every major hydrosystem in Europe and North America also lead dreissenids to become a model species for ecotoxicologists, who identified them as the counterpart of *Mytilus* for freshwaters⁶. Indeed, being abundant, sessile, filter feeder, mussels bioaccumulate contaminants present in the water column, and are good candidates to be included in biomonitoring programs. It allows to evaluate both the presence of contaminants through bioaccumulation measurements⁷,⁸ and their effects through the use of biomarkers, either in the field or in laboratory⁹⁻¹¹.

<footnote n="1">¹Université de Lorraine, CNRS, LIEC, F-57000, Metz, France. ²Université Reims Champagne Ardenne, UMR-I 02 SEBIO, S1097, Reims, France. ³Plate-forme bio-informatique Genotoul, Mathématiques et Informatique Appliquées de Toulouse, INRA, 31326, Castanet-Tolosan, France. ⁴INRA, US 1426, GeT-PlaGe, Genotoul, INRA Auzeville, Castanet Tolosan, Cedex, France. ⁵These authors contributed equally: Romain Péden, Pascal Poupin and Bénédicte Sohm. *email: simon.devin@univ-lorraine.fr</footnote>
A bibliographic research with the keyword “dreissena” performed on the Web of Science the 22th of january 2019 reveals that 3,667 articles were published since 1990, with *ca.* 180 article each year since 2010. In the meantime, these articles were cited 97,000 times, with a h-index of 119. Their distribution among scientific topics (Fig. 1) reveals this dual interest in ecology and ecotoxicology. However, over the eleven Bioprojects present in NCBI for *Dreissena*, most focused on foot, gill or mantle tissues where few concern digestive glands. Moreover, the majority of their studies are not *de novo* RNA-seq but micro-array studies and almost none of them made sequences available in accessible repositories. In the present study, digestive gland was chosen for its importance in ecotoxicology due to its central roles in detoxication functions and energy metabolism¹².

Indeed, genetic studies on zebra and quagga mussel mainly concerns their spread dynamic through a population genetic perspective¹³⁻¹⁵, the identification of genetic markers for their identification, and bivalve phylogeny and phylogeography¹³,¹⁶⁻¹⁹. When the focus of the study was to elucidate toxic effect pathways, only a very limited set of genes were targeted, corresponding to the main processes implied in organism's response to contaminant²⁰,²¹.

Studies from several populations of the two species can pave the way for several new research possibilities. The most promising and useful ones are to make available a reference transcriptome for proteogenomic studies, to identify sequences to develop new sets of biomarkers, and to better understand acclimation mechanisms occuring during the invasion process and implied in tolerance to contamination.

## Methods

**Sampling and RNA extraction.** Samples were performed on 12 sites representative of french dreissenid populations diversity. Sites were chosen according to their variable contamination levels, chemicophysical quality and topography characteristics (see Supplementary Table 1). Five bivalves were sampled by populations and digestive glands were *in situ* dissected and stored in RNAlater (Qiagen) at 4°C. Total RNA extraction was done using RNeasy MinElute kit (Qiagen) according to the manufacturer's protocol with slight modifications. Briefly, tissue disruption was done first with a pestle and then with a beads beater with glass bead (200 µm) for 1 min. at max speed in 1 mL of RLT buffer supplemented with 40 mM DTT. Samples were centrifuged 3 min at 20,000 × g. Qiagen protocol was used for the next steps. Genomic DNA was removed by DNase digestion on column and total RNA was eluted in ultra pure water. RNA purity and quantity were assessed by OD measurements (OD 260 nm and OD ratio 260/280 and 260/230) and RNA integrity was checked using Bioanalyseur 2100 (Agilent, CA, USA). Total RNA samples were stored at -80. Samples were send to GeT-PlaGe core facility in dried-ice.

**RNA sequencing.** RNAseq was performed at the GeT-PlaGe core facility, INRA Toulouse. RNA-seq libraries have been prepared according to Illumina's protocols using the Illumina TruSeq Stranded mRNA sample prep kit to analyze mRNA. The 79 individuals were sequenced separately on 8 different lanes of one flow cell. Briefly, mRNA were selected using poly-T beads. Then, RNA were fragmented to generate double stranded cDNA and adaptors were ligated to be sequenced. Eleven cycles of PCR were applied to amplify libraries. Library quality was assessed using a Fragment Analyser and libraries were quantified by QPCR using the Kapa Library Quantification Kit. RNA-seq experiments have been performed on an Illumina HiSeq3000 using a paired-end read length of 2 × 150 pb with the Illumina HiSeq3000 sequencing kits.

**De novo assembly and annotation.** The RNA-seq libraries read quality was evaluated using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Reads from 45 *D. polymorpha* and 34 *D. rostriformis bugensis* individuals were cleaned, filtered and *de novo* assembled for each species with DRAP (*De novo* RNA-seq Assembly Pipeline, version 1.7)²² using the Oases assembler²³. Contigs were kept if they had at least one FPKM. Contigs were then aligned with NCBI BLAST (version 2.2.26, e-value under 1e-5 parameter) on Refseq, Swissprot and three databases issued from Ensembl (*Crassostrea gigas*, *Lottia gigantea* and *Lingula anatina* sequences) to retrieve corresponding annotations. Contigs were also processed with RNAmmer (version 1.2, standard parameters)²⁴ to find ribosomal genes, with RepeatMasker (version open-4-0-3, -engine crossmatch -gccalc -species *Crassostrea gigas* parameters, http://www.repeatmasker.org) to list contained repeats and with InterProScan (version 4.8, -goterms -pathways parameters)²⁵ for gene ontology and structural annotation. Reads were realigned back to contigs with BWA (version 0.7.12, standard parameters, mem algorithm)²⁶. The resulting sam files were compressed, sorted and indexed with SAMtools (view, sort and index programs, version 1.1, standard parameters)²⁷. Contig expression counts were generated from the bam files with SAMtools (IdxStats program, version 1.1, standard parameters) and merged with UNIX commands (cut, paste). Alignment files were then filtered for duplicates with SAMtools (rmdup program, version 1.1, standard parameters) before variant calling (SNPs and INDELs). Resulting bam files were processed with GATK (version v3.0-0-g6bad1c6, -glm BOTH parameter) following the best practices found on the GATK website²⁸. All the results were uploaded in a RNAbrowse instance²⁹ and can be accessed from the web at http://ngspipelines.toulouse.inra.fr:9014/. Further KEGG annotation were done with KAAS interface (KEGG Automatic Annotation Server) using contig sequences in fasta format.

Biomarker explorations were done using the implemented BLAST in NGSpipeline interface. Biomarker sequences come from close related species (i.e. molluscs). Candidates were selected by the user on score, e-values, identity and length. If several candidates exist, the sequence with highest score was selected.

## Data Records

Raw reads and assemblies were gathered in the same NCBI BioProject (PRJNA507340) which includes all BioSamples used for transcriptome assembly (Table 1)³⁰⁻³². All datasets were also available online on sequencing platform web interface (Table 2)³³. Datasets annotations (contigs, KEGG, GO) are available on Figshare (Table 3)³⁴.
<figure><img src="image_2.png" /><figcaption>Fig. 1 Top 10 fields treemap for *Dreissena* publications.</figcaption></figure>

<table><caption>Table 1. BioProject deposit. The BioProject gathered all BioSamples, SRA³⁰ and TSA³¹,³² related to this Data Descriptor.</caption><thead><tr><th>BioProject</th><th>Datasets (nb)</th><th>Accessions</th></tr></thead><tbody><tr><td rowspan="3">PRJNA507340</td><td>BioSamples (79)</td><td>SAMN10537936 to SAMN10538014</td></tr><tr><td>SRA (79)</td><td>SRR8354718 to SRR8354796</td></tr><tr><td>TSA (2)</td><td>GHIW00000000 and GHIX00000000</td></tr></tbody></table>

<table><caption>Table 2. NGSPipeline deposit. Datasets are fully accessible in a user friendly web interface provided by Genotoul sequencing and bioinformatic platform²⁹,³³. Biomart and BLAST are implemented for quick database interrogation.</caption><thead><tr><th>Data type</th><th>URL</th><th>Implemented softwares</th></tr></thead><tbody><tr><td>Web interface</td><td>http://ngspipelines.toulouse.inra.fr:9014</td><td>BioMart and BLAST</td></tr></tbody></table>

<table><caption>Table 3. Figshare deposit. Supplementary files are available on figshare including annotations (GO, KEGG) and ready to use transcriptomes in fasta format³⁴.</caption><thead><tr><th>Descriptive filename</th><th>Data format</th></tr></thead><tbody><tr><td>Dreissena_polymorpha_contig_sequences</td><td> fasta </td></tr><tr><td>Dreissena_rostriformis_bugensis_contig_sequences</td><td> fasta </td></tr><tr><td>BLAST_annotations_for_Dreissena_polymorpha_contigs</td><td> csv </td></tr><tr><td>BLAST_annotations_for_Dreissena_rostriformis_bugensis_contigs</td><td> csv </td></tr><tr><td>GOterms_annotations_for_Dreissena_polymorpha_contigs</td><td> csv </td></tr><tr><td>GOterms_annotations_for_Dreissena_rostriformis_bugensis_contigs</td><td> csv </td></tr><tr><td>KEGG_annotations_for_Dreissena_polymorpha_contigs</td><td> csv </td></tr><tr><td>KEGG_annotations_for_Dreissena_rostriformis_bugensis_contigs</td><td> csv </td></tr></tbody></table>

## Technical Validation

**Extraction and RNA integrity.** Total RNA purity was assessed with a Nanodrop ND-100 Spectrophotometer (Nanodrop Technologies, Wilmington, USA) and RNA with a 260/280 and 260/230 ratio superior to 1.8 were kept. RNA integrity was evaluated with a Bioanalyzer (Agilent RNA 6000 Nano kit). Due to a non conventional 18S/28S ribosomal ratio in bivalve, sample quality was evaluated with the 18S/28S ratio and on the electropherogram (Supplementary Fig. S1). The absence of degradation fragment in the 5S, Fast and Inter regions was a criteria of selection for sequencing.

**De novo transcriptome assembly validation.** A total of 94,217 contigs were assembled spanning from 200 to 40,000 bp and with an average length of 2,314 bp for *D. polymorpha* and 1,972 for *D. r. bugensis* (Table 4). Assembly validation aims at verifying the correspondence between contigs and assembled reads, between contigs
<table><caption>**Table 4.** Assembly metrics.</caption><thead><tr><th>Type</th><th>D. polymorpha</th><th>D. r. bugensis</th></tr></thead><tbody><tr><td>Number of contigs</td><td>44,538</td><td>49,679</td></tr><tr><td>Total size of contigs (bp)</td><td>103,039,811</td><td>97,982,186</td></tr><tr><td>N50 (bp)</td><td>3,094</td><td>2,674</td></tr><tr><td>Average length (bp)</td><td>2,314</td><td>1,972</td></tr><tr><td>Longest contig (bp)</td><td>39,481</td><td>39,593</td></tr><tr><td>Shortest contig (bp)</td><td>210</td><td>207</td></tr></tbody></table>

<table><caption>**Table 5.** Individual realignment statistics.</caption><thead><tr><th>Type</th><th>D. polymorpha</th><th>D. r. bugensis</th></tr></thead><tbody><tr><td>Total read count average</td><td>69,815,940 ±9,645,652</td><td>73,730,343 ±9,641,714</td></tr><tr><td>Alignment rate</td><td>97.56% ±0.7%</td><td>98.15% ±0.05%</td></tr><tr><td>Properly paired rate</td><td>88.84% ±1.15%</td><td>86.58% ±1.60%</td></tr></tbody></table>

<table><caption>**Table 6.** Crassostrea gigas proteins.</caption><thead><tr><th>Type</th><th>D. polymorpha</th><th>D. r. bugensis</th></tr></thead><tbody><tr><td>Number of C. gigas proteins matching</td><td>4,352</td><td>4,281</td></tr><tr><td>Number of unique Dreissenid contigs mapped to</td><td>3,594</td><td>3,558</td></tr></tbody></table>

<table><caption>**Table 7.** BUSCO analysis.</caption><thead><tr><th>Type</th><th>D. polymorpha</th><th>D. r. bugensis</th></tr></thead><tbody><tr><td>Complete BUSCOs</td><td>921 94.2%</td><td>923 94.4%</td></tr><tr><td>Complete and single-copy BUSCOs</td><td>494 50.5%</td><td>450 46.0%</td></tr><tr><td>Complete and duplicated BUSCOs</td><td>427 43.7%</td><td>473 48.4%</td></tr><tr><td>Fragmented BUSCOs</td><td>5 0.5%</td><td>7 0.7%</td></tr><tr><td>Missing BUSCOs</td><td>52 5.3%</td><td>48 4.9%</td></tr><tr><td>Total BUSCO groups searched</td><td>978 100%</td><td>978 100%</td></tr></tbody></table>

and the proteome of a phylogenetically related species and to check if the contigs host single copy awaited proteins. First, reads were mapped back to contigs for each sample in order to monitor the realignment rates revealing a realignment rate above 95% for both species (Table 5). Second, *Crassostrea gigas* proteins were aligned to contigs using BLAT (v. 35 × 1)³⁵. Proteins mapping at 50% identity and over 50% of the protein length were counted to measure the share of proteins correctly reconstructed by the assembly (Table 6). Last, contigs were processed with BUSCO version 3.0.2³⁶ using the metazoa OrthoDB (v. 9) database to check for awaited proteins (Table 7) which shows that almost 95% of expected metazoa BUSCO sequences were found as complete in our *de novo* transcriptome.

**Annotation quality.** More than the half of assembled contigs found annotation with an e-value under 1e-5 (59% of 44,538 *D. polymorpha* contigs and 57% of 49,679 *D. rostriformis bugensis* contigs). Among top 5 species found as best annotation, *Crassostrea gigas* represent the best hit species matching with approximately 53% of *D. polymorpha* and *D. r. bugensis* sequences (Fig. 2). Among the “others” category, some non bivalves sequences can be present thanks to the analyzed tissue (as digestive gland may contain processing foods). Sequence functional annotations successfully attribute GOterms to 35,4% and KEGG to 21.6% of *D. polymorpha* contigs and attribute GOterms to 34,2% and KEGG to 20.6% of *D. r. bugensis* contigs. Full contigs annotation as well as functional annotations are available in the figshare deposit in csv format for quick and easy reuse (Table 3).

**Usage Notes**

We present here the first whole *de novo* transcriptomes of the two invasive *Dreissena polymorpha* and *D. rostriformis bugensis* species who will be fully available in TSA repository. One of the strength of proposed transcriptome is the conditions in which it was obtained: by investigating individuals coming from several populations, inhabiting contrasted hydrosystems, we encompassed as much as possible the variability of the transcriptome induced by environmental conditions (pollution, biotic interactions, physico-chemistry, climate...), with a mRNA set that is representative of this ecological heterogenity. Our ambitious approach also lead to identification
<figure><img src="image_3.png" /><figcaption>Fig. 2 Most representated species hits. Top 5 best species hits for (A) *Dreissena polymorpha* and (B) *Dreissena rostriformis bugensis*.</figcaption></figure>

<table><caption>**Table 8.** Biomarker best contig hits and closest species hits. *D. polymorpha* and *D. rostriformis bugensis* databases were used to find common ecotoxicological biomarkers. Contig coding sequences were then screened in NCBIr (*blastn*) and closest species (excluding *D. polymorpha* and *D. rostriformis bugensis*) were reported according with E-values and identity percentages. GABARAP: Gamma-aminobutyric acid receptor-associated protein; mTOR: mechanistic target of rapamycin; MRP1: Multidrug resistance-associated protein 1; MDR1: Multidrug resistance protein 1; Cds: coding sequence.</caption><thead><tr><th rowspan="3">Biomarkers</th><th colspan="2">Contig database name</th><th rowspan="3">Species</th><th colspan="4">Closest blastn hit</th></tr><tr><th rowspan="2">D. polymorpha</th><th rowspan="2">D. r. bugensis</th><th colspan="2">E-value</th><th colspan="2">% Identity</th></tr><tr><th>D. p.</th><th>D. r. b.</th><th>D. p.</th><th>D. r. b.</th></tr></thead><tbody><tr><td>Acetylcholinesterase</td><td>Dpatter105324424.2.2</td><td>DbLoc105324424.3.3</td><td>Crassostrea virginica</td><td>2e-8</td><td>4e-12</td><td>68%</td><td>67%</td></tr><tr><td>Beta Actine</td><td>DpAct1.3</td><td>DbAct2.4</td><td>Crassostrea virginica</td><td>0.0</td><td>0.0</td><td>90%</td><td>90%</td></tr><tr><td>Catalase</td><td>DpLoc105339902</td><td>DbLoc105339902.2.4</td><td>Corbicula fluminea</td><td>0.0</td><td>0.0</td><td>77%</td><td>76%</td></tr><tr><td>Superoxide dismutase (Cu-Zn)</td><td>DpSODC</td><td>DbSODC</td><td>Ostrea edulis</td><td>4e-55</td><td>1e-48</td><td>73%</td><td>72%</td></tr><tr><td>Superoxide dismutase (Mn)</td><td>DpLoc101852344</td><td>DbLoc101852344</td><td>Haliotis fulgens</td><td>3e-78</td><td>3e-77</td><td>72%</td><td>72%</td></tr><tr><td>Estrogen receptor</td><td>DpLoc105318922.2.2</td><td>DbLoc105318922.1.2</td><td>Ruditapes philippinarum</td><td>0.0</td><td>0.0</td><td>75%</td><td>74%</td></tr><tr><td>GABARAP</td><td>DpLoc105335545.1.2</td><td>DbLoc105335545.1.6</td><td>Meretrix meretrix</td><td>2e-103</td><td>1e-98</td><td>84%</td><td>83%</td></tr><tr><td>Glutathione peroxidase 1</td><td>DpGPX1.1.2</td><td>DbGPX1.1.2</td><td>Ruditapes philippinarum</td><td>1e-54</td><td>6e-59</td><td>72%</td><td>73%</td></tr><tr><td>Glutathione peroxidase (Se)</td><td>DpLoc106070504.1.5</td><td>DbLoc106070504.6.8</td><td>Meretrix meretrix</td><td>3e-79</td><td>6e-75</td><td>73%</td><td>73%</td></tr><tr><td>Heat Shock Protein70</td><td>DpHSP7D.3.8</td><td>DbHSP7D.4.8</td><td>Corbicula fluminea</td><td>0.0</td><td>0.0</td><td>80%</td><td>81%</td></tr><tr><td>Metallothionein (isoform 1)</td><td>DpMT.6.6</td><td>DbMT.6.6</td><td>Meretrix meretrix</td><td>6e-05</td><td>6e-05</td><td>80%</td><td>80%</td></tr><tr><td>mTOR</td><td>DpLoc105331599</td><td>DbLoc105331599</td><td>Crassostrea gigas</td><td>0.0</td><td>0.0</td><td>72%</td><td>72%</td></tr><tr><td>Na/K ATPase</td><td>DpLoc106058320</td><td>DbLoc106058320.1.2</td><td>Tridacna squamosa</td><td>0.0</td><td>0.0</td><td>76%</td><td>77%</td></tr><tr><td>Succinate dehydrogenase</td><td>DpLoc101864456.1.3</td><td>DbLoc105338659</td><td>Aplysia californica</td><td>9e-14</td><td>2e-16</td><td>66%</td><td>67%</td></tr><tr><td>Acid phosphatase</td><td>DpLotGIDRAFT_139839</td><td>DbContig19914</td><td>Pomacea canaliculata</td><td>1e-11</td><td>0.0</td><td>67%</td><td>95%</td></tr><tr><td>Thioredoxin reductase 1</td><td>DpLoc105322705.2.2</td><td>DbLoc105322705</td><td>Crassostrea spp.</td><td>0.0</td><td>0.0</td><td>71%</td><td>71%</td></tr><tr><td>MRP1 (Abcc1 gene)</td><td>DpLoc105347802.1.2</td><td>DbLoc105347802</td><td>Pomacea canaliculata</td><td>0.0</td><td>0.0</td><td>69%</td><td>69%</td></tr><tr><td>MDR1 (Abcb1 gene)</td><td>DpLoc101858982.2.4</td><td>DbLoc101858982.2.2</td><td>Ruditapes philippinarum</td><td>0.0</td><td>0.0</td><td>68%</td><td>70%</td></tr></tbody></table>

of numerous sequence variants (insertion, deletion and single-nucleotide polymorphism), since studied populations are established for a long time and may have been prone to mutation.

In ecotoxicology, biomarkers implied in responses to oxidative stress, heat shock or xenobiotic exposure are used as indicators of stress. Nowadays, multibiomarker approaches have become a standard, considering that the diversity of contaminants that organisms are exposed to can generate a wide diversity of biological responses. Here, we first wanted to focus on the most frequently used biomarkers in ecotoxicology and which are listed in Table 8. Biomarkers were screened, manually checked and all were found expressed in our both species digestive glands. This table also shows the closest species hits (*blastn*) and biomarkers selected here exhibit relatively high value of homology with the corresponding sequence in other mollusk, supporting our contig assignments. The high sequence conservation levels observed here may facilitate the development and adaptation of further measurements assays from existing assays in close related species. Interspecies sequence alignments were also investigated indicating a high levels of homologies between the two species (Table 9). In the same way, this may lead to the development of biomarker based assays that, because of sequence similarity, could be cross compatible for *D. polymorpha* and *D. r. bugensis*.
<table><caption>**Table 9.** Interspecies sequence homologies. Coding sequence from *D. polymorpha* and *D. rostriformis bugensis* biomarkers contigs were used for interspecies homology analysis. Length of coding sequences were also indicated (in base pair; For acronym significations, see Table 8).</caption><thead><tr><th rowspan="2">Biomarkers</th><th colspan="2">D. polymorpha</th><th colspan="2">D. rostriformis bugensis</th><th colspan="3">Interspecies Cds homology</th></tr><tr><th>Contig accession</th><th>Cds length</th><th>Contig accession</th><th>Cds length</th><th>% Ident.</th><th>E-val.</th><th>Gap</th></tr></thead><tbody><tr><td>Acetylcholinesterase</td><td>Dp_LOC105324424.2.2</td><td>1,599</td><td>Db_LOC105324424.3.3</td><td>1,644</td><td>91%</td><td>0.0</td><td>0 (0%)</td></tr><tr><td>Beta Actine</td><td>Dp_ACT.1.3</td><td>1,128</td><td>Db_ACT.2.4</td><td>1,128</td><td>98%</td><td>0.0</td><td>0 (0%)</td></tr><tr><td>Catalase</td><td>Dp_LOC105339902</td><td>1,515</td><td>Db_LOC105339902.2.4</td><td>1,515</td><td>91%</td><td>0.0</td><td>0 (0%)</td></tr><tr><td>Superoxide dismutase (Cu-Zn)</td><td>Dp_SODC</td><td>459</td><td>Db_SODC</td><td>459</td><td>91%</td><td>0.0</td><td>0 (0%)</td></tr><tr><td>Superoxide dismutase (Mn)</td><td>Dp_LOC101852344</td><td>621</td><td>Db_LOC101852344</td><td>621</td><td>92%</td><td>0.0</td><td>0 (0%)</td></tr><tr><td>Estrogen receptor</td><td>Dp_LOC105318922.2.2</td><td>1,476</td><td>Db_LOC105318922.1.2</td><td>1,476</td><td>85%</td><td>0.0</td><td>111 (3%)</td></tr><tr><td>GABARAP</td><td>Dp_LOC105335545.1.2</td><td>351</td><td>Db_LOC105335545.1.6</td><td>351</td><td>95%</td><td>7e-164</td><td>0 (0%)</td></tr><tr><td>Glutathione peroxidase 1</td><td>Dp_GPX1.1.2</td><td>426</td><td>Db_GPX1.1.2</td><td>426</td><td>94%</td><td>0.0</td><td>0 (0%)</td></tr><tr><td>Glutathione peroxidase (Se)</td><td>Dp_LOC106070504.1.5</td><td>729</td><td>Db_LOC106070504.6.8</td><td>714</td><td>85%</td><td>0.0</td><td>15 (2%)</td></tr><tr><td>Heat Shock Protein70</td><td>Dp_HSP7D.3.8</td><td>1,959</td><td>Db_HSP7D.4.8</td><td>1,965</td><td>92%</td><td>0.0</td><td>14 (&lt;1%)</td></tr><tr><td>Metallothionein (isoform 1)</td><td>Dp_MT.6.6</td><td>219</td><td>Db_MT.6.6</td><td>219</td><td>100%</td><td>6e-117</td><td>0 (0%)</td></tr><tr><td>mTOR</td><td>Dp_LOC105331599</td><td>7,422</td><td>Db_LOC105331599</td><td>7,410</td><td>92%</td><td>0.0</td><td>12 (&lt;1%)</td></tr><tr><td>Na/K ATPase</td><td>Dp_LOC106058320</td><td>3,093</td><td>Db_LOC106058320.1.2</td><td>3,096</td><td>89%</td><td>0.0</td><td>3 (&lt;1%)</td></tr><tr><td>Succinate dehydrogenase</td><td>Dp_LOC101864456.1.3</td><td>504</td><td>Db_LOC105338659</td><td>504</td><td>93%</td><td>0.0</td><td>0 (0%)</td></tr><tr><td>Acid phosphatase</td><td>Dp_LOTGIDRAFT_139839</td><td>1,293</td><td>Db_contig_19914</td><td>1,098</td><td>90%</td><td>0.0</td><td>5 (0%)</td></tr><tr><td>Thioredoxin reductase 1</td><td>Dp_LOC105322705.2.2</td><td>1,788</td><td>Db_LOC105322705</td><td>1,932</td><td>92%</td><td>0.0</td><td>3 (&lt;1%)</td></tr><tr><td>MRP1 (Abcc1 gene)</td><td>Dp_LOC105347802.1.2</td><td>3,510</td><td>Db_LOC105347802</td><td>4,686</td><td>91%</td><td>0.0</td><td>6 (&lt;1%)</td></tr><tr><td>MDR1 (Abcb1 gene)</td><td>Dp_LOC101858982.2.4</td><td>4,032</td><td>Db_LOC101858982.2.2</td><td>4,023</td><td>89%</td><td>0.0</td><td>15 (&lt;1%)</td></tr></tbody></table>

By providing *Dreissena* sequences through a user-friendly interface, we open the way to further explorations of Dreissenids molecular mechanisms by such biomarker assays development, primer design allowing targeted expression analysis or promising proteogenomics studies when coupling with mass spectrometry analysis.

## Code availability

Parameters to involved softwares tools are described in the following paragraph.

**DRAP (De novo RNA-seq Assembly Pipeline):** version 1.7, code available online at http://www.sigenae.org/drap/quick_start.html.

**BLAST:** version 2.2.26, e-value under 1e-5 parameter.

**RNAmmer:** version 1.2, standard parameters.

**RepeatMasker:** version open-4-0-3, -engine crossmatch -gcalc -species *Crassostrea gigas* parameters.

**InterProScan:** version 4.8, -goterms -pathways parameters.

**BWA:** version 0.7.12, standard parameters, mem algorithm.

**SAMtools view, sort and index programs:** version 1.1, standard parameters.

**SAMtools IdxStats program:** version 1.1, standard parameters.

**SAMtools rmdup program:** version 1.1, standard parameters.

**GATK:** version v3.0-0-g6bad1c6, -glm BOTH parameter.

**BLAT:** version 35 × 1, standard parameters.

**BUSCO:** version 3.0.2 using the metazoa OrthoDB (v. 9), standard parameters.

Received: 16 May 2019; Accepted: 11 September 2019;

Published online: 25 October 2019

## References

1. Karatayev, A. Y., Claudi, R. & Lucy, F. E. History of *Dreissena* research and the ICAIS gateway to aquatic invasions science. *Aquat. Invasions* 7, 1–5 (2012).

2. Marescaux, J. et al. Sympatric *Dreissena* species in the Meuse River: towards a dominance shift from zebra to quagga mussels. *Aquat. Invasions*. 10, 287–298 (2015).

3. Akopian, M., Garnier, J., Testard, P. & Ficht, A. Estimating the benthic population of *Dreissena polymorpha* and its impact in the lower Seine River, France. *Estuaries* 24, 1003–1014 (2001).

4. Devin, S., Bollache, L., Noël, P.-Y. & Beisel, J.-N. Patterns of biological invasions in french freshwater systems by non-indigenous macroinvertebrates. *Hydrobiologia*. 551, 137–146 (2005).

5. Pimentel, D., Zuniga, R. & Morrison, D. Update on the environmental and economic costs associated with alien-invasive species in the United States. *Ecol. Econ.* 52, 273–288 (2005).

6. Binelli, A., Della Torre, C., Magni, S. & Parolini, M. Does zebra mussel (*Dreissena polymorpha*) represent the freshwater counterpart of *Mytilus* in ecotoxicological studies? A critical review. *Environ. Pollut.* 196, 386–403 (2015).

7. Mills, E. L., Roseman, E. F., Rutzke, M., Gutenmann, W. H. & Lisk, D. J. Contaminant and nutrient element levels in soft tissues of zebra and quagga mussels from waters of southern Lake Ontario. *Chemosphere*. 27, 1465–1473 (1993).
8. Bourgeault, A., Gourlay-Francé, C., Priadi, C., Ayrault, S. & Tlusseau-Vuillemin, M.-H. Bioavailability of particulate metal to zebra mussels: biodynamic modelling shows that assimilation efficiencies are site specific. *Environ. Pollut.* **159**, 3381–3389 (2011).

9. De Lafontaine, Y. *et al*. Biomarkers in zebra mussels (*Dreissena polymorpha*) for the assessment and monitoring of water quality of the St Lawrence River (Canada). *Aquat. Toxicol.* **50**, 51–71 (2000).

10. Faria, M. *et al*. Separating natural from anthropogenic causes of impairment in zebra mussel (*Dreissena polymorpha*) populations living across a pollution gradient. *Aquat. Toxicol.* **152**, 8295 (2014).

11. Palais, F. *et al*. One-year monitoring of core biomarker and digestive enzyme responses in transplanted zebra mussels (*Dreissena polymorpha*). *Ecotoxicology*. **21**, 888–905 (2012).

12. Viarengo, A., Lowe, D., Bolognesi, C., Fabbri, E. & Koehler, A. The use of biomarkers in biomonitoring: a 2-tier approach assessing the level of pollutant-induced stress syndrome in sentinel organisms. *Comp. Biochem. Physiol. Part C Toxicol. Pharmacol.* **146**, 281–300 (2007).

13. Albrecht, C., Schultheiß, R., Kevrekidis, T., Streit, B. & Wilke, T. Invaders or endemics? Molecular phylogenetics, biogeography and systematics of *Dreissena* in the Balkans. *Freshw. Biol.* **52**, 1525–1536 (2007).

14. Stepien, C. A., Taylor, C. D. & Dabrowska, K. A. Genetic variability and phylogeographical patterns of a nonindigenous species invasion: a comparison of exotic vs native zebra and quagga mussel populations. *J. Evol. Biol* **15**, 314–328 (2002).

15. Penarrubia, L., Sanz, N., Pla, C., Vidal, O. & Vinas, J. Using massive parallel sequencing for the development, validation, and application of population genetics markers in the invasive bivalve zebra mussel (*Dreissena polymorpha*). *PLoS One*. **10**, e0120732 (2015).

16. Gonzalez, V. L. *et al*. A phylogenetic backbone for bivalvia: an RNA-seq approach. *Proc. Biol. Sci.* **282**, 20142332 (2015).

17. Soroka, M. *et al*. Next-generation sequencing of *Dreissena polymorpha* transcriptome sheds light on its mitochondrial DNA. *Hydrobiologia*. **810**, 255–263 (2018).

18. Therriault, T. W., Docker, M. F., Orlova, M. I., Heath, D. D. & MacIsaac, H. J. Molecular resolution of the family Dreissenidae (Mollusca: Bivalvia) with emphasis on Ponto-Caspian species, including first report of *Mytilopsis leucophaeata* in the Black Sea basin. *Mol. Phylogenet. Evol.* **30**, 479–489 (2004).

19. Geda, S. R. *et al*. Multilocus phylogeny of the zebra mussel family Dreissenidae (Mollusca: Bivalvia) reveals a fourth neotropical genus sister to all other genera. *Mol. Phylogenet. Evol.* **127**, 1020–1033 (2018).

20. Hanana, H., Turcotte, P., Andre, C., Gagnon, C. & Gagne, F. Comparative study of the effects of gadolinium chloride and gadolinium-based magnetic resonance imaging contrast agent on freshwater mussel, *Dreissena polymorpha*. *Chemosphere*. **181**, 197–207 (2017).

21. Bultelle, F., Panchout, M., Leboulenger, F. & Danger, J. M. Identification of differentially expressed genes in *Dreissena polymorpha* exposed to contaminants. *Mar. Environ. Res.* **54**, 385–389 (2002).

22. Cabau, C. *et al*. Compacting and correcting Trinity and Oases RNA-Seq de novo assemblies. *PeerJ*. **5**, e2988 (2017).

23. Schulz, M. H., Zerbino, D. R., Vingron, M. & Birney, E. Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. *Bioinformatics*. **28**, 1086–1092 (2012).

24. Lagesen, K. *et al*. RNAmmer: consistent and rapid annotation of ribosomal RNA genes. *Nucleic Acids Res.* **35**, 3100–3108 (2007).

25. Jones, P. *et al*. InterProScan 5: genome-scale protein function classification. *Bioinformatics*. **30**, 1236–1240 (2014).

26. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics*. **25**, 1754–1760 (2009).

27. Li, H. *et al*. The sequence alignment/map format and SAMtools. *Bioinformatics*. **25**, 2078–2079 (2009).

28. McKenna, A. *et al*. The Genome Analysis Toolkit: a MapReduce framework for analyzing next- generation DNA sequencing data. *Genome res.* **20**, 1297–1303 (2010).

29. Mariette, J. *et al*. RNABrowse: RNA-Seq de novo assembly results browser. *PLoS One*. **9**, e96821 (2014).

30. NCBI Sequence Read Archive, https://identifiers.org/ncbi/insdc.sra:SRP173947 (2018).

31. GenBank, https://identifiers.org/insdc:GHIW00000000 (2018).

32. GenBank, https://identifiers.org/insdc:GHIX00000000 (2018).

33. NGSPipeline interface, http://ngspipelines.toulouse.inra.fr:9014 (2017)

34. Péden R. *et al*. Environmental transcriptomes of invasive dreissena: a move towards an upcoming model species. figshare, https://doi.org/10.6084/m9.figshare.c.4483856 (2019).

35. Kent, W. J. BLAT—the BLAST-like alignment tool. *Genome res.* **12**, 656–664 (2002).

36. Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. *Bioinformatics*. **31**, 3210–3212 (2015).

## Acknowledgements

This work was supported by the BIOMICS program (Biomarkers and Omics in the freshwater Bivalve *Dreissena polymorpha*) of the French Agency for Biodiversity/National Agency for Water and Aquatic Environments (AFB/ONEMA) Authors would like to thanks the GeT-PlaGe for performing sequencing and Genotoul platform for transcriptome assembly, annotation and validation.

## Author contributions

B.S., S.D., P.P., L.G. and S.P.-D. conceived and designed the study. B.S., J.F., S.P.-D., M.P. and F.L. realized experiments. C.K. and R.-F.S. performed sequencing and transcriptome assembly and annotation. B.S., R.P. and P.P. proceeded to data analysis. R.P., B.S., S.D., R.-F.S. and C.K. wrote the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41597-019-0252-x.

Correspondence and requests for materials should be addressed to S.D.

Reprints and permissions information is available at www.nature.com/reprints.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_4.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or
format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Cre-
ative Commons license, and indicate if changes were made. The images or other third party material in this
article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the
material. If material is not included in the article's Creative Commons license and your intended use is not per-
mitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the
copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

The Creative Commons Public Domain Dedication waiver http://creativecommons.org/publicdomain/zero/1.0/
applies to the metadata files associated with this article.

© The Author(s) 2019