Isolation by environment (IBE) is a wide spread phenomenon in nature. It is commonly expected that the degree of differences among environments is proportional to the level of divergence between populations in these environments. Consequentially, it is assumed that species’ genetic diversity displays pattern of IBE in the presence of a strong environmental cline if geneflow does not mitigate isolation. We tested this common assumption by analyzing the genetic diversity and demographic history of Pisum fulvum. P. fulvum inhabits very contrasting habitats in the southern Levant and is expected to display only minor migration rates between populations what makes it an ideal test case. Ecogeographic and subpopulation structure was analyzed and compared. Correlation of genetic with environmental distances was calculated to test the effect of IBD and IBE and detect the main drivers of these effects. Historic effective population size was estimated using stairwayplots. Limited overlap of ecogeographic and genetic clustering was observed, and correlation of genetic with environmental distances was statistically significant yet small. We detected a sharp decline of effective population size during the last glacial period. The low degree of IBE may be the result of genetic drift due to the past bottleneck. Our findings contradict the expectation that strong environmental clines cause IBE in the absence of extensive geneflow.

Plant material and habitat characteristics
Naim-Feil et al. (2017) collected 135 P. fulvum accessions spanning its documented ecological amplitude in Israel and the Palestinian Territories. Our collection was comprised of 126 of these accessions (Supplementary Table. S1). Within our study area, comprising Israel and the Palestinian Territories, it occupies a wide range of different habitats. In these habitats average temperature in the main growing season during winter ranges from 4.0°C to 14.3°C. The annual precipitation, which occurs almost entirely during the winter months, spans a range from 295 mm to 914 mm. The echo system of each sampling site, i.e., the climate conditions, soil content, and the vegetative vicinity is tabulated in Supplementary Table S1. For more detailed description of P. fulvum habitats see Abbo et al. (2013). We are unaware of published census size data for wild pea of which is rather difficult to assess. The number of individuals in a population strongly vary over years, and often in certain locations, no plants can be found in some years. Yet, P. fulvum, like most Mediterranean legumes, produces seeds exhibiting strong physical dormancy, and wild populations build considerable soil seed banks (Abbo et al., 2011; Berger, Shrestha, & Ludwig, 2017; Smýkal, Vernoud, Blair, Soukup, & Thompson, 2014). These seed banks allow repopulation after years with low germination rates. Flowers of P. fulvum are usually cleistogamous. Yet, flowers can open in later stages and they have been documented to potentially be prolific after opening (Bogdanova & Berdnikov, 2000). Although Pisum sp. are mostly considered as selfing species, the degree of self-fertilization is debated (Polowick et al., 2002, Smýkal et al., 2018). Due to its selfing nature and its big seeds is can be expected that long range geneflow in P. fulvum is unlikely and populations are well diverged, similar to the situation in Pisum sativum ssp. elatius, the only other wild pea species (Smýkal et al., 2018).

Genotyping by sequencing
DNA samples from individuals of the association panel were digested by MspI/PstI restriction enzymes. Genotyping-by-sequencing (GBS) libraries were prepared by ligating the digested DNA to unique nucleotide adapters (barcodes) using 48-plex, followed by PCR amplification. Libraries were placed in an Illumina NextSeq 500 instrument using 75-base single-end sequencing at the LGC Company. Preprocessing was performed by the LGC Company who generated the GBS sequence reads as follows. All library groups were demultiplexed using the Illumina bcl2fastq (v2.17.1.14) software. No mismatches or Ns were allowed in the inline barcodes, but Ns were allowed in the restriction site. Sequencing adapter remnants were clipped from all raw reads. Reads with final length <20 bases were discarded and reads with 5’ ends not matching the restriction enzyme sites were filtered out. The high quality reads were aligned to the reference genome of Pisum sativum (Kreplak et al., 2019) using the Burrow-Wheelers Alignment tool (Version 0.5.9; Li & Durbin, 2009) with default alignment parameters, the alignment data was processed with SAMtools (version 0.1.19; Li et al., 2009). Sequence variants were identified using freebayes (version 1.1.0-50-g61527c5; Garrison & Marth, 2012). The combined read depth of 10 was used across samples for identifying an alternative allele as a variant, with the minimum mapping and base quality filters of 1 and 3, respectively. Subsequent filtering steps included only biallelic SNPs with genotype calls across the population ≥ 50% and a minor allele frequency (MAF) ≥ 0.01. The QUAL estimate was used for estimating the phred-scaled probability. Sites with a QUAL value less than 20 were removed. The filters were applied using VCFtools (Danecek et al., 2011) and in-house Perl scripts. Afterwards, missing genotype calls were imputed with Beagle (version 5.0; Browning, Zhou, & Browning, 2018; Browning & Browning, 2007) using 30 iterations, 10 burn-in steps and an effective population size of 50,000 because P. fulvum presents rather small population sizes with high inbreeding rates, so the default value was not appropriate. We used individual’s heterozygosity rates to test for potential sample contamination or multilocus mapping. A deviation of more than three standard deviations from the mean heterozygosity rates were considered as too high. All accessions passed this criterion. Linkage disequilibrium (LD) pruning was conducted with PLINK (version 1.90b5.4; Purcell et al., 2007). LD was estimated by r2 in a sliding windows approach. 50 kb windows were shifted by 5 variants in each step and SNPs within one window with r2 ≥ 0.2 were dropped. After application of all filters, we obtained 5144 SNPs in total and 3951 LD pruned SNPs.