DOI: 10.5061/dryad.r7sqv9sdz
title: Genome-wide signatures of synergistic epistasis during parallel adaptation in a Baltic Sea copepod

description: The role of epistasis in adaptive evolution has remained an unresolved
 problem dating back to the Evolutionary Synthesis. This role is now being
 revisited due to its relevance for polygenic adaptation. In the absence of
 epistasis, polygenic adaptation is predicted to result in non-parallel
 evolution, because repeated selection could act on subsets of effectively
 redundant alleles. However, positive epistatic interactions among adaptive
 alleles would make the alleles non-redundant and selection for particular
 allelic combinations could drive parallel evolution. The inability to
 address this fundamental question might arise from traditional approaches
 lacking the power to capture the genomic architecture and dynamics of
 polygenic adaptation. To address this problem, we employed a replicated
 and controlled evolution experiment using the copepod Eurytemora
 affinis to elucidate the evolutionary response architecture to
 rapid salinity decline, a predicted consequence of global climate change
 in higher latitudes. Based on time-resolved pooled whole-genome
 sequencing, we uncovered a remarkably parallel response, despite polygenic
 adaptation involving over 1000 loci across ten replicate selection lines.
 Interestingly, single-nucleotide polymorphism (SNP) frequencies converged
 during the experiment, far beyond expectations, resulting in replicate
 lines sharing 93.1% of selected alleles. Using simulations, we found that
 this polygenic parallelism was consistent with synergistic epistasis among
 alleles responding in concert across replicate lines, a phenomenon that
 may be common for selection on complex physiological traits. Furthermore,
 we found that the same SNPs with signatures of selection in the laboratory
 also exhibited signatures of selection across a natural salinity gradient
 in the Baltic Sea. Our study provides the first experimental evidence that
 polygenic adaptation can actually be highly repeatable at the genomic
 level, given the presence of synergistic epistasis among the loci under
 selection.
