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.

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 de novo 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₁ siblings and their 282 F₂ offspring. Reads were processed in Stacks v1.48 (Catchen et al. 2013) with the module process_radtags to discard low-quality reads and reads with ambiguous barcodes or RAD cut-sites. The remaining reads were demultiplexed and aligned to the Menidia menidia reference genome v1 (Tigano et al. 2021) using Bowtie2 v2.2.9 (-very-sensitive). 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. Variant calling was also performed with pstacks 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₂ generation) with cstacks, and matched progeny against the catalog using sstacks. The populations 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₁ families, and one for the F₂ generation including the three intercross families.

We constructed one female linkage map for each of our three crosses (F1 GAxNY, F1 NYxGA, F2) using Lep-MAP3 (Rastas 2017). In brief, offspring genotypes were called by accounting genotype information of parents (and grandparents in F₂ family) with the ParentCall2 module, markers with high segregation distortion were removed using the distortionLod=1 option in SeparateChromosomes2, 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 OrderMarkers2 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₁ maps, and both maternally and dually informative markers to construct the F₂ map.

We used the female F₁ linkage map for the Georgia population to anchor and order the Atlantic silverside reference genome v1 scaffolds into chromosomes using AllMaps (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 CrossMap v.0.1.4 (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 (M. menidia reference genome v2) following the pipeline outlined in (Tigano et al. 2021). In brief, first we identified and annotated repeats using Repeatmodeler2 (Flynn et al. 2020) and Repeatmasker  (Smit et al. 2015). Next, we used BRAKER2 (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 MAKER 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).