https://doi.org/10.5061/dryad.rxwdbrvjq

Description of the data and file structure

We conducted a laboratory reciprocal transplant of wild F0 Amnicola limosus gastropods sourced from two habitat types (high calcium/invaded and low calcium/refuge) to measure adult survival and fecundity in home and transplant treatments of water calcium concentration (low/high) and round goby cue (present/absent). We then applied pooled whole-genome sequencing of twelve gastropod populations from across the calcium/invasion gradient. The pooled sequencing was supported by a new draft genome for A. limosus. Population genomic analyses (genome scans, analyses of genomic diversity, population structure, demographic modeling) were used downstream after processing the sequencing data with a custom bioinformatic pipeline.

Files and variables

This repository contains several top-level files and folders, as well as subfolders:

1. code_hpc: This folder contains code written to be processed on an HPC cluster.
2. code_rscripts: This folder contains code written to be processed with R on a regular computer.
3. processed_data: This folder contains data files that have been processed and are ready to use. Generally, these files are the result of running data processing scripts on the code_hpc folder on raw data in the 'raw_data' folder.
4. raw_data: This folder contains raw data files, including spreadsheets of raw data.
5. results: This folder contains statistical results, tables, and figures.

In many cases, the files are named after or contain the names of the 12 populations studied, referred to as RAF, PST, PON, PG, PDC, PB, OKA, IPE, IB, HA, GOY and BEA. In the readme, when files or folder names contain one or multilple population names, these were replaced by asterisks * for simplicity (except for the scripts for the Baypass analyses, where they refer to the replicate numbers).

Computing environment

• Bioinformatic analyses were performed using the R programming language and several softwares (fastp, Baypass, popoolation1 and popoolation2, dadi, poolFreqDiff, BUSCO, Jellyfish, samtools_flagstat). Non-HPC analyses were performed on an Asus VivoBook with an Intel Core i5 processor and 12GB of RAM, within a virtual machine (Oracle VM VirtualBox, 6GB of RAM) running on Ubuntu 20.04.5 LTS. HPC analyses were performed on the Digital Research Alliance of Canada cluster Beluga.

Other data

Raw Illumina reads are available on the NCBI SRA portal under the BioProject ID PRJNA1035459. The draft genome of Amnicola limosus has been deposited on the NCBI Genome database (BioProject PRJNA1136450).

List of folders and files

code_hpc:

1_pipeline_SNPcalling:

This folder contains code to run the bioinformatic pipeline with raw paired-end sequencing files as input and a sync file as output. The pipeline is divided into two parts

• trimming_bwa_sambam_part1: 1) trimming of raw fastq files with the function trim-fastq.pl from Popoolation1 v1.2.2 running with perl v5.30.2, 2) align trimmed reads to draft genome with bwa mem function from bwa v0.7.17, 3) filter aligned reads based on a quality of 20 with samtools view function from samtools v1.13, 4) sort bam files with samtools sort from samtools v1.13. This is done separately for each population with dedicated .sh files named pipeline_pt1_populationname.sh
• mpileup_sync_part2: pipeline_pt2_mpileup_sync.sh file to process the sorted bam files from all populations together. 1) Generate an mpileup file with samtools mpileup from samtools v1.13 using the bam_list_all.txt list of sorted bam files and the draft genome as input, 2) filter SNPs with reads having a coverage of 5 or more across all populations using grep -vw, 3) Create .sync file which is the main input files for downstream analyses and is generated with the function mpileup2sync.pl in Popoolation2 v1.10.03

2_popoolation:

This folder contains code to generate diversity indices with the Variance-sliding.pl function from Popoolation1 v1.2.2. For each population, three separate .sh files are used for each diversity index.

• *diversity.sh: ran first, 1) generate an individual pileup file from the sorted bam file using samtools mpileup from samtools v1.13, 2) calculate Tajima's pi with the sliding window method.
• *_diversity_theta.sh: calculate Watterson's Theta with the sliding window method.
• *_diversity_tajima.sh: calculate Tajima's D with the sliding window method.

3_Baypass:

This folder contains code to analyze the input genotyping data files (amnicola_baypass_input.genobaypass.sub1 to amnicola_baypass_input.genobaypass.sub27) with Baypass v2.2, with the input files being generated with poolfstat v2.1.1 (see poolfstat_amnicola.R code in code_rscripts > Poolfstat subfolder). The analysis is run in two parts, first with the core model, and then the STD model. The input genotyping data files correspond to 27 sub-dataset Baypass input files generated with the "thinning" subsampling method, with a sub-sample size of 750,000 SNPs. The Baypass program also uses the amnicola_baypass_input.poolsize (file with pool size for each population, located in processed_data > Baypass_input) as input. For a description of the options used, see the Baypass manual

• Core_model: This folder contains code files named amni_baypass_core_sub12.sh to amni_baypass_core_sub27.sh to run the Baypass core model from pairs of input genotyping data files.
• STD_model: This folder contains code files named amni_baypass_covis_sub*.sh to run the Baypass STD model from individual input genotyping data files. This analysis requires two additional input files located in processed_data > Baypass_input: 1) contrast_goby gives the round goby presence (1) or absence (-1) status population in order, which is needed to run the contrast analysis to compute the C2 statistic (for the round goby binary covariable), 2) cov_data gives the round goby presence/absence (1st line) then the water calcium concentration (2nd line). This folder contains a subfolder, Replicate_runs with code files named amni_baypass_covis_sub*_rep*.sh to run two additional replicate runs of the Baypass STD model (N= 3 runs in total).

