Combining “landscape” and “genetics”
To investigate whether Isolation by distance (IBD) varies among lineages
or regions we fitted a linear model between pair-wise log Euclidean
geographic distance and genetic distance among individuals, calculated
by “1 minus the proportion of shared alleles between each pair of
individuals” using the PropShared function from theadegenet package (Jombart & Ahmed, 2011). The strength of IBD
was measured within each of the taxa (species or intraspecific lineages)
as identified above. The IBD slope is expected to increase with
decreasing local population size and/or migration rate (Rousset, 1997)
but is also reduced at large spatial scales by recent population
expansion (Slatkin, 1994). Mantel tests were performed to calculate a
p-value. We visualized slopes using the R package jtools (Long,
2019), and tested for differences across taxa (skinks, geckos, and
agamids), environments (more below), and between rock specialists vs
habitat generalists. To represent climatic variation we focused on mean
annual rainfall values, the dominant environmental variable of the
region, from WorldClim (Hijmans, Cameron, & Parra, 2006) extracted from
the convex hull polygons created from each species distribution
locations, as the dominant environment variable of the region. To
account for the influence of the precipitation, habitat preference and
lizard group we performed a generalized linear model (glm; as our
variables were non-normal) with the slope value of the IBD analysis and
using 1-standard error as weighting.
To investigate the effect of environmental features on genetic
difference among individuals (Isolation by Environment; IBE), we
downloaded the 19 variables from WorldClim (Hijmans et al., 2006) at
30-s resolution. We tested for correlation among these variables using
the ENMTools R package (Warren et al., 2019), keeping eight
uncorrelated (r < 0.7) variables (Annual Mean Temperature,
Mean Diurnal Range, Isothermality, Mean Temperature of Wettest Quarter
and Driest Quarter, Precipitation Seasonality, Precipitation of Driest
Quarter and Warmest Quarter) and extracting these variables for each
locality per species. The value of each environmental variable per
distribution point was used later as an environmental variable to relate
to genetic divergence using generalized dissimilarity modelling,
discussed above.
To avoid over-sampling our Isolation by Resistance (IBR) models we
excluded occurrence records that were in islands (no-connectivity) or
within the same resolution grid square at ~1km
resolution. The resulting numbers of individuals used for each taxon are
indicated in Table 1. We considered four landscape features to test the
effect of resistance surfaces on species divergence: rivers, slope,
evapotranspiration and potential environmental suitability based on
species distribution models (SDM). Slope and rivers can represent
physical and physiological barriers to dispersal for some of the
species. We obtained a river shapefile from Geoscience Australia (1997),
and considered major rivers as high resistance cells and intervening
regions as low resistance. We extracted slope from altitude (CGIAR-CSI,
2019), using the package raster (Hijmans, 2019). Slope was
obtained from the weighting coefficient for the altitude for each grid
cell, helping to recognize flat areas that could facilitate dispersal of
habitat generalist species or inhibit dispersal for rock-dependent
species. All asci files were normalized from 0 to 1, with greater values
indicating higher resistance. Global potential Evapotranspiration was
obtained from CGIARC-CSI, 2019
(https://cgiarcsi.community/data/global-aridity-and-pet-database) and
represents the water/productivity available for the organisms, which
allows contrasting the mesic and dry conditions across the general
north/south aridity gradient of the AMT (Figure 1A). Finally, we created
SDMs for all species and lineages independently, based on the locality
of genetically identified samples (see complete details in the
Supplementary Material S1). We included all 19 WorldClim variables,
running the Python 2.7 script (available at
github.com/DanRosauer/phylospatial) to automate data preparation and
modeling in ArcGIS 10 (ESRI, 2011) and MaxEnt v3.3predictive approach (Phillips, Dudik, & Schapire, 2012). Following
Elith et al. (2011), we retained even correlated variables, allowingMaxEnt to determine the required predictors for each model. To
avoid over-prediction, background points were limited to a 3 degrees
radius around the location records, with 25 bootstrap replicates for
each SDM, and the result was exported as an asci file and normalized.
The area under the curve (AUC) values for the SDMs vary between 0.95 and
0.71, indicating good model performance, except for the Kimberley
lineage of Diporiphora magna , with an AUC of 0.61 (Table S1).
Both precipitation and temperature related to the seasonality of the AMT
appear to be important variables in most SDMs (Table S1). We used
circuit theory to identify the most effective path by using the four
landscape features described above in the software Circuitscape(McRae, Shah, & Mohapatra, 2013). This looks for all paths and averages
them according to resistance values, where the SDM was used as a
conductance surface (where higher values indicate greater ease of
movement) and rivers, slope and evapotranspiration as resistances.
To investigate how IBD, IBE, and IBR contribute to genetic divergence,
we based our analysis in part on the framework of Myers et al. (2019),
using generalized dissimilarity modelling (GDM; Ferrier, Manion, Elith,
& Richardson, 2007). This method calculates the relationship between
nonlinear models associating the variation in the rate of genetic
dissimilarity among individuals with distance variables. We ran seven
independent tests for each species: first a full model, including
geographic distance, environmental variables, and resistance surfaces;
then, looking for the interaction between two of the datasets; and
lastly looking at each dataset alone (more details in Table 1). We used
genetic divergence as the response variable and the geographic,
environmental and resistance datasets as predictors in a GDM site-pair
table, and quantified the predictor’s importance and significance using
a GDM matrix with 50 permutations.