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

Author/Principal Investigator Information
Name: Brenton von Takach
ORCID: https://orcid.org/0000-0002-7999-3521
Institution: Curtin University
Email: brenton.vontakach@curtin.edu.au

Date of data collection:
2023 to 2025

Geographic location of data collection:
southwestern Australia

Description of the data and file structure:
This repository contains the data and scripts necessary to reproduce the genomic, spatial, and structural diversity analyses of oblong and long-necked turtles (Chelodina oblonga and C. longicollis). It includes raw occurrences, phenotypic metadata, SNP calls, and the custom bioinformatics pipeline scripts used to go from raw sequence data to finalised results.
This project uses reduced representation sequencing (DArTSeq) from a single round of sequencing of 467 individuals. The sequence data has been uploaded to the NCBI Sequence Read Archive. SNP genotypes were called by DArTSeq Pty Ltd, a commercial facility based in Canberra, Australia. These pipelines are proprietary. Most subsequent analyses were conducted in R, with some command-line scripts used to estimate autosomal heterozygosity. Descriptions of all genotyping and analysis commands and scripts are included as part of the file lists below.

RAW SEQUENCE DATA

The raw sequence data for this project is available from the NCBI (BioProject ID PRJNA1484789), via the following link:
http://www.ncbi.nlm.nih.gov/bioproject/1484789

Note that this bioproject will become publicly available on 2027-07-01

DATA & FILE OVERVIEW

File List:

Folder 'data.zip':
This folder contains all the data used as input into the R and command line scripts.

Report_DChelod25-10599_SNP_mapping_2.csv

This file contains the raw DArTseq output of called SNPs, formatted in the SNP 1-Row Mapping Format. It consists of a top header section describing the sequenced samples, followed by a large data matrix containing marker metadata and the genotype calls for every individual.
The first seven rows of the file define the metadata for the individual samples (the columns on the far right of the spreadsheet). An asterisk * indicates that the row does not apply to the marker metadata columns on the left.

• Row 1: Order number (e.g., DChelod25-10599).
• Row 2: DArT plate barcode.
• Row 3: Client plate barcode.
• Row 4: Well row position (A-H).
• Row 5: Well column position (1-12).
• Row 6: Sample comments (e.g., "[Original species: Chelodina oblonga]").
• Row 7: The primary column headers for the data matrix (Genotype names).

In the main data matrix, each row represents a single SNP marker, and the columns under the sample IDs contain the genotype calls:
"0" = Homozygous for the Reference allele.
"1" = Homozygous for the Alternate/SNP allele.
"2" = Heterozygous.
"-" = Double null / null allele homozygote (i.e., missing data or absence of the fragment in the genomic representation of the sample).

The first 20 columns of the dataset (Row 7 downwards) describe the properties and quality metrics of each sequenced SNP marker:

Column Name: Description

• AlleleID: A unique identifier for the specific allele and its mutation (e.g., 100062515|F|0-52:A>T-52:A>T).
• CloneID: A unique numerical identifier for the sequence fragment in which the SNP marker occurs.
• AlleleSequenceRef: The full nucleotide sequence of the Reference allele.
• AlleleSequenceSnp: The full nucleotide sequence of the alternate/SNP allele.
• TrimmedSequenceRef: The nucleotide sequence of the Reference allele with sequencing adapters removed.
• TrimmedSequenceSnp: The nucleotide sequence of the alternate/SNP allele with sequencing adapters removed.
• SNP: The base position and variant details of the defined SNP (e.g., 52:A>T).
• SnpPosition: The exact position (zero-indexed) in the sequence tag at which the SNP variant base occurs.
• CallRate: The proportion of samples for which the genotype call was successfully scored (i.e., not a "-" missing value).
• OneRatioRef: The proportion of samples for which the genotype score is "1" in the Reference allele data.
• OneRatioSnp: The proportion of samples for which the genotype score is "1" in the SNP allele data.
• FreqHomRef: The proportion of successfully scored samples that are homozygous for the Reference allele (scored as "0").
• FreqHomSnp: The proportion of successfully scored samples that are homozygous for the SNP allele (scored as "1").
• FreqHets: The proportion of successfully scored samples that are heterozygous (scored as "2").
• PICRef: The Polymorphism Information Content (PIC) for the Reference allele.
• PICSnp: The Polymorphism Information Content (PIC) for the alternate/SNP allele.
• AvgPIC: The average Polymorphism Information Content (PIC) across both the Reference and SNP alleles.
• AvgCountRef: The sum of the tag read counts for all samples, divided by the number of samples with non-zero tag read counts, for the Reference allele.
• AvgCountSnp: The sum of the tag read counts for all samples, divided by the number of samples with non-zero tag read counts, for the alternate/SNP allele.
• RepAvg: The proportion of technical replicate assay pairs for which the marker score was identical.

