Species‐specific habitat preferences do not shape the structure of a crested newt hybrid zone (Triturus cristatus x T. carnifex)

Abstract Reproductive isolation barriers maintain the integrity of species by preventing interspecific gene flow. They involve temporal, habitat or behavioral isolation acting before fertilization, and postzygotic isolation manifested as hybrid mortality or sterility. One of the approaches of how to study reproductive isolation barriers is through the analysis of hybrid zones. In this paper, we describe the structure of a hybrid zone between two crested newt species (Triturus cristatus and T. carnifex) in the southern part of the Czech Republic using morphological, microsatellite, and mitochondrial (mtDNA) markers. Specifically, we tested the hypothesis that the structure of the hybrid zone is maintained by species‐specific habitat preferences. Comparing the genetic structure of populations with geographical and ecological parameters, we found that the hybrid zone was structured primarily geographically, with T. cristatus‐like populations occurring in the northeast and T. carnifex‐like populations in the southwest. Despite T. cristatus tending to occur in deeper ponds and T. carnifex on localities with more shading, the effect of both ecological parameters on the structure of the zone was minimal. Next, we corroborated that T. carnifex individuals and some hybrids possess mtDNA of T. dobrogicus, whose nuclear background was not detected in the studied hybrid zone. Hybridization between T. carnifex and T. dobrogicus (resulting in unidirectional mtDNA introgression) had to predate subsequent formation of the hybrid zone between T. cristatus and T. carnifex. Populations of crested newts in the southern part of the Czech Republic thus represent a genetic mosaic of nuclear and mitochondrial genomes of three species.

