How the Salamander Example Data Were Generated

Two examples, three kinds of data

Both salamander examples use real southeastern Ohio terrain and land cover, but the animals, genetic data, pond locations, and drying index are simulated.

The statistical examples reuse the landscape but contain no salamander genotypes. They illustrate a method without testing whether its specified relationship describes salamander biology.

Prepare the real landscape

Elevation comes from USGS 3DEP, and land cover from the 2023 Annual NLCD product, Collection 1.0. Inputs were aligned in NAD83 / UTM zone 17N (EPSG:26917), so distances are in meters. The installed provenance files listed below record source requests and file hashes.

The pond analysis retains a 10.2 by 10.2 km landscape at 100 m resolution. Forest fraction is the proportion of finer land-cover cells classified as deciduous, evergreen, or mixed forest, or woody wetlands. Elevation and slope are retained as alternative predictors, although only forest controls pond exchange in the simulation. A surrounding input area reduces edge effects in forest context.

The terrestrial example uses an 840 by 840 m part of the same landscape, with a 300 by 300 m sampling core and a 270 m buffer. Its forest mask uses deciduous, evergreen, and mixed forest. Drying is a hypothetical index:

\[D=0.60(1-W)+0.40X,\]

where \(W\) is the scaled terrain-wetness quantity and \(X\) is scaled exposure. Nonforest native cells receive \(D=1\). These are assumed values, not measured water loss. The 10 m analysis grid interpolates the native 30 m information so nearby individuals can occupy separate graph cells.

library(terradish)
library(terra)
#> terra 1.9.50
data(ohio_ponds)
data(individual_salamanders)
pond_landscape <- unwrap(ohio_ponds$covariates)
files <- system.file("extdata", "individual-salamanders",
                     package = "terradish", mustWork = TRUE)
drying <- rast(file.path(files, "environment_10m.tif"))

local({
  old <- par(mfrow = c(1, 2), mar = c(3, 3, 2, 3))
  on.exit(par(old))
  plot(pond_landscape[["forest"]], main = "Pond forest fraction", cex.main = 0.9)
  points(ohio_ponds$ponds[c("x", "y")], pch = 20, cex = 0.5)
  plot(drying[[1]], main = "Terrestrial drying", cex.main = 0.9)
  samples <- individual_salamanders$data$compensatory_500$samples
  points(samples[c("easting", "northing")], pch = 20, cex = 0.3)
})

Real landscape inputs for the two simulated examples: hypothetical ponds on the forest-fraction map (left) and the assumed terrestrial drying index (right).

How to read the maps. The examples use different extents and sampling units. Points mark sampled ponds or adults, and every cell between them, including nonforest cells, remains available to the graph.

Generate the pond history in SLiM

All 24 hypothetical ponds lie in native forest. Screening required at least 90% forest within 100 m, at least 60% within 300 m, slope no greater than 12 degrees, spacing of at least 600 m, and a 1,500 m domain-edge margin. These were simulation design rules, not habitat thresholds estimated from data.

Forest fraction was smoothed with a Gaussian kernel of standard deviation 500 m. Generating cell conductance was

\[C_i=\exp\{6(F_{i,500}-\overline F_{500})\}.\]

Here \(F_{i,500}\) is smoothed forest fraction on its original 0-to-1 scale. Eliminating the intervening nodes of an eight-neighbor cell graph gives an exchange network among ponds that retains the graph’s effective connections. Total exchange probability per pond is at most 0.025, so most reproductive events use local parents, with a small long-distance tail. The model records reproductive exchange between ponds, not travel tracks.

SLiM 5.2 ran a neutral model for 1,200 nonoverlapping generations. Pond size varies around a target of 120 diploids rather than staying fixed as in a Wright-Fisher model. One parent can come from the exchange distribution, and the other is local. Hydroperiod is a hypothetical pond attribute with no genetic effect in this history, which makes it a negative control for isolation by environment (IBE). The packaged history (seed 79101) was fixed before any analysis.

Generate the terrestrial history in SLiM

This model follows individuals annually for 500 years, with local mating, crowding, repeated breeding, movement, and survival. Two scenarios share parameter values and seed 90401. The control fixes the drying input to the movement functions at 0.5, while the compensatory scenario uses local drying. The scenarios develop different individuals and pedigrees, so their animals are not matched.

Movement functions interpolate linearly between the following drying values:

Heading persistence is the mean cosine of the turning angle, so values near one produce straighter paths. Exposure hazard per minute is \(0.001+0.08\operatorname{mean}(D^8)\), evaluated at four points along each step. Settlement requires forest, and a surviving animal still outside forest at the 60-minute limit dies. Each movement event represents an episodic dispersal opportunity, not all daily foraging or surface activity.

The resulting hypothesis is that quicker, straighter travel can compensate for moderate drying, while very dry terrain remains dangerous. It was motivated by the interpretation in Peterman et al. (2014), not estimated from that study. Because movement, survival, settlement, and later reproduction jointly determine realized gene flow, fitted conductance is not identical to the speed function used here.

Select animals before examining genotypes

