Code.zip

Author: Maria Dance
Date: 2024-10-01

Most code is in .r or .rmd files, which are numbered in the order in which they should be run. .rmd files are best opened in RStudio. Most scripts should run in R version 4.4.1 and RStudio version 2024.04.2. Scripts number 10 and 11 were run in a Unix environment with R version 4.5.1 and relied on the packages admixr v0.9.1 (https://cran.r-project.org/web/packages/admixr/index.html) and ADMIXTOOLS v7.0.2 (doi.org/10.1534/genetics.112.145037), which was installed via Bioconda.

There is also a folder of bash scripts that were run in a unix environment.

---

------------ R scripts ------------

1.0_plot_process_radtags_output

This script makes plots for the Stacks process_radtags per-sample statistics. This R script is best run in a Linux environment using the syntax: $ Rscript process_radtags_stats_radchecks.r.

1.1_gstacks_stats_part1.r

This script plots the per-sample statistics output from the Stacks gstacks module.

1.2_gstacks_stats_part2.rmd

This script calculates stacks-gstacks dataset-wide statistics for the number of loci, mean coverage, and phasing.

2_calc_LD_decay.rmd

This script plots the average LD across different intervals in the genome (calculated using PLINK) to identify the window size where LD decays to a background rate of r=0.1.

3_plot_admixture_output_sampling_loc.rmd

This script plots ADMIXTURE results in a bar plot using linkage disequilibrium-pruned SNPS, which were filtered using sampling sites as populations, and excluding missing populations.

4_making_popmaps_using_admixture_output.rmd

This script makes a new popmap (Stacks populations input file) for the SNP filtering to retain more SNPs in a second stage of SNP filtering (populations run). It uses ADMIXTURE ancestry proportions to assign individuals to a genetic group, which are used as populations in the new populations run, instead of the sample sites.

5.1_plot_admixture_output_subsp.rmd

This script plots new ADMIXTURE results in a bar plot using linkage disequilibrium-pruned SNPS, which were filtered using the first run of ADMIXTURE genetic groups as populations, and excluding missing populations.

5.2_plot_radpainter_output_subsp.rmd

This script plots the output from fineRADstructure output, making coancestry matrices with dendrograms of population clustering. Script adapted from Milan Malinsky (millanek@gmail.com), in turn adapted from a Finestructure R Example by Daniel Lawson.

5.3_isolation_by_distance.rmd

This script tests for isolation by distance (IBD) in the dwarf birch genetic data, implemented using the function dartR::gl.ibd(). This test compares pairwise geographic and genetic distance matrices, using permutations to test for significance. It also tests the association between the matrix of genetic distances and a model matrix of genetic cluster membership, with the matrix of geographical distances as a covariate.

6.1_making_diyabc_input.rmd

This script makes the input genetic data files for DIYABC RF using the output of the populations module of Stacks.

6.2_plotting_LDA.rmd

This script plots the LDA ordination output from DIYABC-RF to help find the best scenario and fit to the data.

6.3_plotting_diyabc_params.rmd

This script plots the posterior distributions of DIYABC RF parameters from the best model (Model set 22).

7.1_Making metadata_for_NMDS.rmd

This script makes metadata for NMDS, and it uses admixture proportions from ADMIXTURE output to assign individuals to a genetic group. This data is then matched with other metadata to make a metadata output file that is useful for e.g. NMDS. It is the same approach as in script #4.

7.2_NMMDS.rmd

This script makes NMDS plots for individuals based on their SNPs.

8_plotting_admixed_groups_ice_sheet_map.rmd

This script makes figure 1 in the manuscript and also other map figures (S24, S33).

9_genetic_summary_stats_for_paper.rmd

This script makes genetic diversity statistics for each sampling site.

FinestructureLibrary.r

This script contains functions that need to be sourced by the 5.2_plot_radpainter_output_subsp.rmd script.

---

The bash scripts were run in a Linux environment.

------------Folder "bash_script_backups"------------

All these shell files can be opened in a basic text editor. For details on the parameters for each software and Unix command, see the online manuals.

2_process_radtags.sh

This script runs the Stacks module process_radtags, which demultiplexes reads on the plates and groups them into individual-level reads. This version includes the option to run radchecks, checking that the rad cut sites are intact.

2b_merge_duplicate_indivs_radcheck.sh

This script merges the duplicated individuals from the sequencing.

3_bowtie2_align_radcheck.sh

This script aligns the reads to a reference Betula nana genome using bowtie2 and sorts the output with samtools.

4_gstacks_new_apr_23.sh

This script runs the Stacks module gstacks, excluding the four populations that were identified as hybrids with another Betula species.

4_gstacks_new_apr_23.sh

This script runs the Stacks module gstacks, excluding the populations with putative Betula hybrids.

5_populations*.sh

The scripts 5_populations_* not only run the populations module of Stacks, but they then prune the SNPs in linkage disequilibrium (LD), then run ADMIXTURE on those LD-pruned SNPs. For output in phased haplotypes (version "B"), the scripts also run the programs RADpainter and finestructure.

Guide to naming convention:

v1 = used popmaps where populations represented sampling locations
v2 - used popmaps where populations represented genetic groups as calculated from admixture proportions from an ADMIXTURE analysis conducted on the v1 output.
A = single SNP per locus output (For ADMIXTURE and DIYABC RF)
B = phased haplotypes with multiple SNPs per locus output (for fineRADstructure/RADpainter only)

"exl_hybrids" excludes the four populations that were identified as hybrids with another Betula species that were not studied.
"w_hybrids" includes those hybrid populations

There are population runs for America, Eurasia, the western Eurasian nana group, and global (not specified in file name).

There was also a population run for calculating D-statistics specifically (5_populations_glandulosa_east_outgroup).

6_process_files_for_admixtools_gland

This script was run to convert population output (vcf) to an input file that was compatible with admixr v0.9.1 (https://cran.r-project.org/web/packages/admixr/index.html) and ADMIXTOOLS v7.0.2 (doi.org/10.1534/genetics.112.145037) which was installed via Bioconda.

ldPruning_relaxed.sh

This script was used for pruning SNPS that were in Linkage Disequilibrium. For details on the parameters, see the manuscript supplementary information.

Example_of_diyabc_script.sh

This is an example of the bash scripts used to run the software DIYABC RF (Collins et al. 2021). It lists the settings used in all DIYABC runs.

data.zip

Author: Maria Dance

Date: 2024-10-01

Guide to data folders:

barcodes

Contains .txt and .xlsx files with the unique DNA read barcodes and the corresponding individual ID. The files with _no_cut_site suffix contain the barcodes without the sfbl restriction enzyme cut site. These files are used in the pipeline to demultiplex the raw sequence reads from plate-level to individual-level reads.

bowtie2_logs

Contains .txt and .xlsx files with logs (summaries) taken from the output of the alignment software Bowtie2. Includes, on a per-individual basis: the number of reads pre-alignment, the number of reads post-alignment, and the number of reads lost in the alignment process. The full_ prefix refers to the dataset containing duplicates from the same individual, the merged_ prefix refers to the dataset once all the duplicate reads have been merged for each individual.

calc_LD_decay

Contains files related to the identification of sliding windows for pruning SNPs in linkage disequilibrium (LD).

Average LD was calculated across intervals in the genome in 100 kb bins using PLINK v1.9 and haplotype variant call format files (VCFs). The LD is represented by an R2 value, and the files in this folder contain R2 values for each locus, which, when visually inspected with a graph, identify a window size of 33 kb where LD decayed to a background rate of 0.1. VCF files were then filtered in windows of 33 kb, sliding by 10 kb, with a threshold of correlation coefficient (r2) > 0.2.

d_stats

Contains some text files that contain the results of D statistics and F4 ratios to test the ancestral origins of the Tundrarum population, and to formally test for admixture. These data were created using the scripts 10_statistics_gland_east_outgroup.r and 11_f4_statistics_gland_east_outgroup.r. The scripts were run in a Unix environment and relied on the packages admixr v0.9.1 (https://cran.r-project.org/web/packages/admixr/index.html) and ADMIXTOOLS v7.0.2 (doi.org/10.1534/genetics.112.145037), which was installed via Bioconda.

diyabc

Contains the input SNP files for DIYABC-RF and the output. Note that only the output for model set 22, the best-supported model, is provided due to large reftable.bin file sizes. The folder with the suffix _ratios contains the output for parameter estimation using ratios of time divided by effective population size. We estimated both types of parameters and found local posterior RMSE estimates to be lower in raw parameters.

The modelset22 folder contains model choice results (cmdline_modelchoice), parameter estimates using 100,000 simulated SNP datasets (param_est_100k), and files created using the DIYABC-RF GUI, which are useful for generating images of the various scenarios.

gstacks_logs

Contains .txt log files (summaries) on gstacks output (gstacks.txt), including per-individual: number of loci, number of reads used to genotype to loci, and mean effective coverage (efective_coverages_per_sample). The per-individual number and fraction of aligned reads and mapping quality (frac_aligned_reads_per_sample). A log file summarising the results of haplotype phasing (gstacks_phasing_results). Also includes a PDF of plots showing coverage, and number/fraction of kept reads on an individual and population basis. Here, log files are only included for the most recent gstacks run, which was used in the paper (new_april_23).

ice_sheet_data

The ice sheet data were used for plotting and comparisons with DIYABC-RF time parameter estimates. Note that due to third-party restrictions, this data is not stored here, but details on how to access the data are in the main manuscript. The shapefile with the world outline used in plotting was also stored here (data from Natural Earth, Free vector and raster map data @ naturalearthdata.com.).

metadata

Contains population metadata, including geographic coordinates (2022_final_betula_pop_coordinates.csv) and collection metadata (Betula_nana_metadata_MED.csv), and other metadata generated during the course of the study.

popmaps

Contains the .txt popmap files (individual id and corresponding population id) used as input for the populations module of stacks for filtering and generating output file formats. "pop_map_merged_w_humilis_gstacks" includes the Betula humilis population. "pop_map_exl_hybrids_" excludes the Betula humilis population, and the four populations which were identified as hybrid with another Betula species were was not studied ("pop_map_w_hybrids_sites" includes those hybrid populations). There are popmaps where the original sampling site is the population "sites", and popmaps where the population was defined using the ADMIXTURE proportions, or genetic groups ("admixture"). There are popmaps for America, Eurasia, the western Eurasian Nana group, and global (not specified in file name). Some popmaps were used for a population run specifically to calculate D-statistics, where Glandulosa East was specified as an outgroup (sample_pops_gland_east and pop_map_gland_east), and the Atlantic nana populations were specified as separate (atlantic_nana_group).

populations_output

Guide to naming convention:

v1 = used popmaps where populations represented sampling locations

v2 - used popmaps where populations represented genetic groups as calculated from admixture proportions from an ADMIXTURE analysis conducted on the v1 output.

A = single SNP per locus output (For ADMIXTURE and DIYABC RF)

B = phased haplotypes with multiple SNPs per locus output (for fineRADstructure/RADpainter only)

There is also population output from the run, which was used to calculate D statistics: gland_east_outgroup.

radchecks

This folder contains the log files, summaries, and plots generated from the process_radtags module of Stacks.

Process_radtags checks that the barcode and RAD cutsite are intact, then demultiplexes the data (grouping reads from plates to certain individuals). Includes the top lines (head) of the reads on each plate (e.g. plate1_head), the log file of the process_radtags process (plate1_process_radtags.radseq_download), a list of the unrecognised RAD sites (plate1_rad_sites_unrecog), and plots summarising the proportion of the library each individual's reads are, the number and proprtion of retained reads per individual and per sampling site. There are plots for the individuals, including duplicates (process_radtags_stats_plate1.pdf) and plots where the duplicates have been merged (process_radtags_stats_plate1_merged.pdf).

rangemaps

This folder contains shapefiles used to create Figure 1 in the manuscript. The shapefiles are species ranges which were digitised in QGIS from range maps in The Flora of North America (Furlow 2004) and the Atlas of North European vascular plants north of the Tropic of Cancer (Hultén & Fries, 1986). ------------------------------------------------------------------------------------------------------------------------------------------------------

There are also several files not in a folder, as they were created during the processing pipeline and accessed multiple times throughout the process:

• 2024_09_16_genetic_sum_stats.csv Genetic summary stats for each sampling site ( calculated in the script 9_genetic_summary_stats_for_paper.rmd).