This repository contains scripts and workflows for population genomic analysis of four invertebrate species from the San Solomon Springs Complex in West Texas. The analysis pipeline includes data preprocessing, assembly of reference sequences (pseudo-reference genomes), read mapping, variant calling, and population genetic analyses using Genotyping-by-Sequencing (GBS) data.

Geographic Context:  The San Solomon Spring Complex lies between the Davis Mountains and the Toyah Basin in the arid Trans-Pecos region of Texas. Springs are primarily fed by water from the Apache Mountains, with some contribution from the Davis Mountains following large rainfall events. The system includes Phantom Lake Spring (6 km west of Balmorhea State Park), San Solomon Springs (at Balmorhea State Park), and the Sandia springs (8 km east of Balmorhea State Park, Nature Conservancy Preserve).

Table of Contents

• Methodology Notes
  • GBS Library Preparation
• Overview
• Scripts Directory Structure
• Prerequisites
  • Software Packages Used
• Labels for Individuals
  • Barcode Key Format
  • Individual ID Format
• Workflow Overview
  • 1. Data Preprocessing
    • Barcode Parsing
    • Create a List of All Individuals
    • Split FASTQ by Individual
    • Quarantine Low-Quality Samples
  • 2. Assembly of Reference Sequences
    • Step 1: Unique Sequence Extraction
    • Step 2: Reference Assembly
    • Remove Contamination
  • 3. Read Mapping
    • BWA Alignment
    • SAM to BAM Conversion
  • 4. Filtering Loci
    • Filter VCF - First Round
    • Collect Depth Statistics
    • Filter VCF - Second Round
  • 5. Genotype Likelihoods
    • VCF to MPGL Conversion
    • Randomly select one variant site per scaffold
  • 6. Running Entropy
    • Prepare Entropy Input Files
    • Create ldak*.txt Files
    • Prepare Last File For Entropy
    • Run Entropy Analysis
  • 7. Entropy Results Analysis
    • Extract Entropy Results
    • Assess MCMC Convergence
    • Create Ancestry Bar Plots
    • Principal Component Analysis
    • Calculate Population Differentiation
    • Genetic Diversity Estimation
      • Run ANGSD Theta Estimation
      • Process and Visualize ANGSD Output
    • Isolation by Distance Analysis
      • Data File
      • Run Analysis
• Acknowledgments

Methodology Notes

GBS Library Preparation

• Whole genomic DNA was extracted using Qiagen DNAeasy Extraction Kits
• Genomic DNA fragmented with restriction enzymes EcoRI and MseI
• Meyer-Kircher adaptor sequences (Meyer and Kircher 2010) ligated with individual-specific barcodes
• PCR was conducted in a 20µL volume using the Illumina PCR primers with the following thermocycler profile: 90°C for 30s; 30 cycles of 98°C for 20s, 60°C for 30s, 72°C for 40s; 72°C for 10 min
• Size selection: 250-350 bp fragments via gel electrophoresis
• Sequencing: Illumina HiSeq 2500 (SR100) at University of Texas Genomic Sequencing and Analysis Facility

Overview

This pipeline processes raw sequencing data from Genotyping-by-Sequencing (GBS) and performs population genomic analyses including:

