Frogs were sampled between the months of March and July in 2013, 2014, and 2015. We collected a total of 191 P. auritus from four sites in Cameroon (Campo Ma’an (CM), N = 29; Ebo forest (EF), N = 24; Ndikiniméki (ND), N = 13; and Takamanda (TM), N = 15), and five sites in Gabon (Gamba Complex (GC), N = 24; Kessala (KS), N = 21; Lopé (LP), N = 18; Monts de Cristal (MC), N = 22; and Minkébé (MK), N = 25) ( Figure 1A ). These sites encompassed a range of forest types, namely: lowland rainforest (CM, MK), sub-montane rainforest (EF, TM), forest-savanna ecotone (KS, MC), coastal rainforest (GC), and mixed lowland-agricultural forest (ND). After frogs were euthanized with MS222 solution, muscle and kidney tissue was placed in 95% ethanol prior to DNA extraction. Male frogs were differentiated from female frogs by the presence of a developed vocal sac, indicated by a dark throat and vocal folds. Frogs that could not be sexed by these secondary sexual characteristics were dissected and sexed by either the type of reproductive organ or by the presence of eggs. All animal handling procedures were carried out according to an approved University of New Orleans Institutional Animal Care and Use Committee protocol 12-008. Specimens were deposited at the Museum of Comparative Zoology, Harvard University.

We used five environmental variables to predict the influence of environmental variation on genomic differentiation given their relevance to amphibian ecology and biology as well as their relevance in previous studies (Ficetola and Maiorano, 2016; Morgan et al., 2020). The five bioclimatic variables were: annual temperature, temperature seasonality, annual precipitation, precipitation seasonality, and precipitation of the coldest quarter from the WorldClim database (Fick and Hijmans, 2017) (www.worldclim.org). Precipitation of the coldest quarter was selected because it represents a unique characteristic of precipitation for this region. It effectively reflects the seasonal inversion that occurs approximately across the equator such that the dry season in central Cameroon coincides with the rainy season in northern Gabon and vice versa (Heuertz et al., 2014; Morgan et al., 2020). We then tested for potential correlations between environmental variables to ensure no two variables had a Pearson’s correlation coefficient > 0.8. All climate variables had a resolution of 30 arc-seconds (approximately 1 km²). Although estimates of vegetation cover may also be important determinants of genomic differentiation, the lack of climate projection maps prevents the use of these variables for modeling under future conditions and were therefore not included in the current study. Future projections of bioclimatic variables were taken from aggregated global climate models (Clee, 2017) for two representative concentration pathways (RCPs) 2.6 and 8.5, projected for year 2080 based on the Intergovernmental Panel on Climate Change (IPCC) 5th assessment report. Together, these two RCPs span the greatest range in predicted radiative forcing under climate change and represent “best” and “worst” case scenarios regarding global mean temperature increases.

We extracted genomic DNA from either kidney or muscle tissue using Qiagen DNEasy Blood and Tissue kit (Qiagen, CA), following the manufacturer’s protocol. A total of 164 individuals had genomic DNA of sufficient quantity (> 50 ng) and quality needed for restriction site-associated sequencing (RAD-seq) as determined through gel electrophoresis (Davey et al., 2011). RAD-seq library preparation followed the BestRad protocol for Illumina sequencing as described in Ali et al. (2016). Briefly, genomic DNA (100 ng) was digested with 4.8 units of SbfI-HF restriction enzyme (New England Biolabs NEB, R3632L) at 37°C for 1 h in a 12 µl reaction volume. Samples were heated to 65°C for 20 min and 4 µl of the indexed BestRad SbfI P1 RAD adapter (10 nM) was added to each sample. Ligation of inline barcoded P1 adaptors to digested genomic DNA was performed overnight at 20°C with 640 units of T4 DNA ligase (NEB, M0202M), then 65°C for 20 min. Following ligation, 10 µl of each sample in each 48 well plate was pooled into a single tube and cleaned using 1x Agencourt AMPure XP beads (A63881; Beckman Coulter). Pooled DNA for each plate was then resuspended in 100 µl low TE and sheared to an average fragment size of 500 base pairs using a Bioruptor NGS sonicator (Diagenode). Sheared DNA was then concentrated to 55.5 µl using Ampure XP beads and used as the template in the NEBNext Ultra DNA Library Prep Kit for Illumina (NEB E7370L; v.1.2). The standard NEBNext protocol for library preparation was followed except that we used custom P2 adaptors which were created by annealing a NEBNext Multiplex Oligo for Illumina (NEB, E7335L) to the oligonucleotide GATCGGAAGAGCACACGTCTGAACTCC AGTCACIIIIIIATCAGAACA*A (where * represents a phosphorothioate DNA base). In addition, instead of the USER^® enzyme step, we used a universal P1 RAD primer (AATGATACGGCGACCACCGAGATCTAC ACTCTTTCCCTACACGAC*G) and a universal P2 RAD primer (CAAGCAGAAGACGGCATACG*A) during final amplification. The final RAD library was cleaned using AMPure XP beads and sequenced at the UC Berkeley QB3 Vincent J Coates Genome Sequencing Laboratory (GSL) on an Illumina HiSeq2500: Rapid Run Mode (Illumina, San Diego, CA, USA) using paired-end 100-bp sequence reads.

