Genetic adaptation to urban living: molecular DNA analyses of wild boar populations in Budapest and surrounding area

Studies of wild boar, Sus scrofa Linnaeus 1758, in urban and suburban areas of Budapest, Hungary, have indicated that these populations do not have continuous contact. Based on the assumption that the city has a discrete population, we hypothesized that the urban wild boar would differ genetically from those in suburban areas. Analysis of single-nucleotide polymorphism (SNP) data using the GeneSeek Genomic Profiler (GGP) Porcine 50 K system (Neogen, Scotland, UK) differentiated three populations: Buda (B) from the Western bank of the Danube; Buda Surrounding (BS); and Valkó (V) from the Eastern bank of the Danube. The coefficient of genetic differentiation (FST) for the B and BS populations was low. The inbreeding coefficients of the populations BS and V were close to zero, while population B had a high positive value reflecting the influence of founders and the inbreeding of the continuous urban population. The genome regions that were most differentiated between the B and BS populations were analyzed based on the FST values of the SNP markers using a mixed linear multi-locus model and BayeScan software. The most differentiated marker, WU_10.2_18_56278226, was found on chromosome 18. The surrounding region contained several candidate genes that could play important roles in adaptations related to human-induced stress. Two of these, encoding the adenylate cyclase 1 (ADCY1) and inhibin beta A chain precursor (INHBA) genes, were sequenced. While IHBA gene did not display variation, the allele distribution of the ADCY1 gene in the B population was significantly different from that of the BS population supporting the parapatric differentiation of wild boar.

Wild boar (Sus scrofa) have frequently been reported to cause conflicts in European cities. These large, prolific mammals (Bieber and Ruf 2005) can cause property damage (Campbell and Long 2009;Saito et al. 2011;Schley et al. 2008) and provoke fear among human populations (Kotulski 1 3 and König 2008). Populations are found on most continents and have adapted to a wide variety of habitats (Baubet et al. 2004;Cahill et al. 2012;Cuevas et al. 2010;Gabor et al. 1999;McGaw and Mitcell 1998;Siemann et al. 2009;Thomson and Challies 1988;Yasin 2011). Conflict has resulted from the spatiotemporal activity of S. scrofa populations overlapping with that of humans (Cahill et al. 2012). The problems have ranged from general inconvenience to serious issues such as animal-vehicle collisions (Cahill et al. 2012;Licoppe et al. 2013;Thurfjell et al. 2015). Wild boar can be a threat to birds (Oja et al. 2017), destroy plants in gardens and parks (Schley et al. 2008;Licoppe et al. 2013;Hamrick et al. 2011), injure domestic and farm animals (West et al. 2009;Schley and Roper 2003), and damage property (Licoppe et al. 2013). Kotulski and König (2008) reported that the main issues with wild boar in Berlin occurred in parks and gardens, although they also caused human injuries (Stillfried et al. 2017;Castillo-Contreras et al. 2018). A related problem is the emergence of zoonotic diseases, which can be a threat to humans, and to domesticated and farm animals (Hygnstrom et al. 2011;Riley et al. 1998;Schley et al. 2008). The published literature has linked wild boar to Escherichia coli infection (Cahill et al. 2012;Licoppe et al. 2013), leptospirosis (Cahill et al. 2012;McGaw and Mitcell 1998;Meng et al. 2009), brucellosis (Hamrick et al. 2011;McGaw and Mitcell 1998;Meng et al. 2009;West et al. 2009), salmonellosis (Hamrick et al. 2011), rabies (Hamrick et al. 2011;McGaw and Mitcell 1998), classical swine fever (Hamrick et al. 2011;Meng et al. 2009), trichinosis (Hamrick et al. 2011;McGaw and Mitcell 1998;Meng et al. 2009;West et al. 2009), toxoplasmosis (Cahill et al. 2012;Meng et al. 2009), African swine fever (McGaw and Mitcell 1998;Meng et al. 2009), pseudorabies (Meng et al. 2009;West et al. 2009), Japanese encephalitis (McGaw and Mitcell 1998;Meng et al. 2009), andtuberculosis (McGaw andMitcell 1998;Meng et al. 2009). These diseases can be a serious health risk not only to other animals but also to humans. These findings indicate that interventions should be considered for managing the wild boar populations, the damage caused by boars, and reducing the human wild boar conflicts in urban areas.
The appropriate intervention will depend on whether the target species is transient or has been settled in the area. In the former case, an effective way to reduce conflict may be to intervene in suburban areas; in the latter case, conflict management should be resolved in part or in whole within the city. Our previous studies have revealed that urban wild boar in Budapest stayed in the city. 89.95% of data points of sows, tracked with a global positioning system transmitter, were within Buda territory during the 265-day study period (Csókás et al. 2020;Sütő et al. 2020). This was supported by the findings of our study of the urban feeding habits of wild boar, which showed that, although they mainly sought natural food sources within populated areas, their primary food sources in this environment differed from those in their natural habitats (Heltai et al. 2016;Katona et al. 2018a, b). In the 12th district of Budapest, Buda, where our study was carried out, the proportion of oak trees was high; therefore, the most advantageous basic food source for wild boar was acorns (Katona et al. 2018b). Based on the assumption that there was no continuous contact between the central and peripheral populations, we hypothesized that there was a breeding population within the city that was genetically different from the suburban population.
Molecular genetic analyses of wild boar have previously been useful in the investigation of their population genetic structure (Hampton et al. 2004), the relatedness of populations under anthropogenic effects (Nikolov et al. 2009), and their genomic diversity related to local pig breeds (Muñoz et al. 2019;Bâlteanu et al. 2019). We therefore chose to test our hypothesis via DNA analyses based on SNP-chip genotype data and sequencing of candidate genes around the top hit provided by genome-wide association analysis.