• Assembly of reference sequences (pseudo-reference genomes)
• Variant calling and filtering
• Population structure analysis using Entropy (hierarchical Bayesian model)
• Principal Component Analysis (PCA)
• Population differentiation (Nei's GST)
• Genetic diversity estimation (Watterson's theta, nucleotide diversity)
• Visualization of results

The study examines population structure and genetic diversity across seven sites in the San Solomon Spring system, encompassing the entire known ranges of these groundwater-dependent invertebrates.

Scripts Directory Structure

scripts
├── 01_data_preprocessing/          # Barcode parsing and quality control
│   ├── GBS_Balmorhea_BarcodeID.csv
│   ├── IDfileMaker2.pl
│   ├── parse_barcodes.pl
│   ├── quarantine_script.sh
│   └── splitFastqMF.pl
├── 02_assembly/                   # Assembly of reference sequences (pseudo-reference genomes)
│   ├── ddocent.sh
│   ├── ddocentPart2.sh
│   └── tapewormRemover.pl
├── 03_read_mapping/              # Read alignment to reference
│   ├── runbwa_aln_samse_forks.pl
│   └── sam2bam.pl
├── 04_filtering/              # VCF filtering and quality control
│   ├── depthCollector.pl
│   ├── vcfFilter_more.pl
│   └── vcfFilter.pl
├── 05_genotype_likelihoods/     # Prepare and format likelihoods
│   ├── gl2genest.pl
│   ├── locusLister.pl
│   ├── sort_mpgl_loci.pl
│   └── vcf2mpgl.pl
├── 06_running_entropy/        # Preparing to run Entropy
│   ├── add_pop_info.pl
│   ├── entropy_IDfileMaker2.pl
│   └── Entropy_infiles.R
├── 07_entropy_results_analysis/    # Population genetic analyses
│   ├── all_species_mantel_test.R
│   ├── All_ThreeSpecies_mantel_test_data.csv
│   ├── angsd_output.R
│   ├── getThetasAngsd.sh
│   ├── GR_and_ESS.R
│   ├── GstCode_entropy.R
│   ├── K_mean_barplot.R
│   ├── pca.R
│   └── run_estpost.sh
└── indexed_phiX_genome/        # Reference genomes and indices
    ├── e_phiX.1.ebwt
    ├── e_phiX.2.ebwt
    ├── e_phiX.3.ebwt
    ├── e_phiX.4.ebwt
    ├── e_phiX.rev.1.ebwt
    └── e_phiX.rev.2.ebwt

Prerequisites

Software Packages Used

• Perl v5.30.0 with modules:
  • Parallel::ForkManager
• Bash v5.0.17
• R v4.2.1 with packages:
  • ggplot2
  • MASS
  • coda
  • reshape2
  • RColorBrewer
  • dplyr
• RStudio with R v4.5.2
  • vegan
  • tidyr
  • dplyr
• Relevant Tools:
  • mawk (v1.3.4) - interpreter for the AWK Programming Language
  • BWA (v0.7.17-r1188) - Read alignment (Burrows-Wheeler Aligner)
  • samtools (v1.9) - SAM/BAM manipulation
  • bcftools (v1.9) - Variant calling
  • cd-hit-est (v4.8.1) - Sequence clustering
  • Bowtie (v1.2.3) - PhiX control sequence removal
  • Entropy (v1.2) - Population structure analysis (Gompert et al. 2014)
  • ANGSD (v1.16) - Genetic diversity estimation

Labels for Individuals

Barcode Key Format

The barcode key file (01_data_preprocessing/GBS_Balmorhea_BarcodeID.csv) should have the format:

barcode_id,barcode,id
1,TGGCTAGAGCCAATTC,GCN_001
2,ACTGGTCGCCAATTC,GCN_002
...

Where:

• barcode_id: Unique barcode identifier
• barcode: Barcode sequence (includes restriction site)
• id: Individual identifier (format: SPECIES_LOCATION_NUMBER)

Individual ID Format

Individual IDs follow the pattern: XXX_### where:

• First 2-3 letters: Species and/or location code
  • G: Gammarus (GCN, GES, GHC, GPH, GSC, GSP, GWS for different locations)
  • P: Pyrgulopsis (PCN, PES, PPH, PSP for different locations)
  • T: Tryonia (TES, TPH, TSP for different locations)
• Number: Individual number (e.g., 001, 002, etc.)

Population codes used in analyses:

• PH = Phantom Lake Spring
• SP = San Solomon Pool
• CN = San Solomon Canal
• HC = Hubbs Ciénega
• SC = San Solomon Ciénega
• WS = West Sandia Spring
• ES = East Sandia Spring

Note: To process Gammarus without Gammarus hyalelloides, remove the individuals from Phantom Lake Spring, or, in other words, the individuals starting with code GPH.

Final datasets:

• Gammarus hyalelloides + G. balmorhea : 136 individuals, 1,576 SNPs (7 sites)
• Gammarus balmorhea only: 127 individuals, 1,589 SNPs (6 sites)
• Pyrgulopsis texana: 57 individuals, 691 SNPs (4 sites)
• Tryonia cheatumi: 26 individuals, 753 SNPs (3 sites)

Workflow Overview

Note: Many scripts use input from the same directory that the script is in, or may need a relative path to the input file in the command. Be mindful of the directories in which you have your files and scripts.

1. Data Preprocessing (01_data_preprocessing/)

Begin by copying the 6 e_phiX files in indexed_phiX_genome/ to your directory containing the raw read data file.

The following command was ran in our bash terminal, using bowtie to assemble the raw reads to the phiX genome.

bowtie -v 3 -a --al aligned_phiX.txt --un unaligned_NOTphiX.txt --suppress 1,2,3,4,5,6,7,8 e_phiX GBS_rawreads.fastq.gz hitFile.txt &

After finishing, rename the output file.

mv unaligned_NOTphiX.txt raw_reads.fastq

Barcode Parsing

perl parse_barcodes.pl GBS_Balmorhea_BarcodeID.csv raw_reads.fastq

This script:

• Removes barcode sequences from reads
• Removes adapter sequences from 3' end
• Sequences not matching any barcodes are written to a separate file

Input:

• GBS_Balmorhea_BarcodeID.csv: Barcode key file with format barcode_id,barcode,id
• Raw FASTQ file(s) with phiX

Output:

• parsed_*.fastq: Parsed FASTQ files
• miderrors_*.fastq: Reads with uncorrectable barcode errors

Create a List of All Individuals

Generate ID files from barcode key:

perl IDfileMaker.pl GBS_Balmorhea_BarcodeID.csv

Input:

• GBS_Balmorhea_BarcodeID.csv

Output:

• IDfile_GBS_Balmorhea_BarcodeID.csv: Sorted individual IDs

Split FASTQ by Individual

Split FASTQ files into individual files:

perl splitFastqMF.pl IDfile_GBS_Balmorhea_BarcodeID.csv parsed_raw_reads.fastq

Input:

• IDfile_GBS_Balmorhea_BarcodeID.csv: Sorted individual IDs, one per line
• Parsed FASTQ file(s)

Output:

• Individual FASTQ files named by ID

Individual FASTQ files can now be separated into different folders for each species code. The following workflow applies to each species, but the included scripts and directions use the species codes for Tryonia. Where appropriate, change the species prefix and location codes from those for Tryonia to those for Pyrgulopsis or Gammarus, as described in the Labels for Individuals section above.

Quarantine Low-Quality Samples

Move the script for quarantining individuals into each respective species folder. In each species folder, create another folder named quarantined_individuals_folder:

mkdir quarantined_individuals_folder

Alternatively, create the output path you'd prefer and modify the script's variable quarantine_dir accordingly.

Identify and quarantine samples with insufficient sequencing depth. We quarantined individuals that had less than 30,000 sequences:

bash quarantine_script.sh

Parameters:

• Minimum sequences: 30,000 (adjustable in script, on line 17)
• Species prefix: T* (adjustable in script, on line 12 where for loop is declared)

Input:

• The script looks at all FASTQ files starting with the specified species letter in the current directory

Output:

• quarantined_individuals_folder/: Quarantined low-quality samples are moved to this folder
• quarantined_individuals.txt: List of quarantined samples

2. Assembly of Reference Sequences (02_assembly/)

Here we do not use a reference for reference-based assembly, but collect homologous sequence reads together to find variable sites. It is adapted from Jon Puritz' dDocent wrapper (http://ddocent.com).

Make a list of fastq files:

ls *.fastq > namelist

Remove '.fastq' from namelist:

sed -i'' -e 's/.fastq//g' namelist

Note: '' after -i is two single apostrophes, not a singular quotation mark

Compress fastq files. One way is to use gzip:

gzip *.fastq

Step 1: Unique Sequence Extraction

Extract unique sequences and assess coverage:

bash ddocent.sh

Requirements:

• namelist: File containing one sample name per line (without .fastq extension)
• All FASTQ files in current directory

Output:

• *.uniq.seqs: Unique sequences per individual
• uniq.seqs: Combined unique sequences
• uniqCperindv: Unique sequences with coverage per individual
• uniqseq.peri.data: Coverage statistics
• zippedFastqs/: Directory with all the zipped fastq files moved into it

Step 2: Reference Assembly

Assemble reference sequences using CD-HIT clustering:

bash ddocentPart2.sh

Parameters:

• Minimum coverage: 4 reads per individual (sequences must have at least 4 copies shared by at least 4 individuals)
• CD-HIT identity thresholds: 90% homology used in final analysis (adjustable in line 9 of script after -c flag where there is value of 0.9)

Input

• uniqCperindv

Output:

• reference.fasta: Assembled reference sequences (pseudo-reference genome)

Remove Contamination

Remove tapeworm contamination from reference:

perl tapewormRemover.pl reference.fasta

Output:

• good_reference.fasta: Cleaned reference sequences

3. Read Mapping (03_read_mapping/)

BWA Alignment

Align reads to scaffolds from previous step. You will need the IDfile_GBS_Balmorhea_BarcodeID.csv file.

Run the following command using Burrows Wheeler Aligner (BWA) to index the reference .fasta file. Change projectName in projectName_ref to your preference.

bwa index -p projectName_ref -a is good_reference.fasta

Output:

• projectName_ref.sa
• projectName_ref.amb
• projectName_ref.ann
• projectName_ref.pac
• projectName_ref.bwt

Copy these files and the runbwa_aln_samse_forks.pl script into your zippedFastqs/ directory.

Edit the runbwa_aln_samse_forks.pl script to replace the letters where the script says YOURNAME and YOURNAME_ref on lines 19 and 20 to be your projectName and projectName_ref, respectively.

Run the script with the following command:

perl runbwa_aln_samse_forks.pl *.fastq.gz

Parameters (see in script):

• Threads: 24 (adjust to available CPUs)
• BWA parameters: -n 4 -l 20 -k 2 -t 20 -R 20 -q 10
• Reference: YOURNAME_ref (update to your reference name)

Output:

• aln*.sai: BWA alignment index files
• aln*.sam: SAM alignment files

SAM to BAM Conversion

Convert, sort, and index SAM files:

perl sam2bam.pl *.sam

Parameters:

• Threads: 24 (adjustable in script)

Output:

• *.bam: BAM files for each individual
• *.sorted.bam: Sorted BAM files
• *.sorted.bam.bai: BAM index files

4. Filtering Loci (04_filtering/)

You will need good_reference.fasta in your directory with the *.sorted.bam files to run the command:

bcftools mpileup -d 8000 -o out1.bcf -O b -I -f good_reference.fasta aln*sorted.bam

Next, run the following:

bcftools call -c -V indels -v -p 0.05 -P 0.001 -o variants.vcf out1.bcf

Filter VCF - First Round

Apply initial quality filters:

perl vcfFilter.pl variants.vcf

Filtering Criteria:

• Minimum coverage: Set to twice the number of individuals sampled for each species. For example, there were 26 Tryonia individuals used in our final analysis, so $minCoverage variable was set to 52 on line 13 of the script. For both species of Gammarus use a value of 272, for Gammarus balmorhea only use 254, and for Pyrgulopsis use 114.
• Maximum missing: 50% of individuals, you can adjust this on line 21 of the script where the $miss variable is declared. Since there are 26 Tryonia individuals used, 50% of that means that $miss is set to 13. For both species of Gammarus use a value of 68, for Gammarus balmorhea only use 63 (127 / 2 was rounded down), and for Pyrgulopsis use 28 (57 / 2 was rounded down).
• Minimum mapping quality (MQ): 30
• Allele frequency range: 0.05 - 0.95
• Maximum rank sum test statistics: 1.96 (BQB, MQB, RPB)
• Remove fixed alternative alleles
• Remove INDELs
• Remove multi-allelic sites

Input:

• variants.vcf file created by bcftools

Output:

• filtered_firstRound_variants.vcf

Collect Depth Statistics

Inspect coverage depth for remaining loci. We aim to filter loci which have extremely high sequencing depth in case they are assemblies of paralogous loci:

perl depthCollector.pl filtered_firstRound_variants.vcf

Output:

• depth_filtered_firstRound_variants.vcf: Depth statistics per variant

Read this file into R. Determine the mean + 2 * standard deviation:

dat<-read.table("depth_filtered_firstRound_variants.vcf",header=F)
sd(dat$V3)
mean(dat$V3) + 2 * sd(dat$V3)

Filter VCF - Second Round

Apply additional depth-based filtering. This script is similar to vcfFilter.pl, but you will need to set the added $maxCoverage variable to be the mean(dat$V3) + 2 * sd(dat$V3) value calculated in R.

perl vcfFilter_more.pl filtered_firstRound_variants.vcf

Additional Criteria:

• Maximum coverage: Based on mean + 2SD from depth statistics

Output:

• filtered_secondRound_filtered_firstRound_variants.vcf

Rename the file:

mv filtered_secondRound_filtered_firstRound_variants.vcf doubleFiltered_variants.vcf

5. Genotype Likelihoods (05_genotype_likelihoods/)

VCF to MPGL Conversion

Convert VCF to Multiple Population Genotype Likelihood (MPGL) format:

perl vcf2mpgl.pl doubleFiltered_variants.vcf

Output:

• doubleFiltered_variants.mpgl: MPGL format file with genotype likelihoods
• indiv_ids.txt: Individual ID mapping file needed later to prepare to run entropy software

Randomly select one variant site per scaffold

perl locusLister.pl doubleFiltered_variants.vcf

Output:

• doubleFiltered_variants.list: a list of all loci with scaffold number with three columns, scaffold number:position, and position

Read the doubleFiltered_variants.list into R and randomly pick one SNP per contig

dat<-read.table("doubleFiltered_variants.list",header=F)
contigs<-unique(dat$V1)
loci<-matrix(NA,nrow=length(contigs),ncol=1)
for (i in 1:length(contigs)){
  sites<-dat[dat$V1 == contigs[i],]
  rand<- sample(1:length(sites$V1),1)
  tmp<-as.character(sites$V2[rand])
  loci[i,1]<-tmp
  }
write.table(loci,file="loci.txt",quote=FALSE,row.names=F,col.names=F)

Get the genotype likelihoods for the specific loci in the list.

perl sort_mpgl_loci.pl loci.txt doubleFiltered_variants.mpgl

Input:

• loci.txt: List of loci to extract

Output:

• loci.txtdoubleFiltered_variants.mpgl: Specific loci genotype likelihoods

Transform the genotype likelihood file into a genotype matrix with the mean genotype (point estimate) for each individual for each variant site.

perl gl2genest.pl loci.txtdoubleFiltered_variants.mpgl

Output:

• gl_loci.txtdoubleFiltered_variants.mpgl

This converts the phred score-scaled likelihoods into mean genotypes. The result is a number between 0 and 2, with 0 being homozygous for the reference allele, 1 (heterozygous), or 2 (homozygous for alternative allele). The values are continuous scores, meaning they reflect if there is uncertainty about whether an individual is homozygous or heterozygous.

6. Running Entropy (06_running_entropy/)

Prepare Entropy Input Files

Create input files for Entropy analysis:

perl entropy_IDfileMaker2.pl indiv_ids.txt

Output:

• entropy_indiv_ids.txt: Entropy ID file
• header_ids.txt: Header file

Concatenate header_ids.txt with the genotype likelihood matrix file

cat header_ids.txt gl_loci.txtdoubleFiltered_variants.mpgl > f_gl_loci.txtdoubleFiltered_variants.mpgl

Use R to transpose the genotype matrix

dat<-read.table("f_gl_loci.txtdoubleFiltered_variants.mpgl",header=F)
tdat<-t(dat)
write.table(tdat,file="tran_loci.txt",row.names=F,col.names=F,quote=F)

Add a column with population names to the transposed loci

perl add_pop_info.pl tran_loci.txt

Create ldak*.txt Files

Modify and run the R script Entropy_infiles.R to make the ldak.txt files needed to run entropy. The script needs to be changed for each species to reflect the number of individuals and the number of loci.

The expected values used in the first three lines in the script for number of individuals and loci (when analysis is done with 90% homology and 50% miss parameters as set in 04_filtering) are below:

The first three lines in Entropy_infiles.R for Pyrgulopsis should be:

gprob<-matrix(scan("names_tran_loci.txt", n = 57 * 693, sep = " ", what = 'character'), nrow=57,ncol=693,byrow=T)
gprobMat<-gprob[,-c(1:2)]
gprobMat2<-matrix(as.numeric(gprobMat),nrow=57,ncol=691,byrow=F)

The first three lines in Entropy_infiles.R for Tryonia should be:

gprob<-matrix(scan("names_tran_loci.txt", n = 26 * 755, sep = " ", what = 'character'), nrow=26,ncol=755,byrow=T)
gprobMat<-gprob[,-c(1:2)]
gprobMat2<-matrix(as.numeric(gprobMat),nrow=26,ncol=753,byrow=F)

The first three lines in Entropy_infiles.R for Gammarus balmorhea only should be:

gprob<-matrix(scan("names_tran_loci.txt", n = 127 * 1591, sep = " ", what = 'character'), nrow=127,ncol=1591,byrow=T)
gprobMat<-gprob[,-c(1:2)]
gprobMat2<-matrix(as.numeric(gprobMat),nrow=127,ncol=1589,byrow=F)

The first three lines in Entropy_infiles.R for *Gammarus hyalelloides + G. balmorhea should be:

gprob<-matrix(scan("names_tran_loci.txt", n = 136 * 1578, sep = " ", what = 'character'), nrow=136,ncol=1578,byrow=T)
gprobMat<-gprob[,-c(1:2)]
gprobMat2<-matrix(as.numeric(gprobMat),nrow=136,ncol=1576,byrow=F)

Output:

• This R script produces 8 ldak*.txt files, where the star will be replaced by a single digit number from 2 to 8.
  Note: Analysis was completed using one greater than the number of subpopulations for each species. Both species of Gammarus will use all 8 files. Gammarus balmorhea only will use ldak2.txt through ldak7.txt. For Pyrgulopsis you will only use ldak files 2-5, and for Tryonia you will only use ldak files 2-4.

To run the script in bash, use:

Rscript Entropy_infiles.R

Prepare Last File For Entropy

Add IDs to common loci genotype file.

Make a header file by adding a line to a copy of previously made file.

cp entropy_indiv_ids.txt headerT.txt

At the very top of headerT.txt, add a new line that has the number of individuals, number of loci, and 1.
For Tryonia this should look like: 26 753 1
For Pyrgulopsis: 57 691 1
For Gammarus balmorhea only: 127 1589 1
For all of Gammarus: 136 1576 1

Run Entropy Analysis

Run Entropy for population structure analysis:

See an example of a command below for Tryonia with 2 clusters and the first (c0) of two chains that will be run. The command will need to be run a second time, after renaming the output file, for the second chain. Then the -k value and ldak*.txt input will need to be changed to specify a different number of clusters before running twice again.

entropy -i entropy.mpgl -o out_Tryonia_k2c0.hdf5 -l 105000 -b 5000 -t 10 -s 20 -e .01 -k 2 -q ldak2.txt -m 1 -w 0 &

Parameters of Note:

• -l: Number of MCMC steps for the analysis (varies by species: Gammarus 505,000; Pyrgulopsis 205,000; Tryonia 105,000)
• -b: Burn-in steps (5,000 iterations)
• -t: Thin MCMC samples by recording every nth value (10 was used for Tryonia and Gammarus, 20 was used for Pyrgulopsis)
• -k: Number of clusters (tested K = 2 to number of subpopulations + 1)
• -q: Specify which ldak file to use. The ldak file corresponding the same value of k should be used for both chains created.
• Other flags or options can be read about in the Entropy documentation, and were kept the same as shown in the example command above.

Output:

• out_*.hdf5: Entropy output files

7. Entropy Results Analysis (07_entropy_results_analysis/)

Extract Entropy Results

Extract posterior probabilities and statistics:

bash run_estpost.sh

Note: you will need to edit the script's DATA and DATA2 variables to provide the path to both relevant chains and edit KVAL for the number of population clusters. This script will need to be run once for each value of k being used for each species.

Output:

• q_k*c*.txt: Population assignment probabilities
• q_k*c*_mcmc.txt: MCMC samples
• gprobk*c*.txt: Genotype probabilities
• devk*c*.txt: Deviance statistics
• DICk*c*.txt: DIC values

Assess MCMC Convergence

Check Gelman-Rubin statistics and Effective Sample Size:

Rscript GR_and_ESS.R

Note: the relevant input files will need to be edited and specified, q_k*c0_mcmc.txt on line 6 of the script and q_k*c1_mcmc.txt on line 11 of the script, where the star is replaced by the k value being analyzed. This script will need to be run once for each value of k being used for each species. Be sure to have the relevant species indiv_ids.txt file in the same directory before running as well.

Convergence criteria:

• Gelman-Rubin (GR) scores: ideally the median ≤ 1.2 across individuals
• Effective Sample Size (ESS): at least a few hundred with a higher value being more optimal

Output:

• Summary statistics for convergence diagnostics

Create Ancestry Bar Plots

Individual ancestry estimates are normalized to correct for small deviations from chain non-convergence or numerical rounding. Generate barplots of population assignments:

Rscript K_mean_barplot.R

Parameters (adjustable in script):

• k: Number of populations (tested K = 2 to number of subpopulations + 1). Change this on line 2 of the script. The script will need to be run once for each value of k for each species.
• Site ordering: Varies by species, change one line 32 of the script. For Tryonia the ordering should be ("TPH", "TSP", "TES"), for Pyrgulopsis use (PPH, PSP, PCN, PES), and for Gammarus balmorhea only use ("GSP", "GCN", "GHC", "GSC", "GWS", "GES") or ("GPH", "GSP", "GCN", "GHC", "GSC", "GWS", "GES") for Gammarus hyalelloides + G. balmorhea.
• Output file name can be changed on line 56.

Output:

• YOUR_OUTPUT_NAME_HERE_k*_means_barplot.pdf: Ancestry barplot showing normalized admixture proportions

Principal Component Analysis

Perform PCA on genotype probabilities:

Rscript pca.R

Parameters (adjust in script):

• ks: specify range of k values to be included, line 10 in the provided script
• chains: keep this at 0:1 if you used 0 and 1 to number your MCMC chains output by Entropy, line 11
• n_ind: number of individuals in the species/analysis, line 12
• n_loci: number of loci being used, line 13
• Codes for mapping the locations individuals were sampled from can be changed on line 56 of the script
• Shapes for individuals from different locations can be specified starting at line 73 in the script
• title of the bar plot: can be changed on lines 93 and 113 of the script
• output file name: can be changed on lines 127 and 128 of the script

For the location codes starting at line 56 we used the following:
-Gammarus balmorhea only

  location_names <- c(
  GSP = "BSP San Solomon Pool",
  GCN = "BSP Canal",
  GHC = "BSP Hubbs Cienega",
  GSC = "BSP San Solomon Cienega",
  GWS = "NC West Sandia",
  GES = "NC East Sandia"
)

-Gammarus hyalelloides + G. balmorhea

location_names <- c(
  GPH = "Phantom Springs",
  GSP = "BSP San Solomon Pool",
  GCN = "BSP Canal",
  GHC = "BSP Hubbs Cienega",
  GSC = "BSP San Solomon Cienega",
  GWS = "NC West Sandia",
  GES = "NC East Sandia"
)

• Pyrgulopsis

location_names <- c(
  PPH = "Phantom Springs",
  PSP = "BSP San Solomon Pool",
  PCN = "BSP Canal",
  PES = "NC East Sandia"
)

• Tryonia

location_names <- c(
  TPH = "Phantom Springs",
  TSP = "BSP San Solomon Pool",
  TES = "NC East Sandia"
)

For the location shapes we used the following:
-Gammarus balmorhea only

location_shapes <- c(
  "BSP San Solomon Pool" = 15,
  "BSP Canal" = 15,
  "BSP Hubbs Cienega" = 15,
  "BSP San Solomon Cienega" = 15,
  "NC West Sandia" = 17, 
  "NC East Sandia" = 17
)

-Gammarus hyalelloides + G. balmorhea

location_shapes <- c(
  "Phantom Springs" = 16, # Circle
  "BSP San Solomon Pool" = 15,
  "BSP Canal" = 15,
  "BSP Hubbs Cienega" = 15,
  "BSP San Solomon Cienega" = 15,
  "NC West Sandia" = 17, # Square
  "NC East Sandia" = 17
)

• Pyrgulopsis

location_shapes <- c(
  "Phantom Springs" = 16, # Circle
  "BSP San Solomon Pool" = 15,
  "BSP Canal" = 15,
  "NC East Sandia" = 17 # Square
)

• Tryonia

location_shapes <- c(
  "Phantom Springs" = 16, 
  "BSP San Solomon Pool" = 15, 
  "NC East Sandia" = 17
)

Output:

• YOUR_OUTPUTNAME_HERE_PCA_PC1_PC2.pdf: PC1 vs PC2 plot
• YOUR_OUTPUTNAME_HERE_PCA_PC1_PC3.pdf: PC1 vs PC3 plot
• Variance explained percentages for each principal component

Calculate Population Differentiation

Calculate pairwise Nei's GST (multiallelic analogue of Wright's FST):

Rscript GstCode_entropy.R

Parameters (adjust in script):

• ks: specify range of k values to be included, line 7 in the provided script
• chains: keep this at 0:1 if you used 0 and 1 to number your MCMC chains output by Entropy, line 8
• n_ind: number of individuals in the species/analysis, line 9
• n_loci: number of loci being used, line 10
• Change name of output files on lines 89, 90, and 114 of the script
• Change title of plot on line 111

Output:

• YOUR_OUTPUT_NAME_pairwise_Gst_matrix.csv: Pairwise GST matrix
• YOUR_OUTPUT_NAME_pairwise_Gst_CI_matrix.csv: Bootstrap confidence intervals (2.5th and 97.5th percentiles)
• YOUR_OUTPUT_NAME_pairwise_Gst_heatmap.pdf: Visualization

Genetic Diversity Estimation

Calculate population genomic summary statistics using ANGSD (v1.16) to estimate genetic diversity metrics including Watterson's theta and nucleotide diversity (π).

Run ANGSD Theta Estimation

Calculate theta statistics for each population using ANGSD, where pop1, pop2, pop3, are replaced by the appropriate population codes:

bash getThetasAngsd.sh pop1 pop2 pop3 ...

Input:

• Population codes as command-line arguments
• *.sorted.bam: Sorted BAM files for each population need to be in the working directory
• good_reference.fasta: The path to the reference genome needs to be specified on line 10

Output:

• POPCODE.thetas.idx: Theta estimates per window
• POPCODE.thetas.idx.pestPG: Summary statistics per window

Process and Visualize ANGSD Output

Process theta statistics from ANGSD output and generate visualization plots:

Rscript angsd_output.R

Parameters (adjust in script):

• population_list: Population codes in order matching .pestPG files (line 18), see the notes on the parameters for K_mean_barplot.R in Assess MCMC Convergence for details on what the population codes are for each species. However, they will need to be alphabetical here
• ordered_populations: Population order for x-axis in plots (line 42). This will be the previously referenced population codes but in your preferred order, instead of alphabetical.
• Name of output files can be changed on lines 81, 105, and 120

Input:

• *.pestPG files: Theta statistics per window for each population (generated by ANGSD)

Output:

• YOUR_SPECIES_thetaW_per_site.pdf: Boxplot of Watterson's theta per site by population
• YOUR_SPECIES_thetaPi_per_site.pdf: Boxplot of pairwise theta (π) per site by population
• YOUR_SPECIES_theta_summary_per_population.csv: Summary table with mean and median theta values per population

Isolation by Distance Analysis

Test for correlation between geographic distance and genetic differentiation across all species using a Mantel test. This analysis uses pairwise GST values from the Calculate Population Differentiation section above and geographic distances between sampling locations.

Data File

The All_ThreeSpecies_mantel_test_data.csv file contains combined pairwise population data for all three species:

• Tryonia: 3 populations, 3 pairwise comparisons
• Pyrgulopsis: 4 populations, 6 pairwise comparisons
• Gammarus: 7 populations, 21 pairwise comparisons

Data Format:
The CSV file contains four columns:

• species: Species name (Tryonia, Pyrgulopsis, or Gammarus)
• pops: Population pair identifier (format: "pop1 to pop2", e.g., "ph to sp", "es to sp")
• Geographic Distance: Distance in meters between sampling locations (calculated using QGIS Distance Matrix tool)
• GST: Pairwise Nei's GST values (from GstCode_entropy.R output)

Population codes used in data file:

• ph = Phantom Springs
• sp = San Solomon Pool
• cn = San Solomon Canal
• sc = San Solomon Cienega
• hc = Hubbs Cienega
• ws = West Sandia
• es = East Sandia

Run Analysis

We performed Mantel tests using the vegan package for all species simultaneously in RStudio using all_species_mantel_test.R.

Requires:

• All_ThreeSpecies_mantel_test_data.csv file in the working directory

Method:

• Normalization: GST values are linearized using the formula GST / (1 - GST)
• Geographic distances: Calculated between sampling locations using QGIS Distance Matrix tool (units: meters)
• Statistical test: Mantel test with Pearson's product-moment correlation for each species separately
• Permutations: 9,999 permutations to assess statistical significance
• Visualization: Combined scatter plot showing all three species with species-specific colors and regression lines

Output:

• Mantel test results printed to console for each species (Tryonia, Pyrgulopsis, and Gammarus) including:
  • Correlation coefficient (r)
  • p-value
  • Number of permutations
• Combined scatter plot showing:
  • All pairwise comparisons from all three species
  • Separate regression lines for each species
  • Geographic Distance (x-axis) vs. Linearized GST [Gst / (1 - Gst)] (y-axis)

Note: The script automatically processes all species in the CSV file and creates a combined visualization. Each species is analyzed separately with its own Mantel test, but plotted together for comparison.

Acknowledgments

Samples were provided through the work of Chad Norris, Marty Kelly, Pete Diaz, Randy Gibson, and Nina Noreika. Sampling was supported by Texas Parks & Wildlife Department contract #495775, Baseline Biomonitoring of Springs Associated with the San Solomon Springs System. DNA work was funded by a UTRGV College of Sciences SEED grant. DNA extraction and library prep assistance by Dominique Alvear.

For full literature citations, please refer to the manuscript.