We used the bioinformatics software pipeline, STACKS v.1.44 (Catchen et al., 2011; Catchen et al., 2013) to process the restriction-site-associated DNA markers (RAD-tags) and generate single nucleotide polymorphism (SNP) datasets. First, we executed the “process_radtags” program in STACKS to demultiplex and trim sequence reads by the P1 barcodes and remove low quality reads (Phred quality score less than 20). After removing PCR duplicates with the “clone_filter” script, the processed reads were used to generate RAD loci without a reference genome using “denovo_map.pl” (parameter settings: m = 3 M = 5 n = 4). We empirically determined these parameters to limit the impact of over-splitting loci (Ilut et al., 2014; Harvey et al., 2015). This involved running the de novo assembly over a wide range of values of M (1–8) with “ustacks”. From these runs, we selected a value of M = 5 since we observed that the percentage of homozygous and heterozygous loci reached a stable value and thus minimized over-splitting of alleles for the final SNP calling.

Stacks calls SNPs (“sstacks”) within RAD loci using a multinomial-based likelihood model that estimates the likelihood of the two most frequently observed genotypes at each site and performs a standard likelihood ratio test using a chi-square distribution (Hohenlohe et al., 2010; Catchen et al., 2011). For SNP inference, we used the default alpha significance level of 0.05. Paralogous loci that stacked together were identified and removed by subsequent quality control steps built into STACKS (max number of stacks per loci (m) = 3 (Ilut et al., 2014; Harvey et al., 2015)). After the preliminary assembly of catalog loci using “denovo_map.pl”, we ran the STACKS correction mode (rxstacks-cstacks-sstacks) using the bounded SNP model with a 0.05 upper bound for the error rate. The “rxstacks” program made corrections to genotype and haplotype calls based on population information, rebuilt the catalog loci and filtered out loci with average log likelihood ratio of < 8.0.

We used three additional filtering steps to generate a set of high-quality RAD loci for downstream population genetic analysis. First, we retained only RAD loci that were present in 80% of all samples. Second, we removed RAD loci that contained more than 40 SNPs, as these likely represented sequencing errors or over-clustering of paralogous loci. Lastly, we used the BLAT alignment algorithm (Kent, 2002) to de novo align the RAD loci and remove those that aligned to multiple positions. The final consensus set of RAD loci comprised SNP data from a total of 139 individuals. Genotypes were called, filtered, and bi-allelic SNPs were exported in VCF format using the STACKS “populations” program. SNPs from the last seven bp of the RAD loci were removed as this part of the locus is likely to contain sequence errors at the 3’ end of the reads. The SNP dataset was further filtered with VCFtools v.0.1.14 (Danecek et al., 2011) to remove SNPs below a minor allele frequency (MAF) of 0.05 cutoff to reduce artifacts of sequence and assembly error. Finally, the dataset was filtered to include only one random SNP per RAD locus for use in all downstream analyses in order to avoid linkage disequilibrium between SNPs within RAD loci.

