A Study of Hybridization Between Marmota baibacina and M. sibirica in Their Secondary Contact Zone in Mongolian Altai

The role of hybridization as one of the factors of speciation in mammals has been underestimated for a long time, but now there is a lot of data on its impact in mammalian evolution. Hybridization of species often occurs in their secondary contact zones, which is a natural model for testing factors that ensure species integrity. Studies of hybrid zones are increasingly revealing the essential role of ecological and behavioral features both in initiating crossbreeding and in maintaining interspecific barriers. We studied the hybridization of two species of marmots Marmota baibacina and M. sibirica in the zone of sympatry in Mongolian Altai Mountains. We used a bioacoustic approach to determine the localization of individuals of different species and their cohabitation sites. Genetic typing with two diploid nuclear markers and one marker each of paternal and maternal lines was used to identify hybrids. Habitat preferences of marmots were studied to understand the conditions for the formation of heterospecific pairs. We found a high proportion of hybrid individuals in boulder screes where conditions for the formation of heterospecific pairs probably exist. Our data indicate the viability and fertility of F1 hybrids and their descendants. We hypothesize that the environmental preferences and behavioral features of both species of marmots are important factors that both create conditions for hybridization and limit hybrid dispersal.


INTRODUCTION
The problem of defining species boundaries, their permeability, and estimating the impact of their disturbance for species delimitation is addressed according to the species concept that is a source of long-term discussions in evolutionary biology (Harrison and Larson, 2014;Zachos, 2016). Within the genetic species concept, maintaining genome integrity is one of the main species criteria (Wu, 2001). The genomic effects of interbreeding of heterospecific individuals are diverse and depend on many factors. On the one hand, interspecific hybridization and the breaking down of interspecific barriers can lead to the blurring of species boundaries and, in some cases, to the extinction of one of the species (Todesco et al., 2016). On the other hand, hybridization may enrich genomes with new variations, multiplies local adaptations, and in some cases, may initiate speciation (Abbott et al., 2013;Eberlein et al., 2019). But also hybridization may not have a significant impact on species genomes. The hybrid zones formed during secondary contact of previously spatially separated species may exist for quite a long time, appear and disappear periodically, and vary in their location depending on environmental changes (Wielstra, 2019). Formation of robust postzygotic barriers and gametic isolation during independent evolution requires long-term spatial isolation, disruption of which by secondary contacts under some conditions may lead to gene introgression (Abbott et al., 2013).
Advanced application of genetic approaches in studying gene flow and the consequences of species genome mixture makes it possible to detect population processes that do not appear on the phenotypic level or have occurred in the past (Chavez et al., 2011;Bono et al., 2018;Fadakar et al., 2020). At the same time, studies of phenotypic variability in hybrid zones allow us to understand the contribution of admixture to the adaptive capacity of a population (Fontsere et al., 2019).
Hybridization is more common in plant species (25% species involved) than in animals (10%) (Mallet, 2005). For mammals, hybridization has been previously recognized as extremely rare (Stebbins, 1959), but intensive researches using genetic approach have changed this view (Mallet et al., 2016). Current estimates suggest that at least 6% of mammal species are involved in interspecific hybridization to different extents (Mallet, 2005). This may vary across mammalian taxa, but numerous examples indicate that gene flow is a common phenomenon in mammals. The evolutionary role of hybridization in mammals has not yet been sufficiently understood. Most of the studied cases of hybridization in mammals are associated with secondary contacts of allopatric species that arose under the influence of cyclic glacial events (Colella et al., 2018). Revealing ancient hybridization in different taxa also broadens understanding of the contribution of hybridization to mammalian evolution (e.g., Marques et al., 2017). In the ground squirrels Marmotini up to a quarter of species are involved in hybridization events of varying intensity, from sporadic hybridization to the formation of hybrid zones and wide introgression according to our estimates of the literature data (Semenova, 1967;Nikol'skii et al., 1983;Smirin et al., 1985;Good et al., 2008;Stangl et al., 2012;Thompson et al., 2013;Ermakov et al., 2015;Frare et al., 2017;Ivanova et al., 2017;Leitner et al., 2017;Kapustina et al., 2018).
Marmots (Marmota Blumenbach, 1779) are large rodents distributed in North America and Eurasia, inhabiting open planes of mountain and lowland landscapes. Like most of the ground squirrels, they are burrowing animals. The distinctive features of marmots are the subdivision of the colony into family groups characterized by stability and territoriality, and the social behavior of most species (Armitage, 2003). Stable marmot family groups consist of a breeding couple and their offspring of different ages (juvenile and subadults) (Bibikow, 1996;Armitage, 2014). Hybridization was found only in Palearctic marmots. The crossing of species in all known cases occurs with the involvement of the gray marmot (Marmota baibacina Kastschenko, 1899) in zones of sympatry with the steppe marmot (M. bobak Müller, 1776), the long-tailed marmot (M. caudata Geoffroy, 1844), and the tarbagan (M. sibirica Radde, 1862; Semenova, 1967;Nikol'skii et al., 1983;Smirin et al., 1985).
To date, not a single hybrid zone of marmots has been studied comprehensively. One of the methods used in preliminary studies was bioacoustic. The advantage of acoustic signals as a model of variability is that they are amenable to visualization and quantification of their parameters using digital technology. Sound signals are widely used in the study of hybrid zones in birds and mammals (den Hartog et al., 2008;Wyman et al., 2011;Campbell et al., 2019). Most often, however, such studies focus on the role of specific signals in ensuring assortative mating by recognizing conspecific mates. Vocal traits can provide important information when studying hybridization in taxa with genetically determined acoustic signals not involved in mating behavior. Alarm calls typical for most ground squirrel species have exactly these properties. Species specificity and genetic determination of alarm calls of marmots and ground squirrels have been shown earlier (Nikol'skii, 1976(Nikol'skii, , 1979Blumstein, 2007). Studies of hybrid zones of four Spermophilus species revealed the formation of peculiar alarm calls in hybrids, which was confirmed by genetic data Nikol'skii and Starikov, 1997;Kuzmin et al., 2011;Titov et al., 2020). Characters of alarm calls were used in M. baibacina-M. bobak secondary contact zone for hybrid identification (Nikol'skii et al., 1983). We initiated the study of the hybrid zone of M. baibacina-M. sibirica to confirm hybridization in marmots by genetic method in combination with bioacoustic and ecological approaches.
Gray or Altai marmot (M. baibacina) and Tarbagan or Mongolian marmot (M. sibirica) are closely related, relatively recent divergent species which have inherited significant morphological similarity (Kryštufek and Vohralík, 2013) despite appreciable genetic differences (Brandler et al., 2010a;Steppan et al., 2011). The chromosomal sets of both species have similar parameters (2n = 38, NF = 70) and minimal differences (Brandler et al., 2008). Both species are quite ecologically supple and inhabit a wide range of mountain and intermountain plains steppe outside the zone of sympatry. They have a social organization typical for marmots, based on family groups which consist of a breeding couple, juvenile and subadult individuals (Bibikow, 1996;Adiya, 2007). Pairs in these species form during the summer season long before reproduction time and are maintained for several years. In Mongolia, copulation in both species occurs in spring at the end of hibernation before the first exit to the surface from burrows (Adiya, 2007).
The secondary contact zone of M. baibacina and M. sibirica in the Mongolian Altai is an example of the limited sympatry of two widely distributed species of territorially conserved social mammals. Morphological hybrids were first discovered within this zone in the upper Ulagchin-Gol river valley during the study of the distribution boundaries of Mongolian marmots in the 1980s (Smirin et al., 1985). The first studies outlined the boundaries of a mixed population, described the peculiarities of fur coloration and craniometric characteristics of the supposed hybrids. Variability in the alarm call of marmots was studied as one of this species characteristics (Formozov and Nikol'skii, 1986). A single signal presumably belonging to a hybrid was identified among species-specific signals. At the same time, a difference in the habitat preferences of M. baibacina and M. sibirica was identified (Smirin et al., 1985;Rogovin, 1992). There is information about the existence of several more mixed settlements in the zone of overlapping areas of M. baibacina and M. sibirica, which to date remain unexplored.
The mixed population of M. baibacina and M. sibirica in the valley of the Ulagchin-Gol river is located within the limited area of a small mountain valley and is relatively isolated from the nearest marmot populations. It is a suitable model for studying the various factors affecting the hybridization process in the natural population. We aimed to study the genetic variability of this population in order to provide genetic evidence for hybridization and estimate the real ratio of parental species and hybrids. Since both species have species-specific vocalizations, we used recorded vocalizations as a potential marker of hybridization. Finally, both, genetic and acoustic data, were used to assess the spatial distribution and habitat preferences of each species in the study zone.

MATERIALS AND METHODS
The mixed population in the zone of sympatry of M. sibirica and M. baibacina was studied in Ulagchin-Gol river valley, which is a tributary of the Bulgan river in the central part of Mongolian Altai. The place of study is situated in Bayan-Ölgii Aimag in 180 km south-southeast from the Ölgii city with geographical coordinates of the central point at 47 • 27 N 90 • 53 E and absolute altitudes from 2500 to 3300 m above sea level (m a.s.l.) (Figure 1). The fieldwork was carried out by the team of the Joint Russian-Mongolian Complex Biological Expedition of the Russian Academy of Sciences and Mongolian Academy of Sciences during the summer seasons (June-July) 2007-2009. In our study, we used bioacoustic, genetic and ecological methods.

Bioacoustics Collecting and Analysis
An acoustic signal from marmots was used for remote species diagnosis of individuals in a mixed population. The main task of sound signal analysis was to study the spatial locations of marmots of different species in the population area and to identify places of immediate contact of heterospecific individuals. The alarm call of marmots has well-defined genetically inherited species-specific characteristics (Nikol'skii, 1976(Nikol'skii, , 1984. It has been shown earlier (Brandler et al., 2010b) that the use of the alarm calls as a diagnostic feature is a convenient, non-invasive method to study the spatial distribution of marmots. During summer molting, which lasts from the end of May to mid-July on average, it is rather difficult to distinguish M. sibirica and M. baibacina in a mixed population by their color and structure of fur. Vocalizations of M. sibirica and M. baibacina, however, are well recognizable by hearing (Supplementary Audios 1, 2), which allows estimating the spatial distribution of individuals of these species during field surveys.
The marmots' alarm calls were recorded in the field in 2008-2009 using a digital recorder "Marantz-PMD 660" ("Marantz, " Japan) and a condenser microphone "Audio-Technica AT897" ("Audio Technica, " Japan) in uncompressed digital format Wave (WAV) with a discretization frequency of 48 kHz. The acoustic response was provoked by the researcher standing at some distance from the animal or moving slowly toward it. The coordinates of the calling animal locations were determined using a GPS navigator. GPS readings were obtained either directly on the burrow in which the calling marmot was hiding, or in case the burrow was located in a hard-to-reach place, by measuring the distance from the researcher to the animal using a laser rangefinder and azimuth counting by a compass, with the later calculation of coordinates by mapping these data on the space image in the program OziExplorer v3.95 (D&L Software Pty., Ltd.) using the researcher position GPS data. Sounds were also recorded during tissue sampling when marmots were in cages (SM1). In total, the records of 148 individuals from the contact zone of M. baibacina and The acoustic characteristics of the alarm calls were analyzed visually using graphical images of oscillograms and spectrograms in the SpectraLab 4.32.17 software ("Pioneer Hill Software LLC, " United States) and Raven Pro 1.6.1 software (Center for Conservation Bioacoustics, 2019).
The earlier defined temporal and frequency characteristics of M. baibacina and M. sibirica alarm calls were used as species-specific features (Nikol'skii, 1976(Nikol'skii, , 1984. The alarm call of M. baibacina differs from M. sibirica by the presence of a pause between the lowfrequency (LF) and high-frequency (HF) components (Figure 2 and Supplementary Figures 1, 2). The characteristic of the HF frequency modulation of M. sibirica alarm call is usually symmetrical, while in a M. baibacina, the final frequency of this element is significantly higher than the initial one. In addition, in most cases, the spectrum of an alarm call of a M. sibirica has sidebands above and below the fundamental frequency caused by the amplitude modulation of the base oscillation (Nikol'skii, 2007). Species identification of an alarm call of individual was identified based on the analysis of these features.
In addition to digitally recording alarm calls for acoustic analyses, we opportunistically recorded the location of individuals of both species when animals vocalized as researchers walked along the river in 2008 (Figure 3). In these cases, we assigned each call to a species based on species-specific differences (Supplementary Audios 1, 2). The recognition distance of species alarm calls was about 300 m. Coordinates of the beginning and the end of this route were recorded on a GPS unit.
Tissues Sampling, DNA Extraction, Amplification, Sequencing, and Restriction Analyses   Table 1). The information on the field species identification of animals by features of hair color, sex, age, GPS coordinates of the capture location was collected for each specimen. Due to difficulties in visual species identification during molting, species identification was made only for individuals with the least molting fur. The others were marked as M. sp. The tissues of 14 specimens of M. sibirica and 11 specimens of M. baibacina from "pure" (or "reference") populations located at distances from 150 to 1,300 km from the Ulagchin-Gol valley outside the zone of sympatry (Supplementary Table 1) kept in the same collection were examined to determine the species specificity of genetic markers.
DNA extraction from tissue samples was carried out using the standard salt method (Aljanabi and Martinez, 1997) or commercial kits for DNA extraction Diatom DNA Prep ("Isogen, " Russia). The DNA precipitate was dissolved in deionized distilled water or a standard tris-EDTA buffer (TE) pH 8.0.
Genetic diagnosis of all specimens was carried out using nDNA markers: i5HoxB-5th intron of HoxB (homeobox containing gene); i13BCR-13th intron of BCR gene (breakpoint cluster region); i8SmcY-fragment of the 8th intron malespecific histocompatibility antigen, as well as mtDNA fragmentcytb (cytochrome b gene). Nucleotide sequences obtained from individuals originating from "pure" populations were used as references to identify the species specificity of each marker (Supplementary Table 1). Additionally, cytb sequences from GeneBank were used (M. sibirica AF143937, AF143938, and M. baibacina AF143915). Species-specific substitutions forming the target sites for specific endonucleases were found in three markers (i5HoxB, i13BCR, cytb). This allowed us to carry out genetic tests of the total mixed population sample using restriction analysis of PCR products of these fragments. The species specificity of i8SmcY was identified by nucleotide sequence analysis for all males.
Control region of mtDNA (CR) nucleotide sequences were studied in 51 individuals from a mixed population (Table 1) in order to evaluate a genetic variability in population samples of both species. Sequences of 3 M. sibirica and 3 M. baibacina individuals from "reference" populations were included in the analysis to identify the species specificity of CR haplotypes (Supplementary Table 1). Table 2) were used for the amplification of the investigated DNA fragments. Amplification reactions were carried out either in 25 µl of the reaction mixture using "DNA amplification" kit ("Silex M, " Russia) covered with mineral oil in thermocycler "Tertsik" ("DNA-Technology, " Russia) under the following conditions: initial denaturation 94 • C 3 min; 30 cycles-denaturation 94 • C 1 min, annealing (temperature in Supplementary Table 2) 1 min, extension 72 • C 1 min; final extension 72 • C 10 min, or in 20 µl reaction mixture Screen Mix ("Eurogen, " Russia) without oil in Veriti Termo Cycler ("Applied Biosystems, " United States) under the following conditions: initial denaturation 95 • C 3 min; 30 cycles-95 • C 20 s, annealing 40 s, 72 • C 40 s (or 1 min for cytb and CR); final extension 72 • C 7 min (or 10 min for cytb and CR).

Specific primers (Supplementary
Species identification of the i5HoxB alleles was carried out using MspI\HpaII restriction enzyme ("Fermentas, " United States). The i13BCR and cytb alleles were tested using PmlI restriction enzyme ("Fermentas, " United States). Enzymatic reactions were carried out in 10 µl of reaction mixture containing 20-40 ng of the PCR product, reaction buffer and 2-3 units of the restriction enzyme (according to manufacturer's protocol). The samples were incubated at 37 • C for either 1.5-3 h or overnight. The resulting DNA fragments were separated by electrophoresis in 2% agarose gel containing ethidium bromide with subsequent visualization using "Bio-Rad ChemiDoc MP System" ("Bio-Rad, " United States).

Molecular Data Analyses
DNA sequences were aligned using ClustalW and edited manually where necessary in MEGA X software (Kumar et al., 2018). The detection of species-specific substitutions in the i5HoxB, i13BCR, cytb sequences, as well as the genetic typing of i8SmcY haplotypes, were carried out visually in aligned nucleotide sequences. The median-joining method (Bandelt et al., 1999) was used to construct a CR haplotype network using Network 10 software ("Fluxus Technology Ltd." at) 1 .
Standard population genetic characteristics were calculated to estimate genetic diversity. Taking into account the high proportion of hybrids in the mixed population (see section "Results") reflecting the high potential of heterospecifics for interbreeding, we made the assumption that the population of two cohabiting marmot species of the Ulagchin-Gol valley is a single panmictic population. Frequencies of species-specific alleles (sibirica and baibacina) for each locus, heterozygosity expected (He) and observed (Ho) for each diploid locus and means for samples for all diploid loci, gene diversity, and genetic distances, Hardy-Weinberg equilibria tests were calculated in Arlequin 3.5.2.2 software (Excoffier and Lischer, 2010). The Wright fixation index (F is ) was calculated in the FSTAT 2.9.3.2 software (Goudet, 2001).

Spatial Analysis and Habitat Assessment
Spatial distribution of the studied marmots in the mixed population was analyzed using satellite images with a resolution of up to 30 m/pixel freely available by Google (United States) and digitized topographic maps at scales of 1:200,000 and 1:100,000, on which points of location of animals with identified genotypes, as well as vocalizing animals, were mapped in accordance with the recorded GPS coordinates in the OziExplorer v3.95 ("D&L Software Pty., Ltd.") and MapInfo 9.5 ("Pitney Bowes Software Inc., " United States) software.
To assess the habitat preferences of marmots with different genotypes, the locations of all genotyped individuals were classified as one of three major habitat types: (1) valley bottom and lower slope, (2) upper slope and watersheds, (3) large stone plots (boulder screes). An analysis of the population density of marmots in each of these habitats and the distribution of both species of marmots and hybrids was carried out. The significance of the association between the allocation of marmots with different genotypes and habitat types was examined using Fisher's exact test in STATISTICA 8.0 software ("StatSoft, Inc."). The p-values from the test were computed for 2 × 2 contingency tables with pairwise grouping of the evaluated factors. We applied the two-tailed Fisher's exact test for all tables except those that included cells with 0, for which the one-tailed version was used. Altitude above sea level was recorded for each animal that was genotyped to understand the impact of altitude on the spatial distribution of marmots.  It is known that marmot population densities are highest in high-quality habitats where environmental needs of individuals are met (Bibikow, 1996). Based on this, we used population density by the number of family groups per unit area (family/km 2 ) to estimate the suitability of the above three major habitat types for marmots. The family group sites were determined based on the relative location of permanent and/or winter burrows following the common approach (Mashkin and Chelintsev, 1987). The specific placement of visually different types of burrows (permanent, winter, temporary) allowed an accurate identification of the family sites' area frontiers. The family sites were accounted totally within the surveyed area in the upper part of the Ulagchin-Gol valley. The area of the surveyed territory and areas with different habitat types was calculated by mapping the GPS data to a digitized physical map of 1:100,000 scale.
All experimental protocols were approved post hoc by the recently formed Ethics Committee for Animal Research of the Koltzov Institute of Developmental Biology RAS (Approval No. 36, 05\03\2020) in accordance with the Recommendations for Laboratory Practice in Russian Federation. Animals were treated according to established international protocols, such as Guidelines for Humane Endpoints for Animals Used in Biomedical Research and Guide for the Care and Use of Laboratory Animals. The study is classified as one of moderate severity and the number of the animals used in this study does not exceed the number needed to meet statistical significance.

Bioacoustics Data
Our data suggested that frequency-temporal features of M. sibirica and M. baibacina alarm calls from reference populations outside the zone of sympatry correspond to the specific features described earlier for these species (Nikol'skii, 1976(Nikol'skii, , 1984. All signals of M. sibirica had no pause between the low and high-frequency components and had a symmetrical high-frequency component with sidebands induced by amplitude modulation (Figure 2A and Supplementary Figures 1A-D). All signals of M. baibacina from all reference populations had a pause between LF and HF and an asymmetrical HF frequency modulation ( Figure 2D and Supplementary Figures 2A-D). No amplitude modulation was found in the last signals. Oscillograms and spectrograms of these alarm calls were used as reference standards for species-specific diagnostics of acoustic signals of marmots from mixed populations.
In total, alarm calls of 148 individuals of marmots were analyzed in the study area; 77 of them were identified as M. sibirica and 71 as M. baibacina. Among the analyzed signals, 5 were recorded from animals in cages. Three of them were M. sibirica and two were M. baibacina based on a genetic test. Sound signals of other captured and genotyped marmots could not be recorded due to the short time of keeping in cages and the absence of their vocalization at this time. All the analyzed individuals had clear corresponding species-specific features of the acoustic signal described above (Figure 2 and Supplementary  Figures 1, 2). The spatial distribution of these individuals is shown in Figure 3. The records of M. sibirica alarm calls were collected mostly in the northern upper part of the valley right up to the headwaters of the Ulagchin-Gol river, where there are more gentle slopes with rich alpine meadow grass ( Figure 3A). In contrast, the marmots producing calls typical for M. baibacina prevailed in the southern lower part of the valley (Figures 3B-D). During the survey of the middle part of the valley by walking along the Ulagchin-Gol riverbed, we heard only M. sibirica alarm calls (21 notices) on the valley bottom (yellow dotted line in Figure 3). We have identified in the upper part of the valley the territory with predominance of boulder screes where both types of acoustic signals were equally common (framed in Figure 3). 17 sibirica and 25 baibacina phonotypes were recorded in this place.
An aberration of HF modulation was found in 12 individuals (17%) with baibacina phonotype (Figure 2F). In contrast to the normal form, the frequency of HF of the aberrant signal initially increased with normal rate, then the frequency rose sharply in a short time interval, followed by a weakly modulated fragment of the signal, including the maximum of the fundamental frequency. The signal ended with a weak frequency decrease similar to the normal alarm call. Only one individual had an alarm call that was different from the normal ones among all M. sibirica signals (Figure 2C), which is a rare single aberration that has not been previously described in this species. Despite the differences between aberrant signals and normal ones, they could be assigned to one of the species based on the above-mentioned speciesspecific features.
Thus, the acoustic characteristics described here showed that the species specificity of the alarm call of each specimen from the Ulagchin-Gol valley was clearly diagnosed, and the intrapopulation variability of the acoustic signal through aberrations of its frequency-temporal properties was detected. Differentiated spatial distribution of individuals with different species alarm calls was detected along the valley.

Molecular Genetic Data
In order to be used as references, nucleotide sequences of i5HoxB, i13BCR, cytb, and i8SmcY were sequenced from individuals of M. sibirica and M. baibacina from reference populations outside the zone of sympatry (Supplementary Table 1 Figures 3-6). No intraspecific variations in autosomal markers were detected. Species-specific substitutions in i5HoxB, i13BCR, and cytb formed recognition sites for some restrictases. The PmlI restrictase cut sequences of i13BCR of M. sibirica and cytb of M. baibacina at positions 395 and 292, respectively, into two fragments of unequal length. The MspI\HpaII restrictase cut the i5HoxB sequence of M. baibacina at position 305. This sequence was also cut by this enzyme at positions 392, 456, and 579, forming a set of fragments identical for both species (Supplementary Figure 7). The obtained data made it possible to use restriction analysis for molecular testing of the total sample of marmots (58 individuals) from a mixed population for genetic typing by i5HoxB, i13BCR, and cytb loci. Individuals from reference populations of both species are included additionally as a control of test results. Genotyping of 31 males from the mixed population was carried out using i8SmcY based on two species-specific replacements (Supplementary Figure 5) by sequencing.

). Nucleotide variabilities of all molecular markers were presented in Supplementary Materials (Supplementary
A species test of the total sample of marmots from the mixed population was carried out on the basis of the analysis of allele composition. The alleles of both species of each genetic marker were found in the studied sample. An individual was associated with one species of marmots if it had alleles of only one species for all markers. In other cases, the individual was found to be a hybrid. The results are presented in Table 1. According to the results of genetic diagnostics, more than half of marmots (62%) of the examined sample may be classified as one of the species (51.72% M. sibirica, 10.34% M. baibacina). The others (38%) were classified as hybrid marmots and had alleles of both species in different combinations ( Table 1). In the total sample of marmots, the frequency of sibirica nDNA alleles was three times higher on average, and the frequency of sibirica mtDNA alleles almost twice higher than that of baibacina ( Table 2).
The prevalence of nDNA sibirica alleles was also observed in the sample of hybrids. At the same time, maternal lines (cytb) of M. baibacina prevailed in hybrids and paternal lines (i8SmcY) of both species and were represented almost equally ( Table 2). Most of the hybrid males (8 of 11) had different species markers of mtDNA and Y-chromosomes. The total ratio of species-specific haplotypes of sex markers in hybrid males was close to equal (41% sibirica and 59% baibacina), but baibacina mitotypes prevailed in their mtDNA (73%). An approximately equal ratio of speciesspecific alleles was also found in their nuclear genome (59% sibirica and 41% baibacina).
Since the bioacoustic analysis had detected a place with a predominance of boulder screes where individuals with alarm calls of different species live together, a sample of marmots from this area was analyzed separately. The analyses showed a high prevalence of hybrids (62%) in this sample. The number of M. sibirica individuals found here was only twice as high as M. baibacina (25 and 13%, respectively). In this sample, the frequency of diploid marker alleles of sibirica was higher than that of baibacina, and their ratio was similar to that observed in the sample of hybrids ( Table 2). At the same time, the ratio of cytb mitotypes was near to 1:1.
The observed heterozygosity (Ho) for each of the diploid markers of nDNA was lower than expected one (He), while the Wright fixation index (F is ) is >0 in the total sample ( Table 3). The same trend was observed in the sample from the boulder screes place. In this sample, the heterozygosity values in the latter were appreciably higher, and the F is value was lower by half.
The spatial distribution of genotyped individuals is shown in Figure 4. "Pure" tarbagans inhabit the entire studied area but are more attracted to the banks of creeks and slopes. "Pure" gray marmots, which are much less numerous, are concentrated in or near large rocky places. One M. baibacina was also identified in the lower part of the valley where this species predominates according to bioacoustic analysis ( Figure 3D). Hybrids, as well as tarbagans, were present throughout the area but concentrated more in boulder screes (Figure 4 and Supplementary Tables 3, 4).  Of the six family groups in which all or part of the individuals were genetically tested, only two consisted of "pure" M. sibirica ( Table 1); only hybrids were found in two family groups; one group contained a hybrid and M. baibacina and another contained both M. baibacina and M. sibirica.
In total, 57 complete nucleotide sequences of the mtDNA control region (998-1,001 bp) were sequenced to estimate the level of general intra-population genetic variability. Fifteen haplotypes of CR (Table 1 and Supplementary Table 1) were found in the studied sample, including individuals from "pure" populations and from the contact zone. All haplotypes formed two clusters on the haplotype evolutionary network (Figure 5). One cluster combined two haplotypes of "pure" M. sibirica with four haplotypes from the mixed population, while two haplotypes of "pure" M. baibacina formed another cluster together with seven haplotypes from the mixed population. Interspecific differences consisted of 72 nucleotide substitutions. The genetic variability indices for this locus in the sample of marmots from the mixed population are presented in Table 4.
Of the 24 marmots for which specific identification was visually determined in the field, only in 15 individuals did the visual identification match by the results of genetic typing ( Table 1). It should be noted that in 2 out of 9 cases of incorrect field identification, alleles of the species to which it was visually assigned prevailed in specimen genotypes. These results showed that the species identification of these marmots in the hybrid population was not reliable, based on external features.
All five marmots for which sound records were made in cages were not genetic hybrids (three M. sibirica and two M. baibacina).
In summary, the genetic analysis revealed a high proportion of hybrids in the population and their high concentration in boulder screes in studied part of the colony. All new DNA haplotype sequences were deposited in GenBank under accession numbers MT412445-MT412473.

Spatial Analysis and Habitat Assessment Data
Half of the total sample of genetically tested marmots inhabited boulder screes, 38% on the lower parts of mountain slopes and at the bottom of the valley, 12% on the upper parts of slopes and on watersheds (Supplementary Table 3 Table 4) indicate no differences in the frequency of occurrence of pure individuals of both species in all habitats, and a non-random distribution of hybrids across habitats.
A survey of 103 family sites of marmots on the area of 5.71 km 2 was carried out to estimate a population density in different habitats in the upper part of the Ulagchin-Gol valley at the plot where genetic and acoustic surveys were conducted ( Figure 3A). The population density was 48.6 families/km 2 in lower parts of mountain slopes and at bottom of the valleys, 2.4 families/km2 in upper parts of the slopes and watersheds, and 6.4 families/km2 in boulder screes.
Our data did not reveal the influence of altitude on the distribution of marmots with different genotypes. The average habitat altitude of M. sibirica was 2,862 ± 16.8 m a.s.l. in the range of 2,710 m a.s.l. minimum to 3,084 m a.s.l. maximum; M. baibacina was 2,849 ± 47.5 m a.s.l., 2,680-2,980 m a.s.l. (min-max); hybrids was 2,942 ± 21.3 m a.s.l., 2,718-3,062 m a.s.l. (min-max).

Population Genetic Features of M. baibacina and M. sibirica Hybridization
Our data on molecular genetic typing of marmots from the area cohabited by M. baibacina and M. sibirica indicated an active hybridization of these two species. The high frequency of hybrid individuals in the population (38% of the total sample) and the presence in their genomes of almost all possible combinations of parent species alleles of studied markers ( Table 1) indicated viability and fertility of F1 hybrids and viability of their offspring.
The predominance of sibirica alleles in the nuclear genome in both the hybrids and the general sample most likely reflected the prevalence of M. sibirica individuals in the studied area that could lead to more frequent reciprocal crossbreeding with this form. The lack of M. baibacina in cohabitation sites may be the cause of the relatively low observed proportion of F1 hybrids. At the same time, the high proportion of hybrids in the population indicated intensive crossbreeding with both parent species and other hybrids in different combinations. Taking into account the presence of genotypes with almost all possible allele combinations FIGURE 4 | Spatial distribution of genotyped individuals in the area of cohabitation of two species of marmots (it is marked by a red dotted frame in Figure 3). Symbols: orange circle-M. sibirica; blue square-M. baibacina; green rhombus-hybrid; white circle-single-family individuals; red oval indicates the area of boulder screes. Family identification numbers are in square brackets (see Table 1).  in the population, we supposed that crossbreeding of individuals with different genotypes occurred without significant limitations, and F1 or backcross sterility may be absent or very poorly expressed. However, the amount of our data was not sufficient to confirm the fertility of backcrosses. The presence of both species' alleles in hybrids in both the mitochondrial genome (maternal lines) and Y-chromosome (paternal lines) ( Table 2) indicated the involvement of males and females of both species in the hybridization. In this case, a partial or total violation of Haldane's rule, i.e., sterility and non-viability of the heterogametic sex in hybrids, may be assumed, however, this requires more in-depth study. Apparently, the independent evolution of these species during about 2.5-3 million years in allopatric conditions (Steppan et al., 2011) did not lead to the development of effective postzygotic isolation. However, there was a heterozygote deficiency in the mixed Ulagchin-Gol valley population (Ho < He, F is > 0 both for each diploid autosomal marker and on average in all samples (Table 3), indicating a lack of panmixia and population subdivision-the "Wahlund effect" (Wahlund, 1928). Although our data are insufficient to demonstrate the absence of fertility limitations in backcrossing with hybrids, they indicate a high genetic compatibility of these species. High potential for panmyxia in the mixed population is also supported by a smaller F is value and low Ho and He differences in the sample from boulder screes, where the increase in heterozygosity can be prevented by a constant flow of "pure" individuals into this area. Taking this into account, it may be assumed that, besides the possible reduced viability of hybrids, one of the source of population subdivision may be premating mechanisms such as behavioral or ecological features that limit free crossbreeding in a mixed population.

Phenotype and Genotype Intrapopulation Variability
The frequency-time structure of the alarm call is a good diagnostic characteristic for the studied pair of marmot species. The alarm call similar to the aberrant signals of M. baibacina that we found (Figure 2F) was considered a hybrid attribute in earlier studies of this population (Formozov and Nikol'skii, 1986). No genetic analysis has been carried out in this research. We cannot say that the found aberrations were indications of the hybrid origin of the individual, although this cannot excluded. Recently, aberrations of the sound signal that differ from those described by us were found in the Tien Shan population of M. baibacina in the contact zone with M. caudata (Nikol'skii and Vanisova, 2020). However, the authors believe that the formation of these aberrations in the result of hybridization is unlikely. Earlier, aberrations of the sound signal with frequency modulation characteristics similar to those we found in the Ulagchin-Gol valley were found in the closely related species M. bobak (Nikol'skii, 2008). They were variants of intraspecific variability, the frequency of which increases in marginal populations under isolation. The population of the Ulagchin-Gol valley is located near the boundaries of distribution of both contacting species, and the found aberrations of the sound signal may be elements of intraspecific variability that arise under the influence of the same factors as in M. bobak. The relationship between the sound signal aberrations and interspecies hybridization in the peripheral populations of M. baibacina requires additional genetic studies.
The high variability of baibacina phonotype structure corresponds to the variability of baibacina control region of mtDNA haplotypes in the investigated population (Table 4 and Figure 5). According to a central-marginal hypothesis (Eckert et al., 2008), marginal populations are more characterized by a low level of genetic diversity formed by the founder effect and/or genetic drift. The variability of sibirica haplotypes, but not baibacina, corresponds to this hypothesis. Moreover, a high level of genetic variability in marginal populations is often associated with the areal oscillation (Vandermeer, 2006), such as in the greater long-tailed hamster Tscherskia triton under the influence of Pleistocene glacial fragmentation of their area of distribution (Sheremetyeva et al., 2017). We may hypothesize that the high level of phonotype and haplotype variability of baibacina either have been maintained during the reduction of the M. baibacina range accompanied by the loss of monogenic marginal populations and the transition of central populations to the species distribution boundary or formed by gene flows due to cyclic fragmentation-defragmentation of marginal parts of the range during Pleistocene climate changes. Although populations in hybrid zones are often characterized by a high level of variability, we do not have enough data to conclude that high levels of variability in M. baibacina acoustic signal and CR is associated with hybridization.
Individuals of hybrid origin were detected significantly more by molecular genetic typing than by comparing spectrograms. However, the predominance of individuals with the alarm call of M. baibacina (60%) in the sample from boulder screes (Figure 3) with the high proportion of hybrids did not correspond to the observed low frequency of baibacina alleles in this site ( Table 2). Probably, intermediate phenotypic traits are not expressed or weakly expressed in M. baibacina-M. sibirica hybrids, in contrast to that found in S. major-S. fulvus and S. major-S. suslicus hybrids (Titov et al., 2005;Kuzmin et al., 2011). The lower incidence of intermediate phenotypic traits compared to the level of hybridization detected with genetic markers is also known in other mammals, for example, in wolf-dog hybrids (Lorenzini et al., 2014). An understanding of the phenotypic expression of mixed genotypes in marmots requires further investigation.

Spatial, Habitat, and Behavioral Background of Hybridization
The spatial analysis of locations of individuals with different species traits (Figure 3) indicates that two marmot species are distributed unequally in the Ulagchin-Gol valley space. M. sibirica prevails in the upper part and M. baibacina prevails in the lower part of the valley. The optimal habitat for both marmot species in the study area seemed to be the bottom of the valley and the lower part of the slopes. The density of family sites found in this habitat type was four times higher than the average population density (11 families/km 2 ) in other places of Mongolian marmot areas with similar environmental conditions . According to this and results of Fisher's exact test (Supplementary Table 4), marmots in the Ulagchin-Gol valley least preferred watersheds, and boulder screes, although more favorable, are also suboptimal. Predominant occupation of valley bottoms by M. sibirica may indicate that it displaced M. baibacina from these habitats. At the same time, it is known that stony habitats are unfavorable for M. sibirica (Badmaev, 2020), and it is reluctant to inhabit such habitats. On the other hand, M. baibacina is a more "mountainous" species (i.e., more adapted to living in composite topography) than M. sibirica (Bibikow, 1996), and therefore it is more tolerant of living in boulder screes. The observed ecological segregation of two species of marmots is consistent with the found "Wahlund effect" and may be one of the factors limiting panmyxia and, consequently, hybridization in a mixed population. Taking into account the fact that both species of marmots inhabit all three types of habitats outside the area of range overlap (Rogovin, 1992;Brandler et al., 2010b), the detected habitat separation of these species may indicate their biotopic reorientation (change in habitat preferences) as a result of competitive relationships in a hybrid population leading to a decrease in the frequency of interspecific contacts similar to that described for mixed settlements of ground squirrels S. major and S. suslicus (Titov et al., 2012).
The spatial distribution of genetically analyzed individuals (Figure 4) is generally consistent with sound signal analysis. The high incidence of hybrid individuals and their concentration in boulder screes (Supplementary Tables 3, 4) apparently indicates that the most numerous contacts between heterospecific individuals and a mixed pair formation occur precisely in these places as well as subsequent crossbreeding of hybrids. In addition, the family group, which consists of a heterospecific couple of adults [144] ( Table 1 and Figure 4), inhabited close to a large rocky area. Nevertheless, hybrids were found not only in boulder screes but, apparently, they were able to penetrate adjacent areas and coexist with both parent species, however, insufficient data are currently available to quantify a dispersion of hybrids in this population. Previous studies using tagged animals have shown that the migration distance for M. baibacina is 1-4.5 km, but only 13.7% of migrants move more than 1 km (Bibikow, 1996). Based on these data, we assume that the dispersion of hybrids declines rapidly as they move away from their birthplace.
One of the possible hypothetical hybridization scenarios for this population is based on the ecological segregation of different species, high occupancy of habitats, and the speciesspecific structure of colonies, basis of which is formed by stable territorially-conservative family groups not allowing nonresidents to enter their area (Bibikow, 1996;Adiya, 2007;Armitage, 2014). Young animals are also generally banished by their parents from the family site after reaching puberty and they become migrants. Considering the high level of intrapopulation aggression and the low number of interfamily regroupings found earlier in this population (Kolesnikov and Svininykh, 2010), we may assume that non-resident migrants are highly exposed to elimination in this area by predators, diseases and injuries inflicted by resident individuals during territorial conflicts as well as in other marmot species in similar conditions (Suntsov, 1981;Armitage, 2003;Mashkin et al., 2010). As a result, animals are more likely to succeed if they settle close to their native site in suboptimal habitats where mixed couples may be formed in conditions of limited partner choice. In this case, conspecific assortative mating is significantly reduced due to specific copulation time. Male and female of these marmot species choose mates and form pair bonds in the fall prior to hibernation, and copulation occurs in the burrow prior to emergence in the spring (Bibikow, 1996). Thus, conditions of limited mate choice may contribute to the formation of hybrids in suboptimal habitats. The high probability of mixed couples and the short distance of successful dispersion of migrants may explain the detection of higher concentrations of hybrid individuals in boulder screes places. Obviously, it is not the only factor determining the spatial features of these species hybridization. Further studies on the migratory behavior of hybrids may provide the necessary data to better understand this process.
Divergence times of the bobak group including M. baibacina and the camtschatica group including M. sibirica were estimated to be 3.5-4 Ma (Steppan et al., 2011), i.e., the event occurred relatively soon after an initial invasion of marmots into Asia (4.6 Ma). Allopatric speciation in these groups developed during the epi-platformal orogenesis of Eurasia (Nikol'skii and Rumiantsev, 2012) under the influence of glacial and interglacial events of Middle and Late Pliocene and Pleistocene (Erbaeva, 2003;Polly et al., 2015). Taking the above dates into account, we can say that the duration of M. baibacina and M. sibirica speciation corresponds to the average time of evolution of hybrid inviability in mammals (2-4 million years) (Fitzpatrick, 2004), which is consistent with relatively high genetic interspecific differences between them (Steppan et al., 2011). At the same time, as our study shows, intrinsic postzygotic reproductive barriers between these species have not yet been formed. It can be assumed that the mountain glaciations of the Altai formed an effective geographical barrier between these species at the end of the Pliocene and the beginning of the Pleistocene, which agrees with the model of reduction of marmots' ranges during cold glacial cycles (Nikol'skii et al., 1999;Polly et al., 2015). In any case, the recent overlapping of M. baibacina and M. sibirica areas might occur only as a result of their expansion after the last spatial isolation. We assumed that ranges of M. baibacina and M. sibirica were isolated during mountain glaciations. The last climatic cooling in the Mongolian Altai that may have a significant impact on mountain-steppe ecosystems and, consequently, on the distribution of marmots corresponds to the Sartan Glaciation and has occurred at about 25 Ka (Devyatkin, 1981). On the other hand, M. baibacina and M. sibirica may have been isolated during an expansion of forests in the mountains of Mongolia the last of which was dated to the Middle Holocene and ended 4,000 years ago (Vipper et al., 1989). In any case, it may be assumed that the overlapping zone of M. baibacina and M. sibirica areas followed the disappearance of isolation barriers and existed for a relatively long time, possibly with some interruptions related to short-term climatic changes.
Comparing the studied pair of marmot species with other hybridizing Marmotini, it can be noted that the contact between M. baibacina and M. bobak in Kazakh upland also has a relatively long duration resulting from the expansion of M. bobak 5-7 Ka (Nikol'skii et al., 1983). In most other known cases of intense hybridization, contact zones were formed relatively recently as a result of an expansion of one or several species under influence of anthropogenic factors that transform the natural environment. Examples of this scenario are the contact zones of S. major with the neighboring species S. suslicus, S. fulvus, and S. pygmaeus (Ermakov et al., 2002) and I. parvidens-I. tridecemlineatus (Stangl et al., 2012).
The level of hybridization we documented between M. baibacina and M. sibirica is lower than what has been reported in Palaearctic Spermophilus, both in the proportion of hybrids in populations and in spatial coverage (Ermakov et al., 2002(Ermakov et al., , 2015Titov et al., 2005;Spiridonova et al., 2006;Ivanova et al., 2017). But in general, hybridization occurs more frequently in secondary contact zones within Palaearctic Marmota and Spermophilus than Nearctic Marmotini (Hird et al., 2010;Stangl et al., 2012;Frare et al., 2017;Leitner et al., 2017). This may be due to the greater divergence of Nearctic species and, consequently, forming more effective interspecific barriers compared with Palearctic species (Herron et al., 2004;Brandler et al., 2010a). However, neither our study nor those addressing ground squirrel hybrid zones have identified adaptive effects of gene introgression as found in other mammalian hybrid zones [e.g., formation of a hybrid swarm in Odocoileus (Haines et al., 2019); introgressive replacement of the mitochondrial genome in Lepus (Marques et al., 2017)].

CONCLUSION
Our study is the first description and characterization of the hybrid zone in marmots using genetic, bioacoustic, and ecological approaches. We revealed active hybridization in the secondary contact zone of two marmot species M. baibacina and M. sibirica. Genetic data detected no limitations for heterospecific crossbreeding and indicated fertility of F1 and their descendants, but also indicated subdivision of the mixed population, the source of which requires more extensive genetic research. The differential spatial distribution found in the valley and the change in habitat preferences of the two marmot species may limit the free interbreeding of heterospecifics. We previously suppose that behavioral features of marmots associated with long-term territory maintenance and family structure of the colony combined with ecological specificity contribute to the appearance of hybrids. At the same time, aggressive elements of intrapopulation behavior in combination with landscape heterogeneity of surrounding area may limit the dispersion of hybrid individuals reduce the impact of gene introgression on species genomes. Further studies of hybrid migration and gene pools of surrounding marmot populations may provide required data for understanding the gene flow between these hybridizing species. The adaptive significance of hybridization remains unclear in the study population, and further studies might shed light on it.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: GenBank NCBI, accession numbers MT412445-MT412473.

ETHICS STATEMENT
The animal study was reviewed and approved by the Ethics Committee for Animal Research of the Koltzov Institute of Developmental Biology RAS.

AUTHOR CONTRIBUTIONS
OB designed the study and prepared this manuscript. OB, SK, VK, and BB collected tissue samples. SK extracted DNA and performed PCRs and restriction. OB and SK performed molecular data analysis, mapped spatial data, recorded sound signals, and contributed equally to the study. AN and OB performed bioacoustic analysis. VK tagged marmots and performed population density assessment. YA provided guidance and field works in Mongolia. All authors contributed to the editing and finalizing of this manuscript.