DOI: 10.5061/dryad.573n5tb7r
title: Transcriptome data: salinity adaptation in Rhithropanopeus harrisii across an estuarine gradient

description: Rhithropanopeus harrisii is a common estuarine crab native to the East and
 Gulf Coasts of North America. Here, it is found along a broad range of
 salinities, spanning from ~1 PSU to ~25 PSU along estuarine gradients. As
 part of a larger study on the species' potential use of low-salinity
 refuges from parasitism, we tested for differences in gene expression with
 salinity using crabs from three distinct estuarine reaches along the
 Pamlico River in North Carolina. We acclimated crabs collected at natural
 sites with 1, 5, and 11 PSU mean salinity to 0.8, 3, and 15 PSU in the
 laboratory for one week and sequenced gill tissue (second posterior gill)
 from each group using mRNA-Seq. Data in this repository includes
 information on sequenced samples (*csv), the de novo transcriptome
 assembly (two versions: with and without expression filtering; *.fasta),
 transcriptome annotation from EnTAP (two files: using the runN and runP
 functions; *.tsv), a spreadsheet of normalized read count data from salmon
 pseudo-alignment (*.tsv), high-quality SNPs identified from the
 transcriptome sequencing with GATK (*.vcf), and three custom scripts used
 in processing the SNP data (*.py and *.R). Raw sequence data is archived
 in GenBank's SRA: BioProject PRJNA730785, BioSamples
 SAMN19241288-SAMN19241333.