Note: DArT pipelines typically use ~20 % of samples processed from DNA to allelic calls as technical replicates to generate the reproducibility score.

Following the metadata columns (Column 21 onwards), each column corresponds to a specific turtle sample (e.g., CHOBL001, CHOBL127). The data within these columns consists entirely of the genotype calls (0, 1, 2, or -) for that specific individual across all identified SNP markers.

turtles-full-metadata_20251107.csv

This file is the master metadata file. It contains all sample IDs, population assignments, carapace lengths, sex, and geographic coordinates. Note that there are more rows (individuals) than were sequenced. Columns and descriptions are as follows:

Column Name: Description

• commonName: The accepted common name of the turtle species.
• species: The scientific taxonomic classification of the individual (e.g., Chelodina oblonga, Chelodina longicollis, or Chelodina sp. for suspected hybrids).
• id: The unique alphanumeric identifier assigned to each individual turtle sample.
• state: The Australian state or territory where the sample was collected.
• region: The broader geographic region corresponding to the collection site.
• council: The specific local government area or city council jurisdiction of the collection site.
• pop: The local population name or specific waterbody where the turtle was sampled.
• lat: The latitude of the collection site, recorded in decimal degrees.
• lon: The longitude of the collection site, recorded in decimal degrees.
• sex: The in-field identification of the sex of the turtle (M for male, F for female, or blank/NA if undetermined).
• date: The date the sample was collected in the field (formatted as DD/MM/YYYY).
• collector: The name of the researcher or field staff member who captured the turtle and collected the sample.
• collectorID: Any original identification code or tag number assigned to the turtle by the collector in the field.
• tissueType: The specific type of biological tissue collected for DNA extraction (e.g., Muscle, Web).
• storageMedium: The chemical preservative used to store the tissue sample (e.g., Ethanol90).
• class: The taxonomic class of the organism (Reptilia).
• family: The taxonomic family of the organism (Chelidae).
• notes: General field observations, circumstances of capture, or secondary morphological measurements (e.g., plastron width).
• labnotes: Internal comments regarding the laboratory processing, subsampling, or DNA extraction quality of the sample.
• carapace: The straight-line carapace length of the turtle, measured in millimeters.

ALA_occurrences.csv

This spreadsheet contains the raw occurrences downloaded from the Atlas of Living Australia (ALA) for Chelodina oblonga for records from the year 1900 onwards. Columns and descriptions are as follows:

Column Name: Description

• recordID: A unique Universal Unique Identifier (UUID) assigned to the occurrence record by the Atlas of Living Australia.
• scientificName: The accepted scientific taxonomic classification for the recorded organism.
• taxonConceptID: A URL or unique identifier linking to the official Australian Faunal Directory taxonomic concept for the species.
• decimalLatitude: The latitude of the occurrence, recorded in decimal degrees.
• decimalLongitude: The longitude of the occurrence, recorded in decimal degrees.
• eventDate: The date and time when the occurrence was recorded or observed.
• occurrenceStatus: An indicator of whether the species was recorded as present or absent at the location.
• dataResourceName: The name of the original database, citizen science platform, or institution that supplied the occurrence record to the ALA.