manifested by increased mortality or reduced fertility of first and subsequent generation hybrids. The study of prezygotic and postzygotic barriers is therefore crucial to understanding the origin and maintenance of species in nature.
One of the approaches for understanding the evolution of reproductive isolation barriers, interspecific gene flow (introgression) and mechanisms of speciation is through the analysis of hybrid zones.
These transition zones are relatively narrow regions of admixed genotypes originating when two diverged populations (e.g., species, subspecies, or chromosomal races) come into contact, mate, and produce hybrid progeny (Hewitt, 2001). They are characterized by abrupt changes in species-specific traits (alleles in the case of genetic markers) along a geographic transect or cline (Barton, 1983(Barton, , 2001. The width of the cline increases with the dispersal rate of the parental genotypes and decreases with selection against intrinsically incompatible hybrid genotypes (endogenous selection). The dispersal of parental genotypes into the hybrid zone center is the main source of nonrandom associations among loci (linkage disequilibrium). Linkage disequilibrium generates interactions among selected loci that increase the effective selection of any one locus. When hybrid zones occur in ecotones, where both parental species are adapted to different habitats, the steepness of the cline and its width might be, moreover, influenced by selection against genotypes in the alternative habitat (exogenous selection; Arias et al., 2008;Barton & Gale, 1993;Shurtliff, Murphy, & Matocq, 2013;Yanchukov et al., 2006). Habitat isolation between the parental species may thus limit interspecific gene flow and determine the nature of species boundaries. In this study, we focused on hybridization between newt species belonging to the Triturus cristatus superspecies group in order to determine whether habitat preferences shape the structure of the hybrid zone and thus can be involved in species isolation. We tested whether there is an association between the genotypic composition of populations and the type of aquatic habitats which newts use for reproduction.  (Mikulíček, Horák, Zavadil, Kautman, & Piálek, 2012;Piálek, Zavadil, & Valíčková, 2000;Reiter & Hanák, 2000;Zavadil, Piálek, & Klepsch, 1994). The most complex population structure was observed in the Znojmo region, where T. cristatus comes into contact with T. carnifex. The majority of T. carnifex individuals and their hybrids with T. cristatus possess an introgressed mitochondrial genome of T. dobrogicus (Mikulíček et al., 2012), the nearest populations of which are, however, located 40 km to the east of the Znojmo transect. The hybrid zone in southern Moravia is several kilometers wide (the shortest straight geographic distance between T. cristatus and T. carnifex populations is ca 15 km), but the mechanisms preventing hybridization between the species and maintaining the width of the zone are unknown.
In the present study, we focus on species-specific habitat preferences between T. cristatus and T. carnifex in a hybrid zone situated in southern Moravia. We assume that if the hybrid zone were shaped by ecological preferences of the crested newt species, genetic composition of their populations should correlate with specific habitat characteristics. Specifically, we aim (a) to describe the structure of the hybrid zone based on morphological, microsatellite, and mitochondrial markers using more extensive sampling than in a previous study (Mikulíček et al., 2012) and (b) to find out whether there is an association between genotypic composition of populations and aquatic (reproductive) habitats of newts.

| Sampling
In this study, we focused on the region of southern Moravia (Czech Republic), where a contact zone between three crested newt species has been documented (Mikulíček et al., 2012). Individuals (n = 300) were caught on 38 sampling sites during the breeding season between the years 2010-2015 (Table 1). Funnel collapsible nylon traps with bait (chicken liver) were used for catching, following the methodology of Bock, Hennig, and Steinfartz, (2009). A toenail clip was removed and stored in 96% ethanol.

| DNA markers and laboratory techniques
Population genetic structure and the degree of hybridization between the crested newts were inferred by two types of genetic markers: biparentally inherited nuclear microsatellites and maternally inherited mtDNA. Sixty-six individuals from the Moravian hybrid zone were analyzed for mtDNA. We amplified the >1,400 bp-long portion of mtDNA comprising the complete ND2 gene, five subsequent transfer RNA (tRNAs) genes and the light-strand replication origin using primers (L3780, H5018) and protocol following Krupa et al. (2002).
The final analyzed stretch contained a 620 bp-long fragment of ND2. The novel sequences were deposited in GenBank under accession numbers: MN394474-MN394541. The ND2 fragment was aligned using the Clustal W algorithm (Thompson, Ling, & Grustein, 1994) as implemented in BioEdit (Hall, 1999). Alignments were checked by eye and low-quality ends were trimmed. Ambiguously aligned regions/ gaps were ignored for the subsequent analysis. No stop codons were detected when the sequences were translated using the vertebrate mitochondrial genetic code in the program DnaSP 5.10. We used a network approach (Posada & Crandall, 2001) French et al., 2014) and the implemented median-joining algorithm.
Seven microsatellite loci (Krupa et al., 2002)  all individuals were assigned to the inferred clusters (from K = 1 to K = 10) without any a priori population information. In a second procedure, a priori population information for individuals from reference populations was used. Individuals from populations in or close to a presumable contact zone were assigned to the clusters without using any a priori information. To evaluate convergence and to estimate optimal genetic clustering, five replicates were run for each K value using an admixture and uncorrelated allele model. All Structure analyses were based on runs of 10 6 iterations, following a burn-in period of 10 4 iterations. The number of populations that best fitted the data set was defined by the ΔK method (Evanno, Regnaut, & Goudet, 2005), as implemented in Structure Harvester (Earl & vonHoldt, 2012). For each individual, we calculated the q values which define the proportion of an individual's genome that originated from cluster K. The q values were calculated also for populations as a mean of individuals' q values.
In Geneland, we analyzed only samples from the presumable contact zone (southern part of the Czech Republic). First, we ran analyses with K free to vary, to infer the optimal value of this parameter. We ran the analysis 10 times to verify the consistency of the results, with the following parameters: 500,000 MCMC iterations, maximum rate of Poisson process fixed to 100, zero uncertainty of coordinates, minimum and maximum K 1 to 10, maximum number of nuclei in the Poisson-Voronoi tessellation fixed to 300, and null allele and uncorrelated allele models. All 10 replicates revealed the maximum posteriori estimate of K = 2. Then, we ran the MCMC 10 times with K fixed to 2 and the same parameter settings. The run with the highest log probability was chosen for postprocess analyses. The posterior probability of population membership for each pixel of the spatial domain and for each individual was then computed with 500 pixels along the X and Y axes. The modal population of each individual, maps of population membership, and maps of probability of population membership were finally computed.
To estimate the level of genetic diversity, we applied the program GenAlEx 6.5 (Peakall & Smouse, 2012) for calculation of the number of alleles, the coefficient of inbreeding (F IS ), and observed (H O ) and expected (H E ) heterozygosity. These parameters were calculated for each population under study and then for three groups in the contact zone, that is "pure" populations of T. cristatus, "pure" populations of T. carnifex and hybrid populations. Assignment of populations to these groups was based on an admixture proportion of each individual (parameter q according to Pritchard et al., 2000, i.e., the proportion of an individual's genome that originates from the T. cristatus and T. carnifex cluster) averaged for a population. Populations were considered as hybrid when they possessed more than 20% of introgressed microsatellite alleles (i.e., average population q value was ≤0.8).

| Morphological characteristics and comparison with microsatellite markers
Morphological differences between T. cristatus and T. carnifex are conspicuous. They differ in the body habitus, and the relative length of the trunk and legs (Arntzen & Wallis, 1999). These differences might be expressed as the "Wolterstorff" index, or WI (Wolterstorff, 1923) which is defined as the ratio between forelimb length (Pa) and interlimb distance (LiE; WI = 100 × Pa/LiE). The values of WI increase from T. cristatus to T. carnifex.
Morphological characteristics were obtained from anesthetized newts (0.8% solution of 2-phenoxyethanol) using a dial-caliper (with an accuracy of 0.5 mm). We measured the following characteristics: L-body length, Lcd-tail length, Ltot-total body length, Lc1-jaw length, Lc2-head length, Ltc-head width, Pa-front limb (on both sides of the body), Pp-hind limb (on both sides of the body), LiEinterlimb distance (on both sides of the body). All 12 characteristics were measured by the same caliper and same person.
We used RDA (Redundancy analysis) to find out which morphologi- Again, we used the q values from Structure to assign individuals to a particular species. An individual was assigned to T. cristatus or T. carnifex when the proportion of its genome (the q value) originating from T. cristatus or T. carnifex cluster was equal or higher than 0.8.

| Association between genotypic composition of populations and habitat characteristics
In order to test preferences of particular species to specific habitats, we compared genotypic composition of populations with environmental variables. For each population, we calculated average q values from Structure, which were entered (after the log transformation) to the analysis as response variables.
Seven habitat characteristics of ponds where crested newts reproduced, plus geographic coordinates (latitude and longitude), were recorded in April 2014. These habitat characteristics are relevant for the distribution of crested newts as was found in previous studies (Maletzky, Kyek, & Goldschmid, 2007). Specifically, we measured pond area (m 2 ), maximum depth in three classes (<30 cm, 30-100 cm, >100 cm), fish presence or absence, origin of pond (natural or artificial), and presence or absence of human use. We also estimated density of submerged vegetation (in 25%-classes) and proportion of shade (in 25%-classes). These characteristics (including latitude and longitude) were tested to collinearity and then entered into the analysis as explanatory variables. is distributed (Figure 1). One haplotype of T. dobrogicus was unique for the study area but diverged by just two mutation steps from the vastly distributed haplotype.

| Microsatellite markers
Out of the seven microsatellite loci, two (Tcri29 and Tcri35) revealed many missing genotypes and difficult-to-interpret patterns in a fragmentary analysis. Therefore, these loci were excluded from population genetic analyses. Structure Harvester using the ΔK method estimated the most likely number of clusters for K = 2 (ΔK = 960.5). Cluster 1 and cluster 2 corresponded to T. cristatus and T. carnifex, respectively. The second highest optimal cluster number, for K = 6, had much lower probability (ΔK = 32.7). Newts from reference populations were assigned to the correct Structure clusters (corresponding to the parental species) with an average probability q = 0.978 for T. carnifex (the only reference locality Matena) and q = 0.837-0.986 for T. cristatus populations (Table 2).
Specimens originating from the contact zone were either assigned to the parental species or showed mixed ancestry between T. cristatus × T. carnifex (Table 2, Figure 2). The hybrid zone forms a sharp and geographically structured cline, with T. cristatus-like populations located in the northeast and T. carnifex-like populations in the southwest.
Geneland corroborated results from Structure and revealed two clusters corresponding to the hybridizing parental species (Table 2).
Cluster 1 and cluster 2 corresponded to T. carnifex-like and T. cristatuslike populations, respectively (Figure 3). The transition of genotypes from one cluster into the other is abrupt, showing that the contact zone between both species is a narrow region where species hybridize.
Parameters of genetic diversity for each population are summarized in Table 2. "Pure" populations of the parental species and hybrid populations from the contact zone revealed a comparable level of genetic diversity (Table 1).

| Morphological characters and their comparison with microsatellite markers
Partial RDA analysis explained 28.60% of variation in samples, the first axis explained 28.39% of variability, the second axis 0.21%, and the third axis explained 65.38% of variation. Test on the first axis

| Association between genotypic composition of populations and habitat characteristics
Genotypic composition of populations based on microsatellites (q values per population in Table 2) was compared with habitat characteristics. Within all geographical and environmental variables, only latitude, longitude, depth, and shading were statistically significant and thus were included into the final model. The model explained 82.5% of variation in the data. The most important factor was latitude, explaining 76.2% of variability; followed by longitude, explaining 5.1%. Effects of depth and shading were minimal but their simple term effects were statistically significant. Triturus cristatus tended to occur in deeper ponds, while T. carnifex occurred on localities with more shading (Figure 4). role of species-specific habitat preferences in the architecture of hybrid zones has not been studied in detail. In this paper, we extended our previous study (Mikulíček et al., 2012)
In this study, we tested whether morphological characteristics discriminate between T. cristatus and T. carnifex. Even though eight F I G U R E 2 Geographical distribution and proportion of admixture (parameter q according to Structure) between T. cristatus (red) and T. carnifex (orange) in south Moravia, Czech Republic estimated on the basis of microsatellite data out of ten morphological traits significantly discriminated between the "pure" T. cristatus and T. carnifex individuals, the Wolterstorff index (WI) as a traditional marker used for crested newt delimitation (Arntzen & Wallis, 1999;Wolterstorff, 1923) failed in the studied hybrid zone. More than a third of T. carnifex and a half of T. cristatus individuals were misclassified according to the WI. This finding corresponds to Arntzen and Wallis (1999), who showed 31% misclassification of newts using WI. On the contrary, Brede et al. (2000), who morphologically evaluated hybrids between introduced T. carnifex and native T. cristatus in UK, were able to recognize parental species and F1 hybrids based on morphological traits. In the Czech hybrid populations, however, the existence of different hybrid categories (mainly backcross hybrids with one or another parental species and F2 hybrids) and extensive introgression (Mikulíček et al., 2012) prevent the use of WI in species identification.

| Species-specific habitat preferences
One of the characteristics that defines species in nature is their ecological differentiation involving habitat preferences. In phylogenetically closely related species, ecological differentiation, however, could be very subtle. We evaluated the habitat preferences of F I G U R E 3 Maps of posterior probability of population membership using Geneland. The northerly distributed cluster 1 corresponds to T. cristatus, the southerly distributed cluster 2 corresponds to T. carnifex. A hybrid zone is shown as a sharp cline.
The axes indicate the longitude (X coordinates) and latitude (Y coordinates) Note: Two characteristics (interlimb distances on both sides of the body, LiE1, and LiE2) were removed from the RDA analysis because of their collinearity. Abbreviations: L, body length; Lc1, jaw length; Lc2, head length; Lcd, tail length; Ltc, head width; Ltot, total body length; Pa, front limb (on both sides of the body); Pp, hind limb (on both sides of the body).

F I G U R E 4
Ordination diagram of RDA model with genotypic composition of populations as response variables and habitat characteristics as explanatory variables. Latitude was the most important factor, explaining 76.2% of variability T. cristatus and T. carnifex comparing the habitat features of the studied localities with the genotypic composition of newt populations.
The authors however did not discriminate between the hybridizing species. They found that shading and density of submerged vegetation had significant effects on pond occupancy (Maletzky et al., 2007). Blab and Blab (1981), besides shading and density of submerged vegetation, identified also pond area as the relevant characteristic for the abundance of crested newts. Other authors showed other key habitat features, including water chemistry and the structure of terrestrial habitat (Cooke, Cooke, & Sparks, 1994;Harper et al., 2018;Sztatecsny, Jehle, Schmidt, & Arntzen, 2004). Arntzen and Wallis (1999)

| Structure of the hybrid zone
The absence of marked species-specific habitat preferences indicates that the structure of the hybrid zone between T. cristatus and T. carnifex in the Znojmo region is not shaped ecologically. The zone is structured primarily geographically with cristatus-like genotypes in the northeast and carnifex-like genotypes in the southwest. We assume that this pattern might have been established when the two species came into contact during their recolonization of Central Europe after the last glaciation. Pleistocene glaciations led to retraction of ranges of many European species which survived in refugia (i.e., regions with suitable climatic and ecological conditions) located southerly. At the end of the Pleistocene, when ice-sheets retreated and climatic conditions became more hospitable, species dispersed out of refugia and colonized the newly available areas further north. In the case of the studied species, T. cristatus had a glacial refugium in the Carpathians (Wielstra, Babik, & Arntzen, 2015) and thus had to colonize the Znojmo region postglacially from the west. Triturus carnifex had to reach Central Europe from the south, while its glacial refugia were probably located in the Adriatic region (Canestrelli, Sacco, & Nascetti, 2011). When the two newt species came into contact, they established a secondary hybrid zone which is currently geographically structured from the northeast to the southwest (Hewitt, 2001(Hewitt, , 2011 to produce other generations of hybrids, but these hybrids reveal higher mortality, disturbed meiosis, and production of dysfunctional gametes Callan & Spurway, 1951;MacGregor et al., 1990 and references herein). The structure of the hybrid zone between T. cristatus and T. carnifex thus might be maintained by a balance between the dispersal of parental genotypes into the center of the zone and selection against hybrids; that is, two mechanisms playing an important role in the tension zone model (Barton & Hewitt, 1985;Key, 1968;Macholán et al., 2007;Szymura & Barton, 1986. The width of the T. cristatus and T. carnifex hybrid zone in southern Moravia was estimated by Mikulíček et al. (2012) at ca 15 km.
The shortest straight geographic distance between T. cristatus-like (locality 18) and T. carnifex-like (locality 11) populations in this study (covering more populations and individuals) was ca 17 km. The fact that hybridization and interspecific gene flow are restricted to a relatively narrow transect depends (besides selection against hybrids) on the limited dispersal abilities of newts (Arntzen & Wallis, 1991;Meiling et al., 2015) and their strong site fidelity (Mori et al., 2017).

CO N FLI C T O F I NTE R E S T S
The authors declare that they have no competing interests.

AUTH O R CO NTR I B UTI O N S
The authors have made the following declarations about their contributions: ZM, MR, and PM conceived and designed the study. ZM, AR, and LJ collected data in field. ZM, DJ, SR, and PM analyzed the data. ZM, DJ, and PM wrote the manuscript. All authors participated in the scientific discussions, revised the manuscript, and read and approved the final manuscript.