4_dadi:

This folder contains the code used to run the demographic analyses with dadi v2.1.0 from python v3.8. The codes are based on the well-documented user guide for dadi https://dadi.readthedocs.io/[.](https://dadi.readthedocs.io/) This folder is divided into three sub-folders named fixed_bottleneck_time_*_** for each population pair, and each of these fixed_bottleneck_time__* folders contains eight sub-folders for each of the seven demographic models analyzed: 1) bottleneck and growth in both populations (most complex model) with defined effective population sizes after the split (nu1 and nu2), followed by a bottleneck in both populations followed by exponential recovery in both populations. TS is the scaled time between the split and the bottleneck and TB is the scaled time between the bottleneck and present, 2) bottleneck and growth only in the invaded population (constant Ne for the refuge population) with uneven migration, 3) only bottlenecks in both populations without recovery and uneven migration, 4) only bottleneck in the invaded population with uneven migration (constant Ne for the refuge population), 5) a simple population split at TS with uneven migration, 6) a population split with symmetric migration and 6) a population split without migration). These eight sub-folders are in the order listed below:

• twopops_bottlegrowth
• *_bottlegrowth
• twopops_splitmigbottleneck
• splitmigbottleneck*
• splitmig_unevenMig
• splitmig_evenMig
• splitnomig

Each of these sub-folders contains a python file .py used to run the optimization of the parameters for each model and a bash .sh script to run the python file on a computing cluster. The Python scripts for the optimization consists of several parts:

1. Data import and size frequency spectrum SFS (fs) generation: pop_ids and ns give the population pair ID and the number of chromosomes included (2X nb of individuals). The fs is masked for the first 5 alleles.
2. Define custom demography model: define model parameters (tailored for each model) and output the model SFS.
3. Parameter optimization: this part runs a local optimizer on the log of the parameters based on the defined grid size, the demographic model defined above, initial parameter values p0, and upper and lower bounds of the parameters. p0, the upper and lower bounds of the parameters are modified after each optimization run until the optimization reaches a plateau (log-likelihood of the three best models within 1% of each other).
4. Log-likelihood of the model calculation: outputs the log-likelihood of the model based on the model SFS with the optimized parameters and the observed SFS.
5. Plot: After each optimization run, outputs a plot of the folded joint SFS for the observed data and the model (optimized parameters) as well as the residuals of the models.

The input data for these analyses are the dadi_input_filtercov__ files located in the processed_data > dadi_input sub-folders, which were generated with the R scripts named script_prep_input_dadi__ located in the code_rscripts > dadi sub-folder.

For each population pair, an additional sub-folder named bootstrapping_uncertainties, contains two .py scripts:

• bootstrapping__**.py: code to generate 100 bootstrapped SFS with a chunk size of 100K SNPs and masking the first 5 alleles in each bootstrapped SFS.
• uncertainty_GIM_*.py: Calculate the uncertainties on the parameters for the best model (to run after finishing the optimizations for all models tested) using the Godambe Information Matrix. This script takes as input the observed SFS generated at the beginning of each of the optimization scripts and the bootstrapped SFS generated with the bootstrapping.py script. It redefines the best model in the same manner as in the optimization scripts and defines the optimized parameters as input. The scripts output the parameters standard deviations from the GIM, the estimated 95% confidence interval from the GIM, and the upper and lower bounds of the 95% confidence interval of the parameters.

5_poolFreqDiff:

This folder contains two .sh scripts to run the poolFreqDiff analysis, one for each environmental covariable analysis (independent analyses for the calcium concentration and round goby presence/absence status). This analysis takes modified .sync files as input, which were reorganized from the initial .sync file output at the end of the bioinformatic pipeline described above, by alternating one population from a low-calcium habitat (or round goby absent) and one population from a high-calcium habitat (or round goby present). The parameters for the scripts are -npops = 6 for six populations in each environmental level (low-calcium vs high-calcium or round goby absent vs round goby present) defined by -nlevels = two levels of environmental characteristics possible, -n = 40 individuals, -mincnt = 2 for a minimum read count of two per alleles, -minc = 5 reads minimum per SNP per population, -maxc = 300 reads maximum per SNP per population, -rescale = neff to rescale the allele counts to an effective sample size neff as defined in (Feder et al., 2012) and -zeroes = 1 to add one to zero count cells.

6_BUSCO:

This folder contains one .sh script to run the BUSCO analysis on the A. limosus draft genome. It uses the --offline and --mode = genome options, and takes as input the mollusca_odb10 lineage dataset which was pre-downloaded from https://busco-data.ezlab.org/v5/data/lineages/ in 2020 with the 2020-08-05 time stamp.

7_fastp:

This folder contains three .sh scripts to run fastp v.0.23.4 on the three 10X paired-end reads used for the assembly of the draft genome, to obtain the insert size, using the -Q option to disable quality filtering.

8_Jellyfish:

This folder contains one .sh script jellyfish_amnicola_illumina_3pops_R1_untrimmed.sh to estimate the genome size with Jellyfish v2.3.0. It uses as input the R1 reads of the three runs of 10X sequencing simultaneously (-F 3) using the options count 21-mers (-m 21), a hash of 2 billion elements (-s 2G), and 20 threads (-t 20). The Jellyfish output file .jf is then converted to a histogram of the k-mer occurrences in a .histo file. The resulting histogram was processed online with GenomeScope http://qb.cshl.edu/genomescope/ to obtain an estimate of the haploid genome length, the % of repeats, and the % of heterozygosity.

