DOI: 10.5061/dryad.76hdr7t3s
title: Data from: Temperature-dependent gene regulatory divergence underlies local adaptation with gene flow in the Atlantic silverside

description: Gene regulatory divergence is thought to play an important role in
 adaptation, yet its extent and underlying mechanisms remain largely
 elusive for local adaptation with gene flow. Local adaptation is
 widespread in marine species despite generally high connectivity and is
 often associated with tightly linked genomic architectures, such as
 chromosomal inversions. To investigate gene regulatory evolution under
 gene flow and the role of inversions associated with local adaptation to a
 steep thermal gradient, we generated RNA-seq data from Atlantic
 silversides (Menidia menidia) from two locally adapted populations and
 their F1 hybrids, reared under two temperatures. We found substantial
 divergence in gene expression and thermal plasticity between populations,
 with up to 31% of genes being differentially expressed. Reduced thermal
 plasticity, temperature-dependent gene misexpression and the disruption of
 co-expression networks in hybrids point towards a role of regulatory
 incompatibilities in local adaptation, particularly under colder
 temperatures. Chromosomal inversions show an accumulation of regulatory
 incompatibilities but are not consistently enriched for differentially
 expressed genes. Together, these results suggest that gene regulation can
 diverge substantially among populations despite gene flow, partly due to
 the accumulation of temperature-dependent regulatory incompatibilities
 within inversions.
description: To improve the contiguity of the Atlantic silverside reference
 genome (Tigano et al.
 2021), we anchored the genome assembly to a RAD-seq based female
 Georgia linkage map (Akopyan et al.
 2022). However, the linkage map used for this assembly differs
 marginally from the one in Akopyan et al.
 (2022),
 as it was constructed from genome-aligned RAD-seq data rather than
 <em>de novo</em> assembled loci. Overall, both linkage maps
 are highly comparable. The linkage map was constructed as described in
 Akopyan et
 al. (2022) with slight modifications outlined here:
             In brief, we used ddRAD sequencing (Peterson et al.
 2012) to identify and genotype single nucleotide polymorphisms (SNPs) for
 linkage map construction from ​​568 individuals across five families,
 including the two founders, 138 F1 offspring, six additional
 F<sub>1</sub> siblings and their 282 F<sub>2</sub>
 offspring. Reads were processed in Stacks v1.48 (Catchen et al. 2013) with
 the module <em>process_radtags</em> to discard low-quality
 reads and reads with ambiguous barcodes or RAD cut-sites. The remaining
 reads were demultiplexed and aligned to the <em>Menidia
 menidia</em> reference genome v1 (Tigano et al.
 2021) using Bowtie2 v2.2.9
 (<em>-very-sensitive</em>). We only retained those reads that
 were uniquely mapped to the reference genome and extracted RAD loci with:
 i) minimum read depth of three, ii) minimum mapping quality of 10, iii)
 and maximum clipped proportion of 0.15<em>. </em>Variant
 calling was also performed with <em>pstacks</em> using the
 default SNP model with a genotype likelihood ratio test critical value (α)
 of 0.05. We built a catalog of all loci using parents (and grandparents
 for the F<sub>2</sub> generation) with
 <em>cstacks</em>, and matched progeny against the catalog
 using <em>sstacks</em>. The <em>populations</em>
 module was used to filter variants to retain only the first SNP per locus
 and generated a VCF file for each of the two F<sub>1</sub>
 families, and one for the F<sub>2</sub> generation including
 the three intercross families. We constructed one
 female linkage map for each of our three crosses (F1 GAxNY, F1 NYxGA, F2)
 using <em>Lep-MAP3</em> (Rastas
 2017). In brief, offspring genotypes were called by accounting
 genotype information of parents (and grandparents in
 F<sub>2</sub> family) with the
 <em>ParentCall2</em> module, markers with high segregation
 distortion were removed using the distortionLod=1 option in
 <em>SeparateChromosomes2</em>, separated markers were merged
 into linkage groups with a logarithm of odds (LOD) score limit of 20 and
 minimum linkage group size of 10 markers using markers informative in
 females only, and we used the <em>OrderMarkers2</em> module to
 compute genetic distances in centimorgan (i.e., recombination rates)
 between all adjacent markers for each linkage group using the default
 Haldane’s mapping function. We used maternally informative markers to
 construct the F<sub>1</sub> maps, and both maternally and
 dually informative markers to construct the F<sub>2</sub>
 map. We used the female F<sub>1</sub>
 linkage map for the Georgia population to anchor and order the Atlantic
 silverside reference genome v1 scaffolds into chromosomes using
 <em>AllMaps</em> (Tang et al.
 2015). Chromosomes were renamed based on synteny with the medaka
 genome. Furthermore, we converted the coordinates of RAD loci for all
 linkage maps from scaffold to anchored chromosome coordinates using
 <em>CrossMap v.0.1.4</em> (Zhao et al.
 2014) and identified inversions between NY and GA by comparing
 the F1 NY-linkage map and F2 linkage map to the GA-anchored reference
 genome (see (Akopyan et al.
 2022) for details).             Lastly, we
 re-annotated the anchored genome assembly (<em>M. menidia</em>
 reference genome v2) following the pipeline outlined in (Tigano et al.
 2021). In brief, first we identified and annotated repeats using
 <em>Repeatmodeler2</em> (Flynn et al.
 2020) and <em>Repeatmasker</em>  (Smit et al.
 2015). Next, we used <em>BRAKER2</em> (Brůna et al.
 2021) with evidence from: i) RNA-seq for diverse Atlantic
 silverside individuals from different populations and at different
 developmental stages (Therkildsen
 and Baumann 2020; Tigano et al. 2021); ii) protein-homology
 evidence from six different teleost species and the UniProtKB (Swiss-prot)
 database. Subsequently, we performed five iterative rounds in
 <em>MAKER </em>as described in (Tigano et al.
 2021). Lastly, functionally annotated the predicted genes using
 Blast2GO in Omnibox v.1.2.4 (Gotz et al.
 2008) using the UniProtKB (Swiss-Prot) database and
 InterProScan2 (Zdobnov and
 Apweiler 2001).
description: # Temperature-dependent gene regulatory divergence underlies local
 adaptation with gene flow in the Atlantic silverside
 [https://doi.org/10.5061/dryad.76hdr7t3s](https://doi.org/10.5061/dryad.76hdr7t3s) ## Description of the data and file structure 1. Mmenidia_refgenome_anchored.all_renamed_v2.fasta: Reference genome scaffolds for Menidia menidia anchored to a consensus linkage map. See Methods section and Supplementary material for detailed methods. 2. mme_annotation_anchored_genome_final_clean.noseq.gff: Gene annotation for the anchored reference genome. Output from Maker. 3. mme_protein_blast2go_annot_topGO.annot: Gene ontology annotation (from Blast2GO) used for GO overrepresentation analysis using topGO 4. mme_blast2go_annotation_descriptions.txt: Functional annotation (from Blast2GO) of genes annotated in Menidia menidia genome.