description: A draft genome for <em>Eurytemora
 affinis</em> complex was constructed from long read and long-range
 sequencing technology. To generate genomic data for assembly, an inbred
 line was generated from 30 generations of full-sibling mating of copepods
 derived from a saline population in Baie de L’Isle Verte, St. Lawrence
 marsh, Quebec, Canada (Atlantic clade, aka <em>E.
 carolleae</em>). Prior to DNA extraction, the culture was treated
 with a series of antibiotics to reduce bacterial contamination, including
 Primaxin (20mg/l), Voriconazole (0.5 mg/l for at least 2 weeks prior to
 DNA extraction), D-amino acids to reduce biofilm (10 mM D-methionine,
 D-tryptophan, D-leucine, and 5 mM D-tyrosine, for at least for 2 weeks
 prior to DNA extraction). DNA was extracted from pooled adult copepods and
 used to generate Pacific Biosciences (PacBio) long-read data (2.65 million
 subreads, 30.2 Gb, read N50 = 19.25 kb), Dovetail
 Chicago<sup>® </sup>(213 million 150 bp read pairs) and Hi-C
 (171 million 150 bp read pairs) Illumina HiSeq X data.
 PacBio reads were assembled into contigs using
 <em>wtdbg </em>v2.5 and polished with Racon v1.4.3. This
 procedure resulted in an assembly of size 575 Mb in 3013 contigs with an
 N50 of 951 kb, an N90 of 106 kb, and maximum contig length of 7.7 Mb.
 Contigs were scaffolded with HiRise<sup>TM</sup> (Dovetail
 Genomics) using Dovetail Chicago<sup>® </sup>and Hi-C data.
 This procedure resulted in a highly contiguous assembly with 90% of the
 genome present in four scaffolds (1706 total scaffolds, 0.02% gaps),
 likely corresponding to the four chromosomes of <em>E.
 affinis</em> observed in karyotype (data not shown). These 4
 scaffolds consisted of 518 Mb of sequence. Gene annotations were generated
 by mapping annotations from the previously published, short-read-assembled
 genome of <em>E. affinis</em> complex using LiftOff
 v1.6.1. European <em>E.
 affinis</em> populations, the subject of the laboratory selection
 experiment, are highly genetically divergent from North American
 populations of the <em>E. affinis</em> complex, from which the
 draft genome was derived. Therefore, to maximize the number and accuracy
 of SNP calls, a pseudo-reference genome was assembled for SNP calling
 using an approach used in a previous study. This approach uses a reference
 transcriptome as an anchor to assemble Pool-seq data around coding
 regions, generating a reference assembly that includes coding sequences
 and surrounding sequence with putative regulatory function.
 To generate RNA-seq data for assembly of the reference
 transcriptome, approximately 100 individuals of all life stages were
 collected from the European <em>E. affinis</em> laboratory
 culture derived from Kiel, Germany. Animals were pooled and subjected to
 a Trizol/Qiagen RNeasy hybrid total RNA extraction protocol. An mRNA
 sequencing library was prepared using the Illumina TruSeq Stranded mRNA
 kit and sequenced on an Illumina HiSeq 4000 sequencer, generating
 approximately 186 million 100 bp paired-end reads. Raw RNA-seq reads were
 processed using BBDuk in the BBtools package
 (https://jgi.doe.gov/data-and-tools/bbtools) to filter adapter sequences,
 low complexity sequences, and low quality (Q &lt; 10) bases using a
 sliding window (ktrim=r k=23 mink=11 hdist=1 qtrim=w trimq=10 minlen=36
 entropy=0.01 entropywindow=50 entropyk=5 tbo). To
 filter sequences that may have originated from microbial contamination
 prior to assembly, reads were mapped against a database of reference and
 representative bacterial, archaeal, and fungal genomes from the NCBI
 RefSeq using Bowtie 2 v2.3.5 and were removed if they mapped concordantly.
 Cleaned reads were assembled using Trinity v2.6.6 in strand-specific mode.
 mRNA expression estimates were made using RSEM v1.3.1 with cleaned reads
 mapped to the Trinity assembly with Bowtie 2. Only the most highly
 expressed Trinity ‘isoform’ per ‘gene’ was retained to obtain the most
 highly supported transcriptome assembly. Transcript sequences were
 clustered to 95% similarity using CD-HIT v4.7 to reduce redundancy
 associated with allelic variation and assembly errors. Transdecoder v5.5
 was used to predict open reading frames and coding sequences. Predictions
 included <em>blastp </em>hits to annotated proteins of
 the <em>E. affinis </em>complex draft genome and HMMER
 v3.2.1 (http://hmmer.org) hits to the PFAM protein database.
 The resulting non-redundant transcriptome assembly was used to
 anchor the de-novo assembly of Pool-seq data around coding regions to
 obtain a pseudo-reference genome. First, the ‘left’ pairs of Pool-seq
 reads from the starting laboratory population were mapped to the coding
 sequences of the transcriptome with BWA-MEM v0.7.17, retaining reads with
 a mapping quality &gt; 20. The corresponding ‘right’ pairs of the
 mapped reads were extracted and the mapped reads were assembled using
 Trinity v. 2.6.6, as Trinity was developed to assemble highly variable
 sequences with potentially uneven coverage. The resulting assembly was
 clustered to 95% similarity using CD-HIT v4.7. The above mapping and
 assembly procedures were repeated one additional time, mapping reads to
 the Trinity-assembled Pool-seq data rather than to the transcriptome. In
 this way, the assembly was extended into genomic regions surround the
 transcriptome-derived coding sequences. Finally, only assembled contigs
 with a significant BLAST hit (E-value &lt; 0.001) to the draft genome
 (see above) were retained to further eliminate assembly errors. The
 pseudo-reference genome ultimately spanned a greater portion of the genome
 with greater contiguity than the transcriptome (Supplementary Table 10).
 These pseudo-reference genome contigs were arranged and oriented into
 scaffolds by BLASTing (<em>blastn</em> E-value &lt; 0.001)
 contigs to the long-read genome assembly and retaining the top hit per
 contig.