description: <strong>Methods reproduced from the associated paper,
 "Invasion of the body snatchers: the role of parasite introduction in
 host distribution and response to salinity in invaded
 estuaries".</strong> <em>Salinity
 Exposure Experiment</em><br> Adult R. harrisii were collected
 in January 2019 from three populations in the Pamlico River: the Estuarium
 (~1 PSU mean salinity), Mallard Creek (~5 PSU mean salinity), and
 Wright's Creek (~11 PSU mean salinity). Crabs were acclimated for
 three days at their collection salinities and held at 20°C. 20 crabs from
 Mallard Creek and Wright's Creek were randomly divided into two
 groups and acclimated to 3 or 15 PSU. Crabs were ramped ±2.5 PSU every 24h
 until reaching final experimental salinities. For the Estuarium, 10 crabs
 each were ramped to 3 and 15 PSU, with a third group initiated at 0.0 PSU
 and transitioned to 0.8 PSU after 5d. After reaching experimental
 salinities, crabs were maintained for 7d in individual chambers of an
 18-well polycarbonate parts-box. Crab position was randomized across
 treatments, and boxes were rotated daily. Crabs were fed every 3d
 (carnivore pellets), and water was changed 24h post-feeding. Following the
 acclimation period, crabs were dissected and the second posterior
 (primarily osmoregulatory) gill from the crab’s left and right sides were
 placed into cryovials containing 1ml RNALater. Preserved gills were stored
 at -80°C until processing. <em>Sample sequencing,
 transcriptome assembly, and mapping for gene
 expression</em><br> RNA was extracted from stabilized gills
 using a modified version of the TRI reagent protocol (Simms et al. 1991),
 with the addition of 10 mg of glycogen to 500 uL of the aqueous phase
 before precipitation followed by 75% ethanol wash steps. After extraction,
 47 samples yielded sufficient DNA to proceed with cDNA library
 construction using an Illumina TruSeq stranded mRNA sequencing kit. Forty
 samples were sequenced across two lanes of 50 bp single-read sequencing on
 an Illumina HiSeq4000 at the Vincent J. Coates Genomics Sequencing
 Laboratory at UC Berkeley; seven samples (randomized across groups) were
 sequenced with samples for another project on a single lane of 150 bp
 paired-end sequencing at Genewiz (South Plainfield, NJ, USA).
 In total, 46 samples yielded good sequence data and were used in
 the genomic analyses, with 4-9 individual replicates per treatment group.
 Raw reads were trimmed with Trim galore! V0.6.4_dev, using cutadapt v2.6,
 to remove bases with Phred &lt;20 and ≥3 bp adaptor overlap with a
 maximum trimming error rate of 0.1
 (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/; Martin
 2011). Reads were discarded if, after trimming, they were &lt;20 bp
 (50 bp single-end reads) or &lt;50 bp (150 bp paired-end reads). We
 assembled a <em>de novo</em> transcriptome with Trinity
 v2.11.0 (Haas et al. 2013), using the single-end reads from all 46 samples
 including the forward reads from the samples sequenced with 150 bp
 paired-end sequencing. For the Trinity assembly, we set the minimum contig
 length to 200, normalized coverage to a maximum of 50 reads, and set the
 minimum kmer coverage to 2. The resulting <em>de novo
 </em>transcriptome had 140,196 contigs and 214,737 isoforms; we
 retained only the longest isoform per gene for downstream
 analysis. The <em>de novo</em>
 transcriptome was annotated using EnTAP v0.9.1 with Uniprot’s swissprot
 and TrEMBL protein databases and NCBI’s non-redundant nucleotide database
 as references (Hart et al. 2020). EnTAP was run twice, with both blastp
 (e-value ≤10-5) and blastx (e-value ≤10-3) alignments. Between the two
 approaches, we were able to annotate 16,894 contigs with a sequence match,
 and an additional 5,613 contigs with a gene family match. Of the contigs
 we could annotate with EnTAP, 1,240 were identified as non-crab
 contamination (e.g. matches to fungi, plants, bacteria, etc). To
 additionally screen for potential contamination by <em>Loxothylacus
 panopaei</em>, a parasitic barnacle which grows tissue throughout
 the host’s body, we ran a local blastn search with BLAST+ against a
 <em>de novo</em> transcriptome for <em>L. panopaei
 </em>(e-value ≤10-3; Camacho et al. 2009; Tobias et al., in review).
 To identify non-messenger RNA, we ran another local blastn search against
 the SILVA rRNA representative database (e-value ≤10-3; Quast et al. 2012).
 We identified 134 and 206 contigs as likely <em>L.
 panopaei</em> and rRNA contamination, respectively. We removed all
 1,467 contigs identified as contamination with any approach from the
 transcriptome, leaving a <em>de novo R. harrisii</em> gill
 transcriptome of 138,729 contigs. For all downstream
 analyses, paired-end reads were converted to 50 bp single-end reads by
 taking only the forward reads and cutting them to 50 bp before trimming as
 above. This was done to ensure that mapping and subsequent analysis would
 not be biased by sequence length or format.<em> </em>For
 differential gene expression analysis, we mapped individual reads back
 against the full cleaned transcriptome with Salmon v1.3.0 (Patro et al.
 2017). <em>SNP Identification and
 Genoptyping</em><br> For sequence-based analyses, we further
 filtered the transcriptome to remove contigs with low sequence support
 (&lt;1 TPM), leaving 43,908 contigs. We mapped individual sample reads
 back to the transcriptome with Bowtie 2 v2.4.1 (Langmead &
 Salzberg 2013), using the default settings. Reads were mapped against a
 transcriptome that included contigs that were putative outside or rRNA
 contamination to reduce bias; these contigs were removed post-mapping but
 before any subsequent analysis. We also ran blastn searches against the
 full mitochondrial genomes of six Portunidae species (e-value ≤10-3), and
 removed 49 contigs as likely mitochondrial after mapping. SNPs were
 identified and called using the Genome Analysis Toolkit v4.1.8.1 (GATK;
 McKenna et al. 2010; DePristo et al. 2011). With GATK we initially
 selected only biallelic SNPs with a quality ≥20 Phred; this yielded
 417,963 raw SNPs. We further filtered these markers using a custom python
 script which retained only SNPs which had high-quality (≥20 Phred),
 high-coverage (≥5 reads) individual genotypes in at least 9 individuals in
 each of the three sources. We removed one sample from Wright’s Creek from
 downstream sequence analysis because of high genotype missingness
 (&gt;20%). After this process, we had a set of 23,303 high-quality
 SNPs in 45 individuals, with each individual genotyped at ≥84% of these
 SNPs. We further screened this SNP panel to remove SNPs in strong linkage
 disequilibrium (R<sup>2</sup> ≥ 0.8). To do this, we first
 calculated within and between-contig LD values using vcftools v0.1.16
 (Danecek et al. 2011), converted these pairwise R<sup>2</sup>
 values to a matrix using a custom R script, and used the R package
 ‘igraph’ v1.2.5 to identify clusters of SNPs in strong LD (Csárdi
 & Nepusz 2006). All but one SNP in each such cluster was removed,
 leaving a set of 16,592 high-quality, high-coverage, independent SNPs for
 downstream sequence analysis.
description: All necessary information is available in the uploaded
 files.