Candidate genes for shell colour polymorphism in Cepaea nemoralis

Sample collection, mRNA extraction, sequencing and assembly

Four juvenile C. nemoralis with different phenotypes (two with brown shell and no banding and two with a yellow shell with multiple dark bands; Fig. S1) were collected at the Van Veldhuizenbos near Dronten, The Netherlands. Total RNA was extracted separately from the mantle and the foot for each of the four C. nemoralis individuals using the NucleoSpin RNA kit (Macherey-Nagel, Düren, Germany), following the manufacturer’s protocol. The remains of the specimens have been deposited in the alcohol collection of Naturalis Biodiversity Center under reference numbers RMNH.5004228-5004231. The total RNA was subjected to polyA selection, converted to cDNA and used to generate sequencing libraries as described in Salazar-Jaramillo et al. (2017). The libraries were paired-end 2 × 100 bp sequenced on an Illumina HiSeq 2000 at the Leiden Genome Technology Center.

De novo assembly, protein prediction and annotation

Sequence reads of all eight samples were combined to create a reference transcriptome assembly using Trinity v2.0.5 with default settings (Grabherr et al., 2011). To reduce redundancy, transcripts were clustered using CD-HIT-Est (Fu et al., 2012) at a cut-off of 95% identity. To remove contamination, we used NCBI’s VecScreen. This tool searches for segments that match sequences in the UniVec database. We used BUSCO (Simao et al., 2015) to assess the completeness of the transcriptome. This tool searches the transcriptome for a reference set of single-copy orthologs. The metazoan set (odb9) was used as a reference set. All genes in the transcriptome were annotated for gene ontology (GO) using the Trinotate part of the Trinity package. Trinotate uses blastp and blastx to compare the predicted peptide sequences and the reference transcriptome to the uniprot_sprot and the uniprot_uniref peptide databases, performs an hmm-scan on the Pfam-A database and assigns GO terms.

Differential expression

To identify genes that were overexpressed in the shell-forming mantle tissue, the Trinity RNA-seq pipeline (Grabherr et al., 2011) was used. First, the reads for both organs in all four individuals were mapped to the reference transcriptome. Second, the mapped reads were counted to visualize the expression of these transcripts, using the estimation method eXpress (Roberts & Pachter, 2013). A high count of mapped reads indicates a high expression rate, while a low count of mapped reads indicates a low expression rate. With the estimated counts of reads, expression profiles were generated using the R-package EdgeR (Robinson, McCarthy & Smyth, 2010). Normalization took place as part of the EdgeR workflow. The profiles were then filtered on fold change and false discovery rate (FDR) with the Trinity default cut-off scores of 4 and 0.001 respectively. The genes that were overexpressed in the mantle were selected for further analysis. A heatmap of differentially expressed genes versus samples was produced using the script analyze_diff_expr.pl within the Trinity package.

Sequence variants

To identify consistent sequence differences between the two phenotypes (yellow/banded and the brown/unbanded), we conducted variant calling following the following protocol. First, reads were mapped to the reference transcriptome using Bowtie2 (Langmead & Salzberg, 2012). Next, GATK’s HaplotypeCaller (McKenna et al., 2010) was run to search for single-nucleotide polymorphisms (SNPs). The default parameters were used with a stand_call confidence of 20.0 and a stand_emit confidence of 20.0. The VariantFiltration tool was used with default parameters to filter out false-positive variants. Clusters were generated when 3 SNPs were present within a window of 35 bases. Variants with a Fisher-strand Score (FS) greater than 30.0 or a quality by depth (QD) value less than 2.0 were filtered out. SNPs with sequencing depth less than 10 in any of the samples were removed. The effect of each SNP on the transcript product was predicted using SnpEff (Cingolani et al., 2012).