Pond sampling began with 30 animals per pond, adding captures where needed to retain 20 after screening out known immediate pedigree relatives within each pond. The saved history contains 780 raw captures and 480 retained animals.

Terrestrial sampling drew 160 adults alive in year 500 from sixteen 75 m strata, 10 per stratum, with 40 in each 150 m quadrant. Every adult is at least five years old, at least 10 m from the others, and in a distinct 10 m cell. Direct and shared recorded parents and grandparents were excluded.

table(ohio_ponds$samples$pond)
#> 
#> P01 P02 P03 P04 P05 P06 P07 P08 P09 P10 P11 P12 P13 P14 P15 P16 P17 P18 P19 P20 
#>  20  20  20  20  20  20  20  20  20  20  20  20  20  20  20  20  20  20  20  20 
#> P21 P22 P23 P24 
#>  20  20  20  20
lapply(individual_salamanders$data, function(x) table(x$samples$stratum))
#> $compensatory_500
#> 
#>  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 
#> 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 
#> 
#> $control_500
#> 
#>  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 
#> 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10

How to read the counts. Each pond has 20 animals, and each terrestrial scenario has 10 adults in each of 16 strata. Both marker types use exactly these animals. Kin screening used the recorded pedigree before genotypes were generated, and in both examples more distant shared ancestry can remain. With field data, relatedness must itself be estimated.

Add ancestry and matched markers

SLiM supplies the recent ecological history, and tree-sequence tools then add ancestry and mutations. A tree sequence stores how sampled genomes share ancestors along a chromosome, so several marker panels can come from one history without resimulating the animals.

For ponds, the simulated chromosome is 50 Mb, with recombination rate \(2\times10^{-8}\) per base per generation. Older ancestry is completed using population size 120 and half the forward exchange matrix as a backward-time approximation, an assumption that does not prove equilibrium. The seeds for completing older ancestry (recapitation), SNP mutation, and microsatellite mutation are 79102, 79103, and 79104.

For terrestrial animals, msprime reconstructs inheritance through the saved SLiM parents for 20 synthetic 2.5 Mb blocks. Recombination is \(2\times10^{-8}\) per base per meiosis. Earlier ancestry assumes a randomly mating population with effective size 5,000. Calendar years are converted to approximate generations using mean parental age among recruits in years 401-500, a constant clock that simplifies overlapping generations. Genomes are those inherited at birth, with no somatic mutation added before sampling.

The panels share sampled animals and ancestry within each scenario, but their mutation processes differ, and SNP spacing does not make all loci independent. The panels omit restriction sites, read depth, missing calls, and genotyping errors, so they are not complete RADseq experiments.

A second terrestrial assay generated new genomic realizations from the same pedigrees after model parameters were frozen. It checks whether predictions persist with new markers but adds no new animals or independent ecological history. The individual guide reports both assays and the mixed microsatellite evidence.

From genotypes to the modeled response

The pond guides aggregate allele counts by pond and use a covariance that weights loci by pooled allele frequency. The individual guide instead uses unstandardized allele-frequency profiles averaged across loci, so the two responses differ in numerical scale and interpretation. vignette("wishart-covariance", package = "terradish") covers the constructors, their mathematical checks, and when to use a covariance or a pairwise distance. Always preserve sample order across genotypes, coordinates, and environments.

The pond demonstrations assume nu = 40, and the individual example uses nu = 1 as a covariance-fitting objective. Neither value is a SNP count.

Reproduce an analysis or regenerate a history

The package contains the genotypes, landscapes, settings, and results needed to reproduce the documented analyses, but not the full ancestry archives needed to regenerate either history.

# Read the installed source records and data dictionaries.
file.show(system.file("extdata", "ohio-ponds", "PROVENANCE.md",
                      package = "terradish"))
file.show(system.file("extdata", "individual-salamanders", "README.md",
                      package = "terradish"))

# Copy an executable analysis source before rendering it.
file.copy(system.file("vignette-source", "getting-started.Rmd",
                      package = "terradish"), "getting-started.Rmd")
rmarkdown::render("getting-started.Rmd")

Routine package builds use precomputed figures and results for long analyses. The individual guide’s opt-in refit.R script saves each attempted fit and can resume. Regenerating the ancestry requires that undistributed archive plus the same software (SLiM 5.2, with msprime 1.4.2 and tskit 1.0.3 for marker processing).

See also

• vignette("getting-started", package = "terradish"): pond genotypes to fitted forest effects.
• vignette("individual-salamanders", package = "terradish"): nonlinear IBR and individual prediction.
• vignette("simulation-design", package = "terradish"): repeated statistical draws for planning, distinct from SLiM histories.
• vignette("wishart-covariance", package = "terradish"): turning genotypes into a covariance or distance response.
• ?ohio_ponds and ?individual_salamanders: the complete data dictionaries.

Reference

Peterman WE, Connette GM, Semlitsch RD, Eggert LS (2014). Ecological resistance surfaces predict fine-scale genetic differentiation in a terrestrial woodland salamander. Molecular Ecology, 23. doi:10.1111/mec.12747.