9_samtools_flagstat:

This folder contains one .sh script, flagstat_amnicola_sortedbam.sh to estimate the percentage of Illumina reads aligned to the draft reference genome. This was done separately for each population using the samtools flagstat function from samtools v1.18, with the _sorted.bam files as input obtained by running the bioinformatic pipeline in the 1_pipeline_SNPcalling folder.

code_rscripts

1_Poolfstat:

This folder contains three R scripts to generate input files for other downstream analyses and to plot/assess patterns of isolation by distance IBD and isolation by environment IBE:

1. poolfstat_amnicola.R: This script generates input files for other downstream analyses by converting the All.sync file output from the bioinformatic pipeline into a pooldata object, using pool sizes of 80 for all populations, a mean read count min.rc = 2, a minimum coverage per pool of min.cov.per.pool = 5, a maximum coverage per pool max.cov.per.pool = 300, a minimum minor allele frequency of min.maf = 0.0125, and removing indels. The scripts then generate input files for Baypass with the pooldata2genobaypass() function, creating 27 sub-datasets of 778,846 SNPs with the thinning method. An allele frequency file freq_amnicola.tsv is generated by dividing the reference allele (determined arbitrarily) read count by the coverage for each SNP, which is needed as input to be converted in a dadi-compatible input format. The script also generates the snpinfo_amnicola_poolfstat.tsv and coverage_amnicola_poolfstat.tsv files needed as input for other data processing scripts. Finally, the script outputs the pairwise fst matrix FSTpoolfstat_matrix.tsv file and the population-specific heterozygosities pop_heterozygosity.tsv files.
2. poolfstat_amnicola_nooutliers.R: This script similarly generates input files for other downstream analyses, but removes the outliers identified with Baypass and poolFreqDiff. It takes as input the Chr_pos file name of the scaffold + position on the scaffold as a single column separated by a dot, next to which is added the numbers 1 to 21,312,700 as an index. Then the input files outliers_baypasscore.tsv, outliers_baypassSTD.tsv, outliers_baypasscontrast.tsv, outliers_poolfreqdiff_calcium.tsv, and the outliers_poolfreqdiff_goby.tsv are imported, and the SNPs with matching scaffold + position on scaffold identified as outliers are removed. Finally, the output files are generated with the pooldata.subset() function.
3. results_IBD_IBE_plots_neutralSNPs.R: This script is used to plot pairwise FST, assess IBD and IBE with Mantel tests, and plot the regression of the genetic distance (linearized FST with the formula fst/(1-fst) as in Rousset 1997) as a function of geographical or environmental distances. It takes as input and generates several processed files located in the processed_data > Poolfstat_input folder. The figure IBD_IBE_amnicola.pdf is output as the plot of the IBD and IBE.
4. FSTpairwise_triangularmatrix: lower pairwise FST matrix
5. FSTpairwise_triangularmatrix_ordered.csv: reordered lower pairwise FST matrix by habitat type
6. FSTpairwise_fullmatrix.csv: output by the poolfstat_amnicola_nooutliers.R script, modified with the melt() function to create the IBD_fst_amnicola.tsv
7. IBD_fst_amnicola.tsv: file generated above, to which the distances calculated along rivers with Google Earth (in m) are added.
8. logdistance_amnicola.csv: full matrix of the distances calculated along rivers with Google Earth (in m) converted to a log (ln) scale.
9. linearizedfst_fullmatrix_amnicola.csv: full matrix output by the poolfstat_amnicola_nooutliers.R script, then converted to linearized FST with the formula fst/(1-fst).
10. env_inputdata_IBE_amnicola.csv: three-column file with the population name, the round goby presence/absence status, and the calcium concentration, which is then converted to a Mahalanobis distance full matrix and saved as env_Mahalanobisdist_amnicola.tsv
11. env_Mahalanobisdist_amnicola.tsv: used to do the Mantel test for IBE along with the linearizedfst_fullmatrix_amnicola.csv input file and to do a Mantel test using the distance_amnicola.csv file of distances calculated along rivers with google earth (in m) to check for a correlation between environmental and geographical distance.
12. IBE_fst_amnicola.csv: generated by adding the Mahalanobis distance to the file IBD_fst_amnicola.tsv mentioned above.

2_Baypass:

This folder contains three R scripts to process and plot the results output by Baypass. For a description of the variables in the dataset output by Baypass, see the Baypass manual