targets_227L3CLT1_1.csv

This file serves as the target manifest for the raw sequence data files, mapping sample IDs to their corresponding genotypes for raw read processing. The raw .FASTQ.gz files are found at the NCBI Sequence Read Archive link provided above. Note on Technical Replicates: DArT pipelines will typically use about 20 % of samples processed from DNA to allelic calls as technical replicates. You may notice some duplicate genotype entries in this file; these represent those replicates used to calculate reproducibility metrics. Columns and descriptions are as follows:

Column Name: Description

• targetid: A unique internal identifier assigned by DArT to each individual sequenced target in the run.
• ordernumber: The overarching client order identifier under which the samples were processed (e.g., DChelod25-10599).
• organism: The broad taxonomic group or genus of the samples (e.g., Chelodina).
• species: The species classification, if specified prior to processing (often blank or "-" if handled downstream).
• genotype: The sample identifier or ID provided by the client (e.g., CHOBL001). This column is crucial for mapping sequenced reads back to the main metadata file.
• tissue: The biological tissue type used for the DNA extraction (e.g., skin).
• aliastargetname: An internal alias or alternative processing name assigned to the target (e.g., P_ad_hp).
• primersname: The specific combination of restriction enzymes and primers used for the DArTseq genomic complexity reduction (e.g., EBPCR1+HpaII).
• extractplatebarcode: The unique barcode or physical identifier of the plate used during the initial DNA extraction.
• extractplatewell: The alphanumeric well position (e.g., A1, B1) on the extraction plate where the sample was placed.
• targetplatebarcode: The unique barcode or physical identifier of the plate used for target sequencing preparation.
• targetplatewell: The alphanumeric well position on the target sequencing plate.
• label: Optional column for additional sample labeling (frequently blank/NaN).
• comment: Optional column for specific laboratory comments or flags related to the individual target (frequently blank/NaN).
• barcode9l: The full nucleotide barcode sequence used to uniquely tag and demultiplex the sample's sequence reads.
• barcode: The trimmed or primary nucleotide barcode sequence for the sample.
• extractid: A unique numerical identifier assigned specifically to the physical DNA extraction event.
• tagcounttotal: The total number of sequence reads (tags) successfully generated for this sample.
• tagcountunique: The total number of unique sequence reads (tags) generated, removing identical duplicate reads.
• overrepcount: The raw count of specific sequence tags that are highly over-represented in the sample's sequencing yield.
• overreppct: The percentage of the sample's total reads that are considered over-represented.
• targetquality: A qualitative assessment provided by DArT regarding the sample's sequencing success (e.g., "good").
• flowcellbarcode: The unique identifier for the specific Illumina sequencing flow cell used for the run (e.g., 227L3CLT1).
• flowcelllane: The specific physical lane (e.g., 1) on the flow cell where the sample pool was sequenced.

Note: blank cells used to represent missing data or "not applicable"

Folder 'scripts.zip':

Subfolder 'commandline':

Several command line scripts handle the automated processing of raw sequence data into finalized variant calls. These scripts were only used to obtain autosomal heterozygosity values from sequence data (as opposed to being calculated from the called SNPs).

1_merge.sh

This script extracts TargetID and Genotype data from a target CSV file, automatically generates a Stacks population map, and merges individual raw fastq files matching the target IDs into combined per-genotype gzipped fastq files.

2_clean.sh

This script runs the process_radtags module to clean uncalled bases, check quality scores, and enforce a strict uniform length by trimming all gzipped reads to exactly 80bp.

3_stacks.sh

This script executes the initial denovo_map.pl assembly and population modules within Stacks to generate a preliminary VCF file from the cleaned fastq files.

4_select_balanced_catalog.py

This Python script parses the master metadata file to select a balanced representation of individuals across sixteen target populations, prioritizing high-depth samples based on file size to output a balanced catalog popmap capped at four samples per population.