We performed a principal component analysis (PCA) using the Bioconductor package SNPRelate (Zheng et al., 2012) (https://www.bioconductor.org/) to summarize population genomic structure. We used the program FastStructure (Raj et al., 2014) to estimate the number of genetically distinct populations within the sampled P. auritus range. We tested a range of K values (where K denotes the number of inferred populations) from 1 to 10. The script “chooseK.py” included in the FastStructure package was used to determine the best estimate of K that maximizes the marginal likelihood. We also calculated pairwise estimates of F_ST (Weir and Cockerham, 1984) among sites and among K populations inferred from FastStructure using VCFtools. To test for an IBD effect, a Mantel test was used to assess the correlation between pairwise F_ST values and geographic distance. Mantel tests were run with 999,999 permutations using the vegan package (Oksanen et al., 2022) in R and are reported using both raw F_ST and transformed F_ST/(1-F_ST) distances, as well as both raw Euclidian geographic distance and log-transformed Euclidean distances (Slatkin, 1995; Rousset, 1997).

We used GDM to compare the importance of IBE to isolation by landscape barriers (IBB) or Pleistocene refugia (IBP) on patterns of genomic differentiation. GDM is a matrix regression technique that evaluates the relationship between site-site dissimilarities in environmental or landscape ‘predictor’ variables and a biotic ‘response’ variable (e.g. pairwise genetic distances). A major advantage of GDM over other modeling methodologies is that it can fit non-linear relationships between environmental variables and the biological response variable through the use of I-spline basis functions (Ferrier et al., 2007). This approach includes straight-line geographic distance as a predictor variable and can also incorporate a range of environmental data layers, and matrices derived from resistance surfaces as different predictors.

Pairwise dissimilarity in genomic composition between sites was modeled using two measures: 1) pairwise F_ST values and 2) a pairwise Bray-Curtis dissimilarity index based on the presence or absence of a SNP at each locus referred to as the Bray-Curtis allele frequency difference or AFD (Sherwin, 2022). IBE was represented by the set of five environmental variables described previously. In addition to these environmental variables, a set of predictor variables were generated to model the effect of IBB and IBP under the Last Glacial Maximum approximately 21,000 years ago (i.e. IBP). IBB represents physical barriers (elevation and rivers) to gene flow. Pairwise resistance distances for IBB were generated by creating raster layers of resistance surfaces based on landscape features, elevation and rivers, using the raster calculator available in QGIS v.2.18. We then calculated pairwise resistance distances from these raster layers with CIRCUITSCAPE 4.0 (McRae et al., 2013). Two IBB matrices were generated, IBB1 and IBB2. For IBB1, resistance values increased with increasing elevation and major rivers were treated as impenetrable. For IBB2, resistance increased with increasing elevation and also with Strahler order, which reflects size and strength of perennial river systems. IBP represents the historical landscape based on modeled available habitat. For IBP, we first projected habitat suitability for P. auritus under climate conditions during the LGM using two global climate models (CCSM and MIROC). We then created resistance surfaces where resistance was considered to be inversely proportional to habitat suitability, and finally, calculated pairwise resistance distances from this raster layer with CIRCUITSCAPE. Further details on how these predictor variables were generated and the resultant distance matrices can be found in Appendix 1 . I-spline turnover functions describing the relationship between the biological response variable (pairwise dissimilarity in genomic composition) and each of the predictor matrices were visualized using rug plots and their significance was tested using 1000 permutations.

We evaluated the importance of environmental variables as predictors of environmentally-associated genomic turnover and spatialized these patterns across the study region using GF modeling within the gradientForest package (Ellis et al., 2012) in R ( Appendix 1 ). Response variables were individual SNP minor allele frequencies within each population. Predictor variables were represented by the same environmental variables that were included in the GDM along with latitude and longitude. GF uses a machine-learning algorithm to divide the biological data into different bins (i.e. different values of allele frequencies), with partitions occurring at several split values along each environmental variable. This binning is performed for every SNP, weighting each SNP individually according to its fit to the model (i.e. R²) before aggregating across all SNPs. GF determines the “split importance” by measuring the amount of biological variation explained by a given split value (e.g. between 26 and 27°C), which is then cumulatively summed along each gradient to construct turnover functions (Fitzpatrick and Keller, 2015). The top three environmental variables in modeling genomic turnover from a total of 2000 regression trees were used to predict and map environmentally-associated turnover across the study region using a random grid of 100,000 sample points. To ensure that our GF model was performing better than random, we shuffled the environmental-predictor matrix to generate 200 randomized datasets and compared the number of SNP loci with R² positive values to the mean R² value across SNP loci using GF models describing variation in the real versus randomized datasets.