1. amnicola_results_baypassouput_core.R: This R script is used to process the output of the core model. It takes amnicore_sub*_*mat_omega.out files as input to produce the heatmap and hierarchical clustering tree plot. The file amnicore_all_summary_pi_xtx.out outputted by Baypass is used as input for the rest of the processing, which was obtained by concatenating all the summary_pi_xtx.out files from each sub-sample and removing all the headers of the sub-sample output files except the first one. It also takes the amnicola_baypass_input.snpdet.core.all file as input, which gives information on the location of the SNP on the genome (scaffold and position on scaffold), the identity of the reference allele (defined arbitrarily), and of the minor allele, and was similarly generated by concatenating all the amnicola_baypass_input.snpdet.sub files output by Baypass for each sub-sample (no header). It outputs the list of outliers for the core model defined with the recalculated p-values of the estimated XtX as p-value <= 0.001
2. amnicola_results_baypassouput_STD_contrast.R: This R script is used to process the output of the Baypass STD model for each independent run and has one part for the results of the association with calcium concentration and a second part for the contrast analysis (round goby presence/absence status).
   • Part 1: takes as input the amniSTDcontrast_rep*_summary_BFis.out output file, which was generated by combining all the *summary_betai_reg.out files from each sub-dataset into one result file, then only the columns keeping the columns COVARIABLE, MRK, and BF.db. It also takes as input the amnicola_baypass_input.snpdet.core.all file. The script calculates the median of the BFis for the three runs and outputs a list of outlier SNPs for the association with calcium when the median BFis is > 20
   • Part 2: takes as input the amniSTDcontrast_all_summary_betai_reg.out and amniSTDcontrast_all_summary_contrast.out output files, which were generated by combining all the summary_betai_reg.out or summary_contrast.out files from each sub-dataset into one result file, and removing all the headers of the sub-sample output files except the first one. It also takes as input the amnicola_baypass_input.snpdet.core.all file. The script applies the FDR procedure with the package q-value on the back-transformed p-values (initially on the -log10 scale), with an alpha level of 0.01. It outputs a list of outlier SNPs for the association with goby predation when the q-value is below alpha.
3. check_convergence_omega_STD.R: This R script is used to assess the convergence of MCMC for the STD model in Baypass by comparing the omega matrices of the three independent runs for each subsample and calculating the Förstner and Moonen distance FMD distance. It takes as input the amniSTDcontrast_sub*_mat_omega.out (rep run #1), amnicovis_sub*_rep1_mat_omega.out (rep run #2), and amnicovis_sub*_rep2_mat_omega.out (rep run #3) files as input and calculates the FMD distance for all the omega matrices, which are then collated in a table. The median FMD distance was calculated from that table.

3_Poolfreqdiff:

This folder contains one Amnicola_results_poolfreqdiff_goby_calcium.R script to process and plot the results of the poolFredDiff analyses for the association of SNPs allele frequency with the calcium concentration and the round goby presence/absence status. It takes as input the amni_poolfrediff_goby.out and amni_poolfrediff_calcium.out files. The script recalibrates p-values with a genomic inflation factor gif (Francois 2016) value of 0.85, then applies a FDR procedure on the recalibrated p-values with the q-value package. The results are saved as data tables in the results_poolfreqdiff_amnicola_goby_qval.tsv and results_poolfreqdiff_amnicola_calcium_qval.tsv files. Two lists of outliers are generated by filtering SNPs with q-value below the alpha level 0.01 and saved as outliers_poolfreqdiff.tsv and outliers_poolfreqdiff_calcium.tsv files.

4_candidates_SNPs_outliers:

This folder contains one R script, results_oultiers_candidate_snps_Baypass_poolfreqdiff.R, used to check the list of outliers output from the R scripts used to process the results of the genome scans analyses (Baypass core and STD models, poolFredDiff analyses), and to produce Venn diagram plots. It takes as input the outliers_baypasscore.tsv, outliers_baypassSTDcalcium_repruns.tsv, outliers_baypasscontrast.tsv, outliers_poolfreqdiff_calcium.tsv and outliers_poolfreqdiff_goby.tsv files.

5_dadi:

This folder contains three R scripts to generate the input data files for dadi named script_prep_input_dadi__.R. These R scripts take as input three files generated with poolfstat from the pooldata object, using the poolfstat_amnicola.R script in the **code_rscripts > 6_Poolfstat subfolder. The allele frequency file freq_amnicola.tsv reduced to two columns (only keeping the allele frequency for the population pair of interest), the snpinfo_amnicola_poolfstat.tsv and the coverage_amnicola_poolfstat.tsv files are needed as inputs. Starting from the initial SNPs dataset output by poolfstat (21,312,700 SNPs), the script includes several filtering steps: 1) SNPs were retained if they fell within the 1st and 3rd quartiles of coverage in both populations, 2) filter out SNPs that were detected as outliers (putatively under selection) in the Baypass (core and aux or STD models) and poolFreqDiff analyses and 3) removed uninformative SNPs (fixed or lost in both populations). To accommodate for large computation time during the optimization, the datasets are also thinned at random to retain a final dataset of ≈ 1 million SNPs per population pair. The dadi_input_pools function from the genomalicious R package with the “probs” parameter in the methodSFS option is then used to transform the allele frequency dataset into the SNP data format from dadi.

6_Figures:

This folder contains two R scripts to generate plots.