5_stacksagain.sh

This script resumes the Stacks denovo_map.pl pipeline by leveraging the balanced catalog popmap to save processing time, subsequently running the populations module to produce an updated VCF.

6_run_batches.sh

This script divides the primary population map into distinct chunks of 93 individuals, loops through each chunk to execute gstacks and the populations module independently, and saves the final isolated VCF and haplotype files for each batch.

Subfolder 'R':

The downstream genomic, spatial, and statistical analyses are executed via custom R scripts.

turtles-range-wide-dart-prep-20250116.R

This script performs rigorous filtering and preparation on the primary genomic dataset, filtering out non-target species and cleaning markers based on read depth, individual and locus call rates, minor allele frequency, reproducibility, linkage disequilibrium, and sex-linked loci.

turtle-ind-het-20251127.R

This script reads the multi-batch VCF (produced by the above Bash and Python scripts) and haplotype data to calculate individual autosomal heterozygosity, generates population-level boxplots, tests correlation with genomic SNP heterozygosity, and builds linear mixed-effects models evaluating relationships with carapace length and sex.

turtle-distribution-polygon-20251124.R

This script filters raw Atlas of Living Australia occurrence data, appends additional verified field records, builds a spatial concave hull (alpha hull), and smooths and masks the shape to the southwestern Australian coastline to export a finalized distribution geotiff (for use in the range-wide mapping script).

turtle-range-wide-struct_div_20251116.R

This script conducts comprehensive population structure and diversity analyses by calculating metric MDS coordinates, executing cross-validation landscape admixture runs in TESS3, computing pairwise population FST matrices with Mantel tests, and estimating within-population relatedness.

turtle-range-wide-mapping-20251117.R

This script builds the primary range-wide sampling maps with custom inset layouts using ggplot2, projects species distribution extents, and runs generalized linear models tracking historical annual rainfall trends from Bureau of Meteorology observations.

turtles-hybridisation-20251118.R

This script focuses on hybrid detection between species by assessing locus filtering behaviors, computing individual hybrid indices using Bayesian MCMC diagnostics, and classifying specific hybrid classes through polarized fixed-difference allele frequencies.

Additional workarounds for complete reproducibility:

Due to repository licensing requirements (CC0), certain third-party proprietary mapping and spatial files used in the original analyses are omitted. Researchers looking to reproduce these exact steps should use the following open-access alternatives.

Australian Boundary Shapefile (IGISmap): The mapping script originally utilized a locally hosted shapefile of the Australian coastline from IGISmap. To reproduce these maps without violating copyright, you can dynamically fetch public domain digital boundaries using the rnaturalearth R package. Replace the local file import step with code that loads the rnaturalearth and sf libraries, fetches the map of Australia dynamically using the ne_states function, and converts it to a spatial polygon.

Perth Rainfall Data (BOM): The mapping script evaluates historical rainfall trends using a file originally downloaded from the Bureau of Meteorology. You can freely download the equivalent raw data directly from the BOM Climate Data Online Portal. Search for the Perth Airport station (009021), download the rainfall data, and save it as rainfall_perth_20251204.csv in your data directory.

Environmental Temperature Raster: The distribution polygon script originally used a local raster file (annualmeantemp.tif) to convert the alpha-hull spatial polygons into a finalized .tif distribution raster. This file simply acts as a high-resolution template to define the grid extent and resolution. You can download any standard open-access bioclimatic raster, such as BIO1 Annual Mean Temperature, from WorldClim at your preferred resolution and substitute it into the script.

External Software Dependencies: The hybridisation and filtering scripts make command-line calls to several external bioinformatics programs. To run these successfully, you will need to download and install Stacks for RADseq processing, PLINK for formatting VCFs, and STRUCTURE for Bayesian clustering analysis. Ensure you update the local file paths in the R scripts to match the location of executables/binaries on your system.