Samples
In total, 81 wild boar muscle samples were collected from the Buda (B) region (N = 28), the Buda Surrounding (BS) area (N = 26), and the Valkó (V) area (N = 27) during hunting events organized by the game-management units (Fig. 1). The V population served as an outgroup in our analyses. Based on the estimates of Pilis Park Forestry Company more than 10% of the populations have been sampled. Samples included all, hunter-picked animals. The ratio of adult:young animals were 2.11, 2.25, and 2.37 in B, BS, and V area, respectively. Approval from the animal research ethics committee was not necessary as the hunting events were organized independently from the research presented here. These hunting events were part of the routine wildlife-management activity based on the hunting law (Law No. LV. of 1996 on the Protection, Management and Hunting of Wildlife) and were inspected by the local government authority. We collected tissue samples from the carcasses after the hunting events. Samples were stored at − 20 °C. The B area is the western, hilly part of Budapest located west of the Danube and has a considerable amount of woodland (Fig. 1). The BS area refers to the forested region beyond the Buda boundary. The V area is located east of the Danube, around 40 km from Budapest. Evidence of the existence of the downtown population of S. scrofa prior to the first scientific report by Bogdán and Heltai (2014) was collected via Google using the keywords 'wild boar' and 'Budapest' (Online Resource 1).

Quality control
Data evaluation followed the methods described in our previous surveys (Bâlteanu et al. 2019;Zsolnai et al. 2020). DNA typing was performed using a GGP 50 K Porcine SNP chip (Neogen, Scotland, UK). Quality control procedures were conducted on the raw data using SNP and Variation Suite (SVS) software v.8.8.1 (Golden Helix, Bozeman, MT, USA). SNPs with a call rate < 95% or with a minor allele frequency (MAF) < 0.05 were eliminated from the dataset. Duplicated samples and/or closely related animals identified by an identical by state (IBS) value > 0.95 and individuals with a genotype call rate < 95% were also removed. After filtering steps, the final dataset included 77 animals (N B = 25, N BS = 25, and N V = 27) and 47,765 SNPs.