Lastly, we predicted future environmentally-associated genomic variation based on projecting GF models under climate change for the years 2050 and 2080 for RCPs 2.6 and 8.5, representing “best” and “worst” cases, respectively. To map predicted changes of genomic variation associated with environment, we calculated the Euclidean distance between models based on current and future climate conditions. Areas where environmentally-associated genomic variation changes the least are considered to have low genomic vulnerability whereas areas where they change the most are considered to have high genomic vulnerability.

RAD-sequencing of 139 P. auritus samples generated a total of 838,425,400 paired-end reads across both plates (308,255,191; 530,170,209) after filtering out low-quality samples and reads. The number of raw sequencing reads per sample ranged from 133,185 to 16 million. The mean coverage depth ranged from 5x to 26x across individual samples (mean = 8x, median = 7x, Appendix 2 ). From these reads, we assembled 2,979 high-quality RAD loci and a total of 32,966 SNPs that were present in 80% or more samples. Using a minor allele frequency cutoff of 5%, we retained 1631 RAD loci encompassing 3,092 SNPs, and from which we selected one random SNP per RAD locus.

The PCA identified significant population structure across the sampled range of P. auritus. PC1 explained 25.89% of the variation and separated the three northern sites (EF, ND, TM) from the remaining six sites (CM, MC, MK, LP, KS, GC). PC2 explained 9.26% of variation and separated GC, the southernmost coastal site, from all other sites ( Figure 1B ). FastStructure analyses also revealed a pattern of population structure that is organized latitudinally into five distinct populations: 1) EF, ND, and TM, 2) CM, 3) MC and MK, 4) KS and LP, and 5) GC ( Figure 1C , Appendix 3 ).

Pairwise F_ST values between sites ranged from 0 to 0.438 (mean = 0.229; Appendix 4 ), indicating low to moderate levels of genomic differentiation between sites. We found a significant correlation between pairwise F_ST and geographic distances between the sites (Mantel r = 0.6728; mantel simulated p-value = 0.001), suggesting a strong pattern of IBD.

The GDM based on pairwise F_ST values and AFD explained 80% and 81% of the variation in the data, respectively. Precipitation of the coldest quarter was the only significant variable (p = 0.02) and the most important variable in the model with F_ST values. Geographic distance was the only significant variable in the model based on AFD (p = 0.01), however precipitation of the coldest quarter was the most important variable, followed by geographic distance ( Table 1 ). The I-spline plot for precipitation of the coldest quarter shows that observed compositional turnover in genomic differentiation increases with greater levels of precipitation ( Supplementary Figure S1 ). Importantly, neither of the two resistance matrices (IBB1, IBB2) modeling the effects of riverine barriers nor the two matrices modeling the distribution of suitable habitat since the Pleistocene (IBP-CCSM, IBP-MIROC) had any significant effect on the model (i.e. all matrices had a coefficient = 0 except IBP-MIROC).

We used a GF approach to determine associations between SNPs and environmental variables and map environmentally-associated genomic turnover across the total study area. A total of 458 SNPs (28% of all SNPs) had R² values > 0 (average = 0.28) ( Figure 2A ). When testing model performance, the number of SNPs with R² values > 0 for all of the randomized datasets fell below the number observed for the real data ( Supplementary Figure S2 ) and the mean R² value generated for the real dataset fell within the upper 95% quartile of values generated for the randomized datasets ( Supplementary Figure S2 ), both indicating that the GF model shows a stronger association between environmental and genomic variation for our dataset relative to the set of randomized datasets. Precipitation of the coldest quarter, latitude, and precipitation seasonality were the most important environmental predictors of genomic turnover ( Figure 2B ). Projected associations between allele frequencies and these three predictor variables revealed areas of pronounced genomic turnover throughout the Cameroonian highlands (green to orange), forest-savanna ecotone of south-central Cameroon (orange to green), across the equator (green to blue), and from the coast to the interior of Gabon (purple to blue) ( Figure 2C ).

Predictions of environmentally-associated genomic turnover under future climate change projections showed similar patterns of genotype-environment associations across the landscape relative to current predictions. When we subtracted the current prediction of environmentally-associated genomic turnover from the future predictions under each of the four climate change projections (RCP 2.6 & RCP 8.5 for the years 2050 and 2080), we found a number of areas with high genomic vulnerability across the landscape ( Figure 2D ; Supplementary Figure S3 ). Areas with high genomic vulnerability (i.e. greater than 50% difference in environmentally-associated allele frequencies) occur to the north and south of the Cameroon highlands and throughout the southwest region of Gabon.