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.