1. bitplot_baypass_poolfreqdiff.R: This script takes as input the output files from the amnicola_results_baypassouput_STD_contrast.R and the Amnicola_results_poolfreqdiff_goby_calcium.R scripts, which generated the results_amniSTDcalcium_repruns.tsv, results_amniSTD_contrast_goby.tsv, results_poolfreqdiff_amnicola_calcium_qval.tsv and results_poolfreqdiff_amnicola_goby_qval.tsv data tables. This scripts generate the biplots results with a) the q-value from poolFreqDiff for the calcium vs the BFis from STD analysis for calcium and b) the q-value from poolFreqDiff for the round goby predation status vs q-values from the Baypass contrast analysis for the goby round goby predation status. It outputs the biplot_baypass_poolfreqdiff.png plot.
2. snps_outliers_allele_frequency.R: This script is used to transform the column M_P (mean of the posterior distribution of the αij parameter which is closely related to the frequency of the reference allele) from the concatenated Baypass core model output amnicore_sub_summary_yij_pij.out files (removing all the headers of the sub-sample output files except the first one) into a data frame giving the allele frequency for all populations and SNPs (sub-datasets combined). Then this allele frequency table is combined with the SNP information from the amnicola_baypass_input.snpdet.core.all file. It takes as input the outlier files from the STD model in Baypass and the poolFreqDiff analyses: outliers_baypasscore.tsv, outliers_baypassSTDcalcium_repruns.tsv, outliers_baypasscontrast.tsv, outliers_poolfreqdiff_calcium.tsv and outliers_poolfreqdiff_goby.tsv. The scripts plot the allele frequency of outlier SNPs in common between poolFreqDiff and Baypass for the calcium or the round goby predation status, and uniquely associated with each covariable, as well as the allele frequency of the non-outliers (not outlier for both the poolFreqDiff and Baypass analyses and thinned). Plots are generated for each outlier SNPs depending if the allele frequency is increasing or decreasing with the covariable (water calcium concentration gradient and goby abundance).

7_Popoolation:

This folder contains one R script, results_diversity_popoolation.R, used to process the genomic diversity indices (Tajima's pi, Watterson's Theta, Tajima's D) output by Popoolation1 and of the heterozygosity output by poolfstat. It takes as input the .pi, *.theta, and the *.D files for each population the and pop_heterozygosity.csv file in the **processed_data > Diversity_indices_input* folder. The *.pi, *.theta, and the *.D data frames are then filtered to remove windows with zero coverage (V3 column) and a minimum fraction covered < 0.05. The data frames from all populations are then combined for each diversity indices, and characterized by the habitat type: goby present or goby absent. Hedge's G for different sample sizes is computed for each diversity index and habitat type to obtain the effect size of the goby predation status. Violin plots for each diversity index and colored by habitat type are generated. For the population's heterozygosity, effect sizes of the goby predation status are computed as Cohen's D and plotted by habitat type. All the plots are then combined and output as the popoolation_heterozygosity_diversity.pdf file.

8_Transplant_experiment:

This folder contains one R script, script_amnicola_exp.R, used to process the results of the reciprocal transplant experiment. It takes as input the file Amnicola_experiment.csv located in the subfolder raw_data > transplant_experiment. This script is used to remove data from combo water treatments as mortality was too high (mentioned in the paper and presented in the supplementary material, results for combo water plotted at the end of the script), calculate survival and fecundity, and prepare the data table for plotting. The script then presents the statistical analysis for the survival (GLMM) and the fecundity (GLM). Raw data and results of the models are plotted with ggplot and combined to output Figure 2.

processed_data

1_Baypass_input:

This folder contains three files that are used as input for the Baypass analysis with the core and STD model. amnicola_baypass_input.poolsize gives the pool size for each population, contrast_goby gives the round goby presence (1) or absence (-1) status population in order, and cov_data gives the round goby presence/absence (1st line) then the water calcium concentration (2nd line).

2_dadi_input

This folder contains three files that are used as input for the demographic analyses with dadi and are in the SNP data format https://dadi.readthedocs.io/en/latest/user-guide/importing-data/. The files dadi_input_filtercov__ are for each population pair analyzed, and give the allele count for each allele. These input files were generated with the script_prep_input_dadi__.R scripts located in the **code_rscripts > 5_dadi subfolder

3_Diversity_indices_input:

This folder contains the files .pi, *.theta, and the *.D for each population and the pop_heterozygosity.csv file, which were output by the *_diversity.sh, *_diversity_theta.sh and *_diversity_tajima.sh scripts in the **code_hpc > 2_popoolation **subfolder and the *poolfstat_amnicola.R script in the code_rscripts > 1_Poolfstat subfolder. They are used as input for the R script results_diversity_popoolation.R located in the code_rscripts > 7_Popoolation subfolder. The files *.pi, *.theta, and the *.D have 5 columns:

• Column 1: Scaffold
• Column 2: Position on the scaffold (sliding window of 50Kb)
• Column 3: Number of SNPs in the sliding window
• Column 4: fraction of the window covered by a sufficient number of reads (higher than min coverage and lower than max coverage)
• Column 5: diversity index

4_Poolfstat_input:

This folder contains 14 files which are used as input and are outputs of the results_IBD_IBE_plots_neutralSNPs.R script in the code_rscripts > 1_Poolfstat subfolder:

• distance_amnicola.csv: full matrix file of distances calculated along rivers with google earth (in m)
• FSTpoolfstat_matrix_neutralSNPs.ods: spreadsheet document with the pairwise FST matrix for neutrals SNPs (generated with the poolfstat_amnicola_nooutliers.R script in the code_rscripts > 1_Poolfstat subfolder) from which multiple input files found in this subfolder were generated.
• FSTpairwise_fullmatrix.csv: output by the poolfstat_amnicola_nooutliers.R script, modified with the melt() function to create the IBD_fst_amnicola.tsv.
• FSTpairwise_triangularmatrix: lower pairwise FST matrix
• FSTpairwise_triangularmatrix_ordered.csv: reordered lower pairwise FST matrix by habitat type
• linearizedfst_fullmatrix_amnicola.csv: full matrix output by the poolfstat_amnicola_nooutliers.R script, then converted to linearized FST with the formula fst/(1-fst).
• linearizedfst_fullmatrix_amnicola_nooutliers: file used for the Mantel test where the effect of removing the pairs PB-IPE and GOY-BEA considered as outliers is tested
• logdistance_amnicola.csv: full matrix of the distances calculated along rivers with google earth (in m) converted to a log (ln) scale.
• logdistance_amnicola_nooutliers: file used for the Mantel test where the effect of removing the pairs PB-IPE and GOY-BEA considered as outliers is tested
• IBD_fst_amnicola.tsv: file generated above, to which the distances calculated along rivers with Google Earth (in m) are added.
• env_inputdata_IBE_amnicola.csv: three-column file with the population name, the round goby presence/absence status, and the calcium concentration, which is then converted to a Mahalanobis distance full matrix and saved as env_Mahalanobisdist_amnicola.tsv
• env_Mahalanobisdist_amnicola.tsv: used to do the Mantel test for IBE along with the linearizedfst_fullmatrix_amnicola.csv input file and to do a Mantel test using the distance_amnicola.csv file of distances calculated along rivers with google earth (in m) to check for a correlation between environmental and geographical distance.
• IBE_fst_amnicola.csv: generated by adding the Mahalanobis distance to the file IBD_fst_amnicola.tsv mentioned above.

5_Candidates

This folder contains files of the outliers of the Baypass and poolFreqDiff generated with the amnicola_results_baypassouput_core.R, amnicola_results_baypassouput_STD_contrast.R scripts in the code_rscripts > 2_Baypass subfolder and the Amnicola_results_poolfreqdiff_goby_calcium.R script in the code_rscripts > 3_poolFreqDiff subfolder: outliers_baypasscore.tsv, outliers_baypassSTDcalcium_repruns.tsv, outliers_baypasscontrast.tsv, outliers_poolfreqdiff_calcium.tsv and outliers_poolfreqdiff_goby.tsv. These files are input for the results_oultiers_candidate_snps_Baypass_poolfreqdiff.R script in the code_rscripts > 4_candidates_SNPs_outliers subfolder.

raw_data

1_transplant_experiment:

This folder contains the raw data file Amnicola_experiment.csv for the reciprocal transplant experiment on Amnicola limosus snails sampled in the wild, which gives for each population and treatment tested the survival and number of eggs found at different time points throughout the experiment. In that data table, the following variables are recorded:

• Pop: name of the populations of origin of the snails
• origin_water: rivers where the populations were sampled. SL refers to the St. Lawrence River (high calcium and round goby present, except for PDC which is high calcium but round goby absent) and OR to the Ottawa River (low calcium and round goby absent, except for RAF which is low calcium but round goby present)
• Pot: refers to the population of origin and the water treatment. BW is for "brown water", which was sampled in the Ottawa River, GW is for "green water", which was sampled in the St. Lawrence River, CW is for "combo water" reffering to the artificial growth medium Combo without the addition of calcium, HCW is for "high combo water" reffering to the artificial medium Combo with the addition of calcium, GP is for "goby present", which correspond to the addition of goby cues to the water treatment. When GP does not appears, it means goby cues were not added.
• Traitement: treatment water, see explanation above.
• goby_cue: goby cue treatment, for GP, see above, NOGP means goby cues were not added.
• initial_nb_ind: initial number of individuals at the starting date of the experiment.
• Prise_1 to Prise_4: X note if data were recorded in a given pot at a given date.
• date1 to date4: date at which the data was recorded.
• NB_IND1 to NB_IND4: number of individuals found alive in the treatment pot at the recording date.
• NB_EGGS1 to NB_EGGS4: number of eggs found in the treatment pot at the recording date.

results

1_Baypass:

This folder contains two subfolders for the results of the Baypass core and STD models.

1. Core: This folder contains all the omega matrix files amnicore_sub*_mat_omega output by the core model, and three plots generated with the amnicola_results_baypassouput_core.R script in the code_rscripts > 2_Baypass subfolder.
2. STD_BFis_contrast: This folder contains one plot generated with the amnicola_results_baypassouput_STD_contrast.R script in the code_rscripts > 2_Baypass subfolder.

2_Candidates:

This folder contains the Venn diagram plots generated with the results_oultiers_candidate_snps_Baypass_poolfreqdiff.R script in the code_rscripts > 4_candidates_SNPs_outliers subfolder, and the allele frequency plot generated with the snps_outliers_allele_frequency.R script in the code_rscripts > 6_Figures subfolder.

3_dadi:

This folder contains the results files output by running the demographic analyses with dadi. Within the subfolder for each population pair analyzed, there is an information file info_data_input on the input dataset for dadi for each population pair analyzed, as well as the uncertainties_twopopsbottlegrowth__ file generated with the uncertainty_GIM_.py script located in the **code_hpc > 4_dadi* subfolder and subfolders within. This folder is divided into three sub-folders named fixed_bottleneck_time__**__* for each population pair, and each of these fixed_bottleneck_time_*_* folders contains eight sub-folders for each of the seven demographic models analyzed (see above).

• twopops_bottlegrowth
• *_bottlegrowth
• twopops_splitmigbottleneck
• splitmigbottleneck*
• splitmig_unevenMig
• splitmig_evenMig
• splitnomig

Each of these sub-folders contains the output files generated by running the optimization of the parameters for each model and population pair, as well as the accompanying plots of each optimization run. A .ods file summarizes the results of all the optimization runs in each model and population pair subfolders.

4_Figures_combined_Inkscape:

This folder contains seven .svg files for various figures from the main manuscript and the supplementary material which were combined or modified with Inkscape.

5_genome_processing:

This folder contains the output files obtained by running the scripts in the code_hpc > 6_BUSCO, code_hpc > 7_fastp, code_hpc > 8_Jellyfish, and the code_hpc > 9_samtools_flagstat subfolders.

• 1_BUSCO: contains the output file from the BUSCO analysis
• 2_fastp: contains the output files from running fastp on the three 10X paired-end reads used for the assembly of the draft genome, to obtain the insert size, using the -Q option to disable quality filtering.
• 3_fastqc_illumina_reads: contains two subfolders, raw and trim, which each contain the fastqc.html file output by fastqc for the raw and trimmed Illumina read for each population and R1/R2 reads.
• 4_flagstat: contains the output files obtained by running the flagstat_amnicola_sortedbam.sh script. There are 12 .csv output files for each population, and one .ods file with the combined information for all the populations
• 5_jellyfish: contains the subfolder 10X_allruns_R1, in which is located the amnicola_10X_allruns_R1.histo output file of Jellyfish, and the plot and results output from running the GenomeScope program online.

6_poolFreqDiff:

This folder contains two subfolders, calcium, and goby_predation, which each have two .pdf plot files output by running the Amnicola_results_poolfreqdiff_goby_calcium.R script in the code_rscripts > 3_Poolfreqdiff subfolder.

7_poolfstat:

This folder contains two plot .pdf files generated by running the results_IBD_IBE_plots_neutralSNPs.R script in the code_rscripts > 1_Poolfstat subfolder.

8_popoolation:

This folder contains one .pdf plot file obtained by running the R script, results_diversity_popoolation.R script in the code_rscripts > 7_Popoolation subfolder.

9_Transplant_experiment:

This folder contains the plots generated by running the script_amnicola_exp.R in the code_rscripts > 8_Transplant_experiment subfolder

Code/software

Code/software

code_hpc:

1_pipeline_SNPcalling:

Popoolation1 v1.2.2
perl v5.30.2
bwa v0.7.17
samtools v1.13

• trimming_bwa_sambam_part1: 1) trimming of raw fastq files with the function trim-fastq.pl from Popoolation1 v1.2.2 running with perl v5.30.2, 2) align trimmed reads to draft genome with bwa mem function from bwa v0.7.17, 3) filter aligned reads based on a quality of 20 with samtools view function from samtools v1.13, 4) sort bam files with samtools sort from samtools v1.13. This is done separately for each population with dedicated .sh files named pipeline_pt1_populationname.sh
• mpileup_sync_part2: pipeline_pt2_mpileup_sync.sh file to process the sorted bam files from all populations together. 1) Generate an mpileup file with samtools mpileup from samtools v1.13 using the bam_list_all.txt list of sorted bam files and the draft genome as input, 2) filter SNPs with reads having a coverage of 5 or more across all populations using grep -vw, 3) Create .sync file which is the main input files for downstream analyses and is generated with the function mpileup2sync.pl in Popoolation2 v1.10.03