Diversity statistics and genome-wide screening
The genetic distance (F ST ) and the of inbreeding coefficients (F IS ) were calculated to describe the studied populations. The genetic distance of the markers (F ST_markers ) was calculated to look at the differences at genomic level. Determination of the previously mentioned F values was performed using SVS. The inbreeding coefficient was estimated as: (observed-expected number of homozygous markers)/(number of genotyped markers-expected number of homozygous markers). The expected number of homozygous markers was calculated as where L was the number of genotyped markers, p and q were the frequencies of the alleles, and T was twice the number of non-missing genotypes for marker j. F ST was estimated according to Weir and Cockerham (1984).
To display the relative positions of the sampled animals relative to each other, we used the PLINK software v1.9 (-mds-plot 2 and -cluster options) to build a multidimensional scaling plot using a genome-wide IBS pairwise distance matrix (Purcell et al. 2007). The IBS pairwise similarity between any two individuals was calculated as ((number of markers sharing two alleles + 0.5 * number of markers sharing one allele)/number of markers). Identities of frequency distributions of the pairwise IBS values of the studied populations were tested by Kruskal-Wallis probe (Kruskal and Wallis 1952).
Tests of the proportion of mixed ancestry and population structure were performed via the ADMIXTURE software v.1.3 using the -cv option (cross-validation procedure was performed fivefold) for determining the most probable cluster number (K) from the value of cross-validation error in each Ks (Alexander and Lange, 2011). The cross-validation procedure was performed five times, and the algorithm was terminated when the log-likelihood increased by less than 10 −4 between iterations. Before analysis, the alleles of the SNP loci were recoded to numerical values 1 and 2 by PLINK using the -recode12 switch as required by ADMIX-TURE. The file format was also changed to the binary type by the -make-bed option to facilitate faster computation by ADMIXTURE.
To gain insight into the possible differences of the internal structure of runs of homozygosity (ROH) of the populations, calculations were carried out with the PLINK v.1.9 software using the -homozyg command (Purcell et al. 2007). The minimum length of an ROH was set at 1 Mb to minimize the detection of spurious ROH. To make sure that ROH length was not affected by low SNP density, the minimum number of SNPs that constituted an ROH was set to 50 considering the following calculation method (Lencz et al. 2007): where n s was the number of SNPs per individual, n i was the number of individuals, α was the percentage of false-positive ROH (set to 0.05), and het was the mean SNP heterozygosity across all SNPs. The density of SNPs was set to one l = log e n s ×n i log e (1 − het) , SNP for each 100 kb and a maximum distance of 1000 kb was allowed between two consecutive SNPs. The scanning window contained 50 SNPs, and the maximum number of missing SNPs per window was set to five with an allowance for one heterozygous SNP. The ROH was classified based on physical length into four size categories: 1 to ≤ 5 Mb; 5 to ≤ 15 Mb; 15 to ≤ 30 Mb; and > 30 Mb (Bâlteanu et al. 2019;Zsolnai et al. 2020). For each ROH category, the mean of the ROH per breed was calculated by summing the lengths of all ROH in a given individual for each one of the categories under consideration and divided by the number of animals. The inbreeding coefficient derived from the ROH genomic coverage (F ROH ) was calculated by dividing the total ROH length per individual by total genome length (2444 Mb) for each individual.
To look for differentiated loci in the comparison of B and BS, the following approaches were tested: (1) the genetic distance of the markers (F ST_markers ) was calculated by the SVS program, (2) BayeScan software (Foll and Gaggiotti 2008) was used to determine the departure of each SNP loci from neutrality, where the total number of iterations was 100,000, and (3) genotypic association tests were performed using a multi-locus mixed-model (MLMM) algorithm implemented in SVS. A genomic kinship matrix was used in MLMM for the correction of population structure (Bâlteanu et al. 2019;Segura et al. 2012). The applied model was: y = Xβ + Zu + e, where y denoted the population (B and BS were recoded to one and zero) of the sampled animals, X was the matrix of fixed effects composed of SNPs, Z was the matrix of random animal effects, e was the residual effects, and β and u were vectors represent the coefficients of fixed and random effects, respectively. Nearby genes of the identified SNPs were identified via S. scrofa Assembly Build 10.2 database.

Sequencing
Target regions, identified by genome-wide screening as described above, adenylate cyclase 1 (ADCY1), and inhibin beta A chain precursor (INHBA) encoding genes were selected for sequencing. The polymerase chain reaction (PCR) thermal cycle of the reactions was 95 °C for 1 min followed by 30 cycles consisting of 95 °C for 15 s, T °C for 15 s, and 72 °C for 35 s, where T was 65, 62, and 59 in cases spanning exon 2 and exons 6-7-8 of the ADCY1 gene and exon 1 of the INHBA gene, respectively. PCRBIO HiFi Polymerase (PCR Biosystems, London, UK) was used in the reactions according to the manufacturer's instructions (PCR Biosystems, Wayne, PA, USA).
The raw sequence data were processed by bwa 0.7.17-r1188 (Li and Durbin 2009); variant calling was performed by samtools 1.9 and bcftools 1.9 using htslib 1.9 (Li 2011) following the instructions in the respective manuals. Briefly: raw sequences were aligned by bwa to the bwa-indexed reference sequence. Sorting and pileup were accomplished by samtools. View, stats, and call options of bcftools were used to create variant files carrying the detected sequence variations. After allele calls, the distributions of alleles among the populations were compared via the Chi-square test (IBM SPSS 25 Statistics Premium Campus Edition, US).