The resulting set of SNPs was filtered for consistency with the phenotypes using GATK and custom R-scripts. The inheritance of the polymorphism is well understood (Cain & Sheppard, 1954; Murray, 1963; Jones, Leith & Rawlings, 1977). The brown ground colour is dominant over the yellow ground colour and absence of banding is dominant over the presence of multiple bands. The candidate SNPs were therefore filtered so that the brown, unbanded snails were either heterozygous or homozygous and the yellow, banded snails were homozygous. Mantle-enriched transcripts contained at least one phenotype-consistent sequence variant were blastn-searched (Altschul et al., 1990) against the non-redundant nucleotide database.

RAD-tags

Richards et al. (2013) identified eleven anonymous RAD tag sequences that lie near the colour-polymorphism supergene. In the original research, a cut-off of 96 base pairs was used for the tags. Because of this, the tags are shorter than the full assembly of the RAD tag reads. The overlapping parts of the RAD tag reads were used to generate a consensus for each tag that was longer than the original consensus RAD tag. These extended tags were blast-searched against the reference transcriptome using the megablast algorithm within the Blast+ (Camacho et al., 2009) command line tool. Hits were filtered by E-value less than 10⁻¹⁰.

Transcriptome assembly and annotation

Eight RNAseq libraries (two for each individual) were constructed and sequenced, yielding between 63.9 × 10⁶ and 139.5 × 10⁶ paired reads per sample (Table 1). These reads were assembled into contigs representing 171,051 putative transcripts deriving from 33,109 putative genes (Table 2). The frequency of the number of transcripts per gene is shown in Fig. S2. Clustering using CD-HIT reduced the number of transcripts to 150,380. On average, 71% of the raw sequence reads aligned at least once to the reference transcriptome (Table S1). After removing vector contamination, 147,397 transcripts remained.

The completeness of the transcriptome was assessed using BUSCO. Of the 978 single-copy orthologs that were searched for, 920 were found completely. Of these, 765 (78.2%) were found single-copy, 155 (15.8%) orthologs were found duplicated. A further 37 (3.8%) of the orthologs were found fragmented and 21 (2.2%) were not found in the transcriptome.

A total of 111,416 (75.6%) transcripts were annotated by Trinotate. Protein sequence comparison to the uniprot_sprot database yielded annotations for 89,386(60.6%) transcripts, to the uniprot_uniref database for 87,045 (59.1%) transcripts and to the Pfam database for 30,405 (20.6%) transcripts. After annotation, Trinotate assigned GO-terms to 19,658 genes, covering 26,039 transcripts (17.7%; Data S1).

Differential expression

Hierarchical clustering of the overall expression data clearly separated the mantle and foot tissue samples (Fig. 1). EdgeR identified 1,961 transcripts as upregulated and 1,260 as downregulated in the mantle relative to the foot (Fig. 2, Data S2).

Sequence polymorphisms

A total of 73,817 SNPs passed our filtering steps (Table S2). This included 12,273 synonymous and 12,451non-synonymous variants. The remaining 49,092 SNPs fell outside predicted open reading frames. A total of 569 SNPs were found in transcripts that were overexpressed in the mantle and showed phenotype-consistent allelic variation (Table 3). A total of 197 mantle-enriched transcripts contained at least one phenotype-consistent sequence variant. This set was subjected to more detailed annotation. We obtained database matches for 98 transcripts, of which 38 were functionally annotated (Table S3). Annotations that indicate putative functions related to shell formation are summarized in Table 4. Prominent among these annotations are sequences that are thought to play a role in the production of the nacreous layer in mollusc shells, including dermatopin and nacroperlin genes (marked in Table 4).

RAD-tags

A set of 12 transcripts had a match with at least one of the elongated tags that lie in close proximity to the supergene (Richards et al., 2013) (Table S4). None of these transcripts were found in our list of differentially expressed transcripts with phenotype-consistent variants. Eight of the twelve transcripts (matching two of the RAD tags closest to the supergene) had a blast hit to the nr and nt databases (Table S4, Data S3). All of these hits were the Helix pomatia homolog of metallothionein. This same sequence was also found in the variant analysis (six of 197 transcripts = 3%). None of the transcripts in these two sets overlap. The percentage of metallothionein hits among the entire reference transcriptome was 6,1% (2,506 of 40,748 transcripts), suggesting that this gene is not overrepresented in our variant analysis.