2_popoolation:

Popoolation1 v1.2.2

• *_diversity.sh: ran first, 1) generate an individual pileup file from the sorted bam file using samtools mpileup from samtools v1.13, 2) calculate Tajima's pi with the sliding window method.
• *_diversity_theta.sh: calculate Watterson's Theta with the sliding window method.
• *_diversity_tajima.sh: calculate Tajima's D with the sliding window method.

3_Baypass:

Baypass v2.2

• Core_model: code files named amni_baypass_core_sub12.sh to amni_baypass_core_sub27.sh to run the Baypass core model from pairs of input genotyping data files.
• STD_model: code files named amni_baypass_covis_sub*.sh to run the Baypass STD model from individual input genotyping data files. This analysis requires two additional input files located in processed_data > Baypass_input: 1) contrast_goby gives the round goby presence (1) or absence (-1) status population in order, which is needed to run the contrast analysis to compute the C2 statistic (for the round goby binary covariable), 2) cov_data gives the round goby presence/absence (1st line) then the water calcium concentration (2nd line). The subfolder Replicate_runs has code files named amni_baypass_covis_sub*_rep*.sh to run two additional replicate runs of the Baypass STD model (N= 3 runs in total).

4_dadi:

dadi v2.1.0
python v3.8

Seven demographic models were analyzed with dedicated scripts for each population pair:

• amnicola_timefixed_twopops_splitmigbottlegrowth.py: model of bottleneck and growth in both populations (most complex model) with defined effective population sizes after the split (nu1 and nu2), followed by a bottleneck in both populations followed by exponential recovery in both populations. TS is the scaled time between the split and the bottleneck and TB is the scaled time between the bottleneck and present,
• amnicola_demo_dadi_splitmigbottlegrowth*.py: model ofbottleneck and growth only in the invaded population (constant Ne for the refuge population) with uneven migration
• amnicola_timefixed_twopops_splitmigbottleneck.py: model ofonly bottlenecks in both populations without recovery and uneven migration
• amnicola_timefixed_splitmigbottleneck*.py: model ofonly bottleneck in the invaded population with uneven migration (constant Ne for the refuge population)
• amnicola_dadi_splitmig_unevenMig.py: model of a simple population split at TS with uneven migration
• amnicola_dadi_splitmig_evenMig.py: model ofa population split with symmetric migration
• amnicola_dadi_splitnomig.py: model of population split without migration
• bootstrapping__.py: code to generate 100 bootstrapped SFS with a chunk size of 100K SNPs and masking the first 5 alleles in each bootstrapped SFS.
• uncertainty_GIM__.py: Calculate the uncertainties on the parameters for the best model (to run after finishing the optimizations for all models tested) using the Godambe Information Matrix. This script takes as input the observed SFS generated at the beginning of each of the optimization scripts and the bootstrapped SFS generated with the bootstrapping__**.py script. It redefines the best model in the same manner as in the optimization scripts and defines the optimized parameters as input. The scripts output the parameters standard deviations from the GIM, the estimated 95% confidence interval from the GIM, and the upper and lower bounds of the 95% confidence interval of the parameters.

5_poolFreqDiff

poolFreqDiff (no version)

poolfreqdiff_amnicola_calcium.sh and poolfreqdiff_amnicola_goby.sh scripts to run the poolFreqDiff analysis, one for each environmental covariable analysis (independent analyses for the calcium concentration and round goby presence/absence status). This analysis takes modified .sync files as input, which were reorganized from the initial .sync file output at the end of the bioinformatic pipeline described above, by alternating one population from a low-calcium habitat (or round goby absent) and one population from a high-calcium habitat (or round goby present). The parameters for the scripts are -npops = 6 for six populations in each environmental level (low-calcium vs high-calcium or round goby absent vs round goby present) defined by -nlevels = two levels of environmental characteristics possible, -n = 40 individuals, -mincnt = 2 for a minimum read count of two per alleles, -minc = 5 reads minimum per SNP per population, -maxc = 300 reads maximum per SNP per population, -rescale = neff to rescale the allele counts to an effective sample size neff as defined in (Feder et al., 2012) and -zeroes = 1 to add one to zero count cells.

6_BUSCO:

Busco v5.2.2

script_amnicola_busco to run the BUSCO analysis on the A. limosus draft genome. It uses the --offline and --mode = genome options, and takes as input the mollusca_odb10 lineage dataset which was pre-downloaded from https://busco-data.ezlab.org/v5/data/lineages/ in 2020 with the 2020-08-05 time stamp.

7_fastp:

Fastp v0.23.4

fastp_Amnicola_Limosa_run333_Lane2.sh, fastp_Amnicola_Limosa_run399_Lane2.sh, fastp_Amnicola_Limosa_run399_Lane3.sh to run fastp on the three 10X paired-end reads used for the assembly of the draft genome, to obtain the insert size, using the -Q option to disable quality filtering.

8_Jellyfish:

Jellyfish v2.3.0

jellyfish_amnicola_illumina_3pops_R1_untrimmed.sh to estimate the genome size with Jellyfish v2.3.0. It uses as input the R1 reads of the three runs of 10X sequencing simultaneously (-F 3) using the options count 21-mers (-m 21), a hash of 2 billion elements (-s 2G), and 20 threads (-t 20). The Jellyfish output file .jf is then converted to a histogram of the k-mer occurrences in a .histo file. The resulting histogram was processed online with GenomeScope http://qb.cshl.edu/genomescope/ to obtain an estimate of the haploid genome length, the % of repeats, and the % of heterozygosity.

9_samtools_flagstat:

Samtools v1.18

flagstat_amnicola_sortedbam.sh to estimate the percentage of Illumina reads aligned to the draft reference genome. This was done separately for each population using the samtools flagstat function from samtools v1.18, with the _sorted.bam files as input obtained by running the bioinformatic pipeline in the **1_pipeline_SNPcalling* folder.

code_rscripts

R v4.1.2 and v4.4.0

R packages:

Poolfstat v2.1.1
Data.table v1.14.2
Stringr v1.4.0
Tidyr v1.2.0
Dplyr v1.0.9
Ggplot2 v3.3.6
Cowplot v1.1.1
Reshape2 1.4.4
Ggpubr v0.4.0
Effsize v0.8.1
Ecodist v2.0.9
Vegan v2.6-2
ade4 v1.7-19
ape v5.5
ggtree v3.0.2
stats 4.4.0
RColorBrewer v1.1-3
R.utils v2.12.3
Corrplot v0.92
Gplots v3.1.3.1
Qvalue v2.26.0
VennDiagram v1.7.3
Genomalicious v0.4
Ggrepel v0.9.5
Hmisc v5.1-3
lme4 v1.1-35.4
lattice v0.22-6
MASS v7.3-6.1
Bbmle v1.0.25.1
MuMIn v1.48.4
gridExtra v2.3
faraway v1.0.8\
sjPlot v2.8.16
sjlabelled v1.2.0
sjmisc v2.8.10

Access information

Other publicly accessible locations of the data:

• The draft genome of Amnicola limosus has been deposited on the NCBI Genome database (BioProject PRJNA1136450): https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1136450
• The raw sequencing reads have been deposited in the National Center for Biotechnology Information Sequence Read Archive SRA repository (BioProject PRJNA1035459) and the accompanying metadata are also stored in the SRA (BioProject PRJNA1035459), using the Eukaryotic water MIxS package https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1035459