Population analysis
To supplement the first scientific literature (Bogdán and Heltai 2014) on the existence of the downtown population data were gained via Google search (Online Resource 1). The first reports on the appearance of wild boar were made in 1993, when 0.2% of the human population used the Internet. The first report on human-wild boar conflict dated back to 2001 (at which time there was 27.7% Internet penetration). Since 2004, the reports were distributed throughout the year. The first video report was created in 2010 (at which time there was 65% Internet penetration).
The inbreeding coefficients of the populations BS and V were close to zero (F IS = 0.020 and 0.017, respectively), while population B had a high positive value (F IS = 0.120) ( Table 2). The 95% confidence intervals (CIs) of the F ST values did not overlap. The coefficients of genetic differentiation (F ST ) for the pair B and BS were 0.015. For the populations BS and V, the F ST was 0.049, while the distance of population B from V was 0.087. The 95% CIs of the F ST values did not overlap (Table 3).
The characterization of the length and distribution of the ROH (Fig. 2) showed that the estimates for the B group differed from those in the BS and V populations. The means of the ROH (Fig. 3) in the first (0-5 Mb) and the third (15-30 Mb) classes were similar regarding the BS and V groups. In all ROH categories (0-5, 5-15, 15-30,  To investigate the genetic relationship between the sampled populations, we built a multidimensional scaling plot calculated with PLINK (Fig. 4). Eigenvalues of first and second components were above one. Third, fourth, and fifth components have eigenvalues 0.753, 0.653, and 0.610, respectively. The first component (eigenvalue = 2.151, explained 39.06% of total variance) separated population V from B and BS. Partial overlapping was observed between the B and BS groups. The second component (eigenvalue = 1.341, explained 24.34% of total variance) divided both the B and BS groups. In the ADMIXTURE analysis of B + BS population (Fig. 5A), the most probable cluster number was two (at the lowest cross-validation error, K = 2). When B + BS + V populations were examined the difference disappeared between B and BS, the most probable K value was also two (Fig. 5B).
The pairwise IBS values of individuals within each population were visualized (Fig. 6) to check their distribution. The identity of the frequency distribution of pairwise IBS values of the three populations was tested by the non-parametric test (Kruskal and Wallis 1952). The p values of the pairs were p BS-B = 0.000, p BS-V = 0.000, and p B-V = 0.123.

Genome-wide screening and sequencing
The genetic distances of the markers (F ST_markers ) were calculated between populations B and BS where the maximum F ST_markers value (0.405) was harbored by the SNP marker WU_10.2_18_56278226 on chromosome 18, in agreement with the result of the mixed linear multi-locus model (-log 10 p value = 20.7, false discovery rate = 7.9e−17). Selection signs (p BayeScan = 0.78-0.86) were detected by the BayeScan program on chromosomes 5 and 13, and a strong diversifying selection sign (p BayeScan = 0.96, alpha = 1.58) was detected on chromosome 18, on the same marker as mentioned above.
The ADCY1 and INHBA genes were selected for further investigation around WU_10.2_18_56278226. Successful amplifications were achieved spanning exons 2 and 6-7-8 of ADCY1 and on the first exon of INHBA (Table 1). Variation among the samples was not detected on the first exon of INHBA, while ADCY1 displayed several differences (Table 4). Chi-square analysis was performed on the allele distributions of the ADCY1 gene between the studied populations (Table 5). In relation of the population pairs B-V, B-BS, and B-V, the detected differences at position 17,850 and at positions ranging from 51,087 to 51,993 were significant at p < 0.05 level.

Discussion
The suburban populations, BS and V, were found to have been indigenous and continuous within their investigated territory with significant population increases (Csányi et al. 2019). The first publication about the downtown population, B, dated back to 2014 (Bogdán and Heltai 2014) and was followed by further studies (Heltai et al. 2016;Deutsch and Heltai 2017;Csókás et al. 2020). Internet posts about wild boar appearances in the downtown area, B, prior to 2012 suggested that the population has been established since 1993; however, Tari et al. (2016) reported their first observation to have been in 2000. While the Internet penetration was relatively constant, at < 4%, for many years, there was a noticeable increase in the number of reports, which jumped to 178, in 2004 and 2005. Subsequently, the observations increased nearly co-linearly with the Internet penetration value. In 2005, the S. scrofa population appeared to reach a saturated state in terms of numbers (Online Resource 1), which might explain the abovementioned trends. Over such a long period of time, the genome of the downtown population could have undergone changes that were driven by the founder effect, stochastic changes, or human-induced selection. Samplings of the studied populations were random, animals are considered as representative ones for their groups.
Looking at the genome of the animals sampled, the inbreeding coefficient was high (F IS = 0.120), indicating that individuals within group B were more closely related than would be expected. The inbreeding coefficient differed from those of the suburban populations BS and V ( Table 2). The suburban population BS, which lived in close vicinity to the downtown group, and V, which lived on the other side of the Danube and was separated from B and BS by another barrier, the capital, did not differ in their inbreeding coefficients, indicating that their populations had similar genetic statuses.
The coefficient of genetic differentiation, F ST , differed for each studied pair (Table 3). B and BS displayed the lowest value (F ST = 0.015). Genetic differentiation was also low for the suburban populations BS and V (F ST = 0.049). The highest value, reflecting moderate genetic differentiation, was found between the B and V populations (F ST = 0.087). The The relative positions of animals based on their SNP loci were plotted (Fig. 4). The control or outgroup (V) was separated from populations B and BS. The partial overlapping of B and BS might reflect occasional territory swapping of the animals between regions B and BS with possible mating events in the other population. Both the B and BS groups were separated into two sets and it was assumed that these had different founders with different genetic backgrounds.
The lengths of the ROHs (Fig. 2) in the studied animals revealed long segments (400-800 Mb) that were seen almost exclusively in population B, in accordance with its higher inbreeding value F IS ( Table 2). The animals in populations BS and V had similar ROH lengths; indeed, according to the comparative analysis of the animals from the three regions (Fig. 3) in which the ROH categories and their 95% CIs were determined, the BS and V populations did not differ in the low category (0-5 Mb) and the third category (15-30 Mb). The downtown population, B, displayed the highest mean values in each ROH category. The difference in B versus BS or V and the similarity of BS and V were also confirmed by their genomic coverage values of runs of homozygosity (F ROH ). The F ROH_B value of population B was twice as high (0.165) as that of BS or V, which also indicated an inbred downtown population.
The results of the ADMIXTURE analysis of B and BS (Fig. 5A) provided at K = 2 revealed seven animals from population BS (similar to the MDS analysis in Fig. 4) that were sampled in the territory of the downtown population. We assumed that they were explorers from the suburban population, BS, as the same behavior has been observed for other species in which there are both transient and resident city populations at the same time (Iossa et al. 2010). Opposite pattern also occurs with individuals assigned to the B genetic cluster sampled in the territory of BS subpopulation. Distinction of B and BS populations disappears in higher context, when V population was included into the examined set (Fig. 5B), indicating that the separation of the populations B and BS is in the beginning phase. Nevertheless, it is questionable if separation will be stronger as time goes by or remains as reported here, since we can notice individuals in B and BS populations (Figs. 4, Fig. 5A) which are not clustered into their original sampling populations. This phenomenon indicates that members of urban and suburban populations occasionally swap territories.
The pairwise IBS frequency distribution values (Fig. 6) differed significantly regarding population pairs B-BS and BS-V (p BS-B = 0.000; p BS-V = 0.000). The frequency distribution of B and V populations did not contrast according to the test (p B-V = 0.123); however, visual inspection revealed that B had a bimodal frequency distribution, which differed from population V. The observed bimodal distribution of the downtown population suggested that a part of the group Vertical axis is the number of occurrences within a given category was more inbred. Bimodal distribution is expected when individuals with different genetic backgrounds exist in the same population. As reported, individuals from Buda spent 89.95% of their time in the urban environment (Bogdán and Heltai 2014) owing to the enhanced availability of their natural food, maize, offered at baiting sites and the managed lawns (Katona et al. 2018a, b). The wild boar of Buda live in close proximity to humans, which is expected to alter their behavior. Such behavioral changes can be learned (Changizi et al. 2002) or can be acquired via altered gene function (Jorgensen et al. 2019).
To locate the region of the genome in which the downtown animals differed most greatly from that the suburban wild boar, several methods were utilized between the B and BS populations. Three different algorithms exposed the most differentiated marker to be on chromosome 18. The locus WU_10.2_18_56278226 became the center of our search for candidate genes. The selection of candidate genes around this locus was based on the working hypothesis that downtown animals differed in personality from suburban individuals and this would be manifested in genes involved in brain functioning.
The nearby genes of WU_10.2_18_56278226 were identified on the S. scrofa genome assembly 10.2. Some of the genes around this marker have proven or postulated roles in the brain and neural processes, such as phosphoglycerate mutase 2 (PGAM2), drebrin like (DBNL), ADCY1, calcium/ calmodulin-dependent protein kinase II beta (CAMK2B), and INHBA. PGAM2 can be found in different tissues and forms: the slow-migrating muscle (MM) isozyme and fast-migrating brain (BB) isozyme (Jorgensen et al. 2019) are related to meat and production traits in domestic pigs (Fontanesi et al. 2008). DBNL plays a part in neuron morphogenesis (Inoue et al. 2019). ADCY1 is primarily expressed in the brain (https:// www. prote inatl as. org/ ENSG0 00001 64742-ADCY1/ tissue), plays a part in the memory process (Wang et al. 2004) and brain development (Sethna et al. 2017), and is associated with deafness in humans (Santos-Cortez et al. 2014). ADCY1 may affect the feeding efficiency of pigs via the cyclic adenosine monophosphate (cAMP) signaling pathway (Xu et al. 2018); its expression was found to be altered in a comparison of thyroid glands of Jinhua and Yorkshire pig breeds (Shen et al. 2016). CAMK2B is   (Küry et al. 2017). INHBA is associated with follicle and testis size, and plays a role in eye, tooth, and testis development (https:// www. ncbi. nlm. nih. gov/ gene? cmd= Retri eve& dopt= full_ repor t& list_ uids= 3624), and axon regeneration (Omura et al. 2015). ADCY1, which was one of our candidate genes selected for sequencing, displayed variations among the studied populations (Table 4). Based on Chi-square analysis of alleles of the ADCY1 gene, the downtown population B displayed a significant difference from the other populations BS and V (Table 5). Among the single-nucleotide variations in the ADCY1 region, four of the mutations were in exons, two of which were missense mutations G/A and T/C causing amino acid changes Alanine/Threonine and Tryptophan/Arginine on exon 6 and 7, respectively. Alanine is a strong helix former, while Threonine is a strong β-sheet former (Chou and Fasman 1974). Alanine/Threonine substitutions usually are on the surface of the 3D structures facilitating interactions, aggregations of proteins known in human (Podoly et al. 2010). Arginine and Tryptophan highly differ, e.g., in their hydrophobicity. The latter also plays important role in glycan-protein interactions (Yamaguchi et al. 1999). Mutations described above could alter structure of ACDY1 protein which could affect phenotypic traits or developmental processes (Xu et al. 2018;Sethna et al. 2017).

Conclusions
SNP-chip analysis revealed genetic differences between downtown and suburban wild boar populations. Since the urban environment requires altered behavior from the animals, we hypothesized that the downtown B population might carry well-differentiated markers in its genome. Such signs were revealed by different algorithms. The genomic site identified was found on chromosome 18. This region carries several genes connected with neural processes which was in line with our working hypothesis that downtown animals differ in personality from suburban individuals and such differences should be manifested in genes affecting the brain. Out of the candidate genes, INHBA and ADCY1, the latter displayed variations.
Sequence analysis revealed alterations in the investigated groups. PCR products, spanning exons 2, 6, 7, and 8 of the ADCY1 gene, carried mutations. Analysis of the allele distribution of the studied populations showed a significant difference in population B compared to populations BS and V. Further investigation of the structure and expression of the ADCY1 gene might yield further insights (Reed et al. 2020) into the manifestation of stress tolerance in wild boar living in the urban environment and could also confirm the parapatric differentiation of wild boar, S. scrofa Linnaeus 1758, as a consequence of partial isolation (Gavrilets et al. 2000).