A Microsatellite Genotyping-Based Genetic Study of Interspecific Hybridization between the Red and Sika Deer in the Western Czech Republic
Abstract
Although inter-species hybrids between the red and sika deer can be phenotypically determined only exceptionally, there is the eventuality of identification via molecular genetic analysis. We used bi-parentally inherited microsatellite markers and a Bayesian statistical framework to re-examine the proportion of hybrids in the Czech red and sika deer populations. In total, 123 samples were collected, and the nuclear dataset consisted of 2668 allelic values. The number of alleles per locus ranged from 10 (BM1818) to 22 (BM888 and T193), yielding the mean of 16 alleles per locus across the deer. The mean allelic diversity of the red deer markedly exceeded that of the Japanese sika deer. Interspecific hybrids were detected, enabling us to confirm the genetic introgression of the sika deer into the red deer populations and vice versa in western Bohemia. The mean hybrid score equaled 10.6%, with 14.3% of the hybrids being among red deer–like individuals and 6.7% among sika-like ones. At two western Bohemian locations, namely, Doupovské hory and Slavkovský les, the total percentages of hybrid animals equaled 18.8 and 8.9, respectively. No red deer alleles were detected in the sika populations of the subregions of Kladská, Žlutice, and Lány. The NeighborNet network clearly separated the seven red and sika deer sampling populations according to the geography. The knowledge gained from the evaluated data is applicable in hunting management to reduce hybridization with the European deer.
Full text
animals Article A Microsatellite Genotyping-Based Genetic Study of Interspecific Hybridization between the Red and Sika Deer in the Western Czech Republic Lenka ŠtohlováPutnová1, Radek Štohl 2,* , Martin Ernst 3and Kateˇrina Svobodová3 Citation: Štohlová Putnová, L.; Štohl, R.; Ernst, M.; Svobodová, K. A Microsatellite Genotyping-Based Genetic Study of Interspecific Hybridization between the Red and Sika Deer in the Western Czech Republic. Animals 2021,11, 1701. https://doi.org/10.3390/ani11061701 Academic Editors: Koichi Kaji and Dimitrios Bakaloudis Received: 24 March 2021 Accepted: 1 June 2021 Published: 7 June 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Department of Animal Morphology, Physiology and Genetics, Faculty of AgriScience, Mendel University in Brno, Zemˇedˇelská1, 613 00 Brno, Czech Republic; [email protected] 2Department of Control and Instrumentation, Faculty of Electrical Engineering and Communication, Brno University of Technology, Technická12, 616 00 Brno, Czech Republic 3Department of Forest Protection and Wildlife Management, Faculty of Forestry and Wood Technology, Mendel University in Brno, Zemˇedˇelská3, 613 00 Brno, Czech Republic; [email protected] (M.E.); [email protected] (K.S.) *Correspondence: [email protected].cz Simple Summary: The sika deer is a very flexible invasive species, capable of living dynamically in both large forests and mixed environment characterized by a prevalence of agricultural land. The Japanese sika deer was introduced to the Czech Republic at the end of the 19th century. The success of an introduced alien species may consist in their hybridizing with closely related taxa. Where few barriers to gene flow exist, rapid introgression of genetic traits from one species into another frequently occurs. The current Czech sika populations embody the most abundant and expanding group in continental Europe. In western Bohemia, we confirmed the interspecific hybridization with the native red deer. In this context, the red deer gene pool is endangered. The animals proliferate steadily in all directions and will most probably spread all over the Czech Republic if no major, timely changes in game management are adopted. Abstract: Although inter-species hybrids between the red and sika deer can be phenotypically determined only exceptionally, there is the eventuality of identification via molecular genetic analysis. We used bi-parentally inherited microsatellite markers and a Bayesian statistical framework to reexamine the proportion of hybrids in the Czech red and sika deer populations. In total, 123 samples were collected, and the nuclear dataset consisted of 2668 allelic values. The number of alleles per locus ranged from 10 (BM1818) to 22 (BM888 and T193), yielding the mean of 16 alleles per locus across the deer. The mean allelic diversity of the red deer markedly exceeded that of the Japanese sika deer. Interspecific hybrids were detected, enabling us to confirm the genetic introgression of the sika deer into the red deer populations and vice versa in western Bohemia. The mean hybrid score equaled 10.6%, with 14.3% of the hybrids being among red deer–like individuals and 6.7% among sika-like ones. At two western Bohemian locations, namely, Doupovskéhory and Slavkovskýles, the total percentages of hybrid animals equaled 18.8 and 8.9, respectively. No red deer alleles were detected in the sika populations of the subregions of Kladská, Žlutice, and Lány. The NeighborNet network clearly separated the seven red and sika deer sampling populations according to the geography. The knowledge gained from the evaluated data is applicable in hunting management to reduce hybridization with the European deer. Keywords: genetic structure; introgression; microsatellite variability; hybridization; sika-red deer hybrid 1. Introduction Hybridization is a common phenomenon in plants, birds, fish, and many other taxa [1–3] . When, however, hybridization threatens native species, it becomes imporAnimals 2021,11, 1701. https://doi.org/10.3390/ani11061701 https://www.mdpi.com/journal/animals
Animals 2021,11, 1701 2 of 16 tant to understand the evolutionary mechanisms and develop sensible management and conservation measures [4–6]. The human translocation (or illegal introduction) of an exotic non-native species into an ecosystem [ 7 ] often results in anthropogenically induced hybridization between the introduced species and the related native genera [ 8 – 10 ]. In this context, we focused on the red deer (Cervus elaphus Linnaeus 1758), which is native to Europe, and the sika deer (Cervus nippon Temminck 1838), an originally eastern Asian species introduced to diverse parts of the world outside its native range [ 11 , 12 ]. While red deer and sikas differ in body size, a variety of other phenotypic traits, and chromosome number (red 2N = 68; sika 2N = 64, [ 13 ]), the hybrids are fertile, and, therefore, can backcross with their parental species both in captivity and in the wild [ 8 , 12 , 14 – 17 ]. The offspring of such multiple crosses are very difficult to identify according to the morphological traits; essentially, their morphological differentiation from the standard sika/red deer is unfeasible [ 12 , 18 , 19 ]. However, modern molecular genetic techniques help us to detect the proportion of hybrid individuals [ 15 – 17 , 20 – 26 ]. At present, in the Czech Republic, molecular genetic analyses are nevertheless applied only scarcely within this topic, especially due to being costly and logistically problematic. This paper was therefore conceived to examine the presumed crossbreeding of sikas with native red deer by exploiting a genetic panel of microsatellite markers. Microsatellites are popular markers for the study of geographical structure and gene flow because they are co-dominant, multi-allelic, and abundant with a wide genome coverage [ 22 , 27 ]. To date, the hybridization of the deer in free-ranging areas has been explored mainly in the UK, where the morphometric and subsequent genetic analyses proved that the actual existence of the hybrids constitutes not a merely potential threat, but rather a real, unexpectedly widespread problem [ 15 , 16 , 20 , 23 ]. The presence of sika-red deer hybrids has also been confirmed in Poland [ 21 ], Germany [ 28 ], and Lithuania [ 24 ] via modern molecular methods. Importantly, pure species (native and introduced) may no longer exist: the evolutionary trajectory of a species may be radically changed through the introgression of large numbers of alleles, a process that effectively endangers that species’ genetic integrity [17]. The European red deer ranges among the largest representatives of its family: After the European elk (Alces alces), it is the largest wild ungulate in Bohemia and Moravia. Besides the roe deer (Capreolus capreolus), the European red deer embodies one of the original even-toed ungulate species in Bohemia and Moravia. In some areas, the density dropped to zero, with the population wiped out and subsequently recovering through incoming animals of foreign origin [ 29 ]. Considering the Czech Republic, the red deer, originally a forest steppe species, now live mainly in large forests at medium and high altitudes. A distribution map of the red deer can be found in a dedicated paper by Andˇera [ 30 ]. While the number of individuals hunted down in the 1960s and 1970s reached approximately 8000, in 1988, the total count already exceeded 20,000, a figure second only to that of almost 27,000 in 1993 (these variations being due to political changes in the country). Then, there followed a drop back to around 20,000 individuals; in 2019/2020, the indicator equaled 29,017 [29,31]. The Japanese sika deer is considered an exotic invasive alien species, introduced to the Czech Republic only at the end of the 19th century; the first import took place in 1891 [ 32 ]. In the late nineteenth and early twentieth centuries, sikas were frequently delivered to populate deer parks [ 11 , 12 ]. Between the early 1920s and late 1940s, namely, during a period of major political turbulence in the country, some deer parks were destroyed, and the sikas escaped [ 33 , 34 ]. The species then spread farther, especially from the parks in the Bohemian region of Pilsen and the Moravian subregion of Bouzov. This process is reflected in the current distribution of sikas, which have become particularly prominent within western and northwestern Bohemia and, to a slightly lesser extent, northwestern Moravia. The constantly very high population size and density of the species increase annually, and the overall occupied area expands rapidly.
Animals 2021,11, 1701 3 of 16 Evidence from the 1950s indicates that the largest population was hosted by Plzeˇnská pahorkatina (The Pilsen Uplands) [ 35 – 37 ] and that the effort to cull the animals to quickly attain the standardized game stock eventually prompted the sikas to spread beyond the area in a relatively rapid manner. The species then gradually invaded ˇ Ceskýles (The Upper Palatinate Forest), Slavkovskýles (Slavkov Forest), Doupovskéhory (The Doupov Mountains) and, partially, Brdy (The Brdy Hills), Slavkovskýles and Doupovskéhory have included sika populations since the early and middle 1960s. In the latter region and its subareas, the species did not constitute a major concern initially, thus gaining dynamic potential for a future surge in the game stocks and transforming the region into a refuge sought after by surplus animals. Over two decades, the originally attractive game became a virtually unsolvable problem [ 35 ]. Although the necessity to adopt convenient measures without delay to prevent undesired territorial expansion and stock growth was emphasized by Wolf and Vavrunˇek [ 38 ] already almost half a century ago, the situation has clearly not improved: In sharp contrast, sikas have still moved farther to increase in number across western Bohemia, and the trend apparently continues [37,39]. Between 2003 and 2010, the territorial spread of sikas locally distributed across the Czech Republic involved, on average, a further 55 thousand hectares every year [ 40 ]. At present, this species is found over almost a third of the country. In 2016, compared to 2006, the territory had enlarged by 19%, and the culling rate had risen by an incredible 130%, from 6200 to 14,400 animals [ 41 ]. The current Czech sika populations embody the most abundant and expanding group in continental Europe [ 12 , 37 , 42 ]; such conditions are made possible mainly by ample food sources and shelter. In this context, however, it also has to be emphasized that a very significant factor consists in the applied game management approach: in terms of hunting, sikas are widely preferred, prevailing on many hunting grounds to improve the shooting opportunities and to intensify the related benefits. An updated, sika-related map illustrative of the current situation in the Czech Republic was presented by Andˇera [ 43 ]. The regional count of the deer increased exponentially during the last several decades, with the real population numbers probably exceeding those estimated and reported in the yearly hunting statistics [ 39 ]. Currently, the Czech Republic’s sika population, based on hunters’ reports for 2019/2020 includes 17,535 individuals [ 31 ]. Sikas may generate considerable environmental impact due to their capability of causing significant damage in forestry, agriculture, and habitat structures [ 36 , 44 , 45 ]. With regard to their anatomical and behavioral features, sikas seem to compete successfully with local autochthonous deer species [ 38 , 39 , 46 , 47 ]. Concerning the Czech Republic, the interspecific hybridization has been documented via morphological [ 12 , 18 , 46 ] and genetic [ 22 ] analyses. Even though the spontaneous hybridization with the red deer has been known since the second half of the 19th century [ 12 , 48 ], it was long marginalized, belittled, or even purposely ignored [39]. Where the admixture degree is not conspicuously visible, it becomes difficult to differentiate between a standard sika/deer type and a hybrid. In such conditions, microsatellites and SNPs constitute strong DNA markers for identifying an individual or species. Regarding this fact, it should be emphasized that the main objective of our study was to quantify the genetic diversity and structure of locally adapted deer populations in western Bohemia by utilizing modern genetic techniques. For this purpose, we first collected the population data of the nuclear DNA (nDNA) markers to describe the genomic variation of the animals originating from different sites within the investigated region. Then, we searched for possible genetic admixture by using Bayesian clustering analysis. Finally, we also examined the phylogenetic analysis and constructed the NeighborNet dendrograms from the genetic distances between the sampling sites. 2. Materials and Methods 2.1. Individuals in the Study We evaluated 123 animals from the red (n= 63) and sika (n = 60) deer populating different western and central Bohemian localities (Figure 1) within the Karlovy Vary and
Animals 2021,11, 1701 4 of 16 Central Bohemian administrative regions; in the latter case, the samples were collected in the Lány deer park. The sika deer sampling sites were as follows: the Lány deer park (n = 13); Žlutice (n = 8); the Kladskáforestry enterprise (n = 9); Slavkovskýles (n = 20); and Doupovskéhory (n = 10). The red deer sampling took place in Slavkovskýles (n = 25) and Doupovskéhory (n = 38). The samples originated from legally hunted animals. The species were determined by the hunters, based on the morphological traits (e.g., body size, antlers, and coloration). The tissue samples were collected into plastic bags and frozen. The actual collection was carried out in accordance with the laws and ethical guidelines established and accepted in the Czech Republic. Animals2021,11,x4of17 2.MaterialsandMethods 2.1.IndividualsintheStudy Weevaluated123animalsfromthered(n=63)andsika(n=60)deerpopulating differentwesternandcentralBohemianlocalities(Figure1)withintheKarlovyVaryand CentralBohemianadministrativeregions;inthelattercase,thesampleswerecollectedin theLánydeerpark.Thesikadeersamplingsiteswereasfollows:theLánydeerpark(n= 13);Žlutice(n=8);theKladskáforestryenterprise(n=9);Slavkovskýles(n=20);and Doupovskéhory(n=10).ThereddeersamplingtookplaceinSlavkovskýles(n=25)and Doupovskéhory(n=38).Thesamplesoriginatedfromlegallyhuntedanimals.Thespe‐ ciesweredeterminedbythehunters,basedonthemorphologicaltraits(e.g.,bodysize, antlers,andcoloration).Thetissuesampleswerecollectedintoplasticbagsandfrozen. Theactualcollectionwascarriedoutinaccordancewiththelawsandethicalguidelines establishedandacceptedintheCzechRepublic. Figure1.AmapoftheCzechRepublicshowingthesamplinglocationsintheredandsikadeer. 2.2.DNAExtractionandGenotypeDetermination ThetotalgenomicDNAwasisolatedfromthebloodortissuesamplesbyusingthe QIAamp®Blood/TissueKit(Qiagen,Valencia,CA,USA)accordingtotheprotocolhand‐ book.TheextractedDNAwasvisualizedona1%TAEagarosegelstainedwithethidium bromide.Thesampleswerethendenaturedat95°Cfor7minandkeptat−20°Cuntilthe genotypingstarted.Weoptimizedapanelof11autosomalmicrosatellitemarkers.The allelesofthemicrosatellitelociBM888,BM1818,ETH225,RM188,OarFCB5,RT1,RT13, T26,T156,T193,andT501[49–54]wereexaminedviamultiplexpolymerasechainreaction (Veriti®Thermalcycler;LifeTechnologies,Carlsbad,CA,USA)andbyapplyingfragment analysis(ABIPRISM310TMGeneticAnalyzer;LifeTechnologies,Carlsbad,CA,USA).A numericalnomenclaturewasemployedfortheallelesizedesignation(inbp)inaccord‐ ancewithourinternalvalues.Wealsousedaninternalcontroltocalibratetheallelesizes ateveryrunastheISAG(InternationalSocietyforAnimalGenetics)comparisontestsare notavailabletogenotypenonmodelorganismsinthestandardization. 2.3.StatisticalAnalysesoftheGeneticData Thebasiclocus‐specificdiversitymeasuressuchasthenumberofallelesandtheir frequencies,polymorphicinformationcontent(PIC),observedheterozygosity(HO),and Figure 1. A map of the Czech Republic showing the sampling locations in the red and sika deer. 2.2. DNA Extraction and Genotype Determination The total genomic DNA was isolated from the blood or tissue samples by using the QIAamp ® Blood/Tissue Kit (Qiagen, Valencia, CA, USA) according to the protocol handbook. The extracted DNA was visualized on a 1% TAE agarose gel stained with ethidium bromide. The samples were then denatured at 95 ◦ C for 7 min and kept at − 20 ◦ C until the genotyping started. We optimized a panel of 11 autosomal microsatellite markers. The alleles of the microsatellite loci BM888,BM1818,ETH225,RM188,OarFCB5,RT1, RT13,T26,T156,T193, and T501 [ 49 – 54 ] were examined via multiplex polymerase chain reaction (Veriti ® Thermal cycler; Life Technologies, Carlsbad, CA, USA) and by applying fragment analysis (ABI PRISM 310TM Genetic Analyzer; Life Technologies, Carlsbad, CA, USA). A numerical nomenclature was employed for the allele size designation (in bp) in accordance with our internal values. We also used an internal control to calibrate the allele sizes at every run as the ISAG (International Society for Animal Genetics) comparison tests are not available to genotype nonmodel organisms in the standardization. 2.3. Statistical Analyses of the Genetic Data The basic locus-specific diversity measures such as the number of alleles and their frequencies, polymorphic information content (PIC), observed heterozygosity (H O ), and expected heterozygosity (H E , or gene diversity), were calculated separately for each species and locus by using PowerMarker version 3.25 [ 55 ]. As the groups being compared are species rather than true populations, they may not necessarily meet the Hardy–Weinberg expectation. The program FSTAT 2.9.4 [ 56 ] was employed to estimate the Wright’s F-
Animals 2021,11, 1701 5 of 16 statistics and the allelic richness (AR) by using the rarefaction method to correct differences in the population sizes. The testing of the linkage disequilibrium between all pairs of loci in each species was conducted with the same software. The frequencies of the null alleles (F Null ) in each microsatellite locus were estimated via Genepop 4.2.1 [ 57 ]. The genetic distances (D A [ 58 ]; or D R [ 59 ]) between the studied sampling sites were computed with PowerMarker and visualized via a NeighborNet dendrogram in SplitsTree 4.13.1 [60]. We investigated the population structure and allocated the individual hybrid scores by utilizing a Bayesian clustering algorithm implemented in the software package STRUCTURE 2.3.4 [ 61 ]. The most likely number of populations in the dataset (K) was estimated through five independent replicates of K= 1–5. The model was run by utilizing a burn-in period of 5 × 10 4 and a cycle of 15 × 10 4 Markov chain Monte Carlo steps (MCMCs), under the standard admixed ancestry model and the correlated allele frequency model with the default parameter ( λ = 1) to analyze the dataset. To categorize the animals into pure or hybrid sika and red deer, we determined whether the 95% confidence intervals for the Qscores overlapped at 0.01 or 0.99 in the pure sika and red deer, respectively. The research involved more than 30 observations, and the data followed an approximately normal distribution (a bell curve), meaning that we can use the z-distribution in the test statistics. The alpha value of p< 0.05 was employed to ensure statistical significance. In a two-tailed 95% confidence (or credible) interval (CI), the alpha value amounted to 0.025, and the corresponding critical value equaled 1.96. Thus, to calculate the upper and lower bounds of the confidence interval, we can take the mean ± 1.96 × standard deviations from the mean. Using STRUCTURE HARVESTER v0.6.94 [ 62 ], we calculated the mean likelihood, ln P(K); the standard deviation for each value of K; and ∆ K, the second-order rate of change of the likelihood with respect to K[ 63 ]. The results were then software processed to generate input files to be used with the CLUMPP application. CLUMPP 1.1.2 [ 64 ] then estimated the maximum rate of similarity between the Q-matrices during the ten replicate runs. The assignment bar plots were generated by DISTRUCT 1.1 [ 65 ]. We employed the Circos visualization software [ 66 ] and exploited a circular ideogram to facilitate global distribution of the STRUCTURE prediction power detected in each deer population, arising from a comparison of the sika-deer microsatellite data. 2.4. Ethical Approval The entire sampling was performed post-mortem, involving only samples collected (independently of our research) from legally hunted red and sika deer in the Czech Republic. All applicable international, national, and/or institutional guidelines relating to the care and use of animals were followed. The procedures were carried out according to the Hunting Act (Law No. 449/2001), Decree No. 343/2015 Coll., on hunting periods for the individual game species and on detailed conditions governing hunting; for the given purpose, no specific ethical approval was required. 3. Results 3.1. Microsatellite Genetic Diversity and Population Differentiation All biparentally inherited microsatellite loci analyzed in the red and sika deer were polymorphic. A total of 176 distinct alleles were observed at eleven loci over the complete dataset (n = 123; the Cervus genus). The number of alleles at the individual loci ranged from 10 (BM1818) to 22 (BM888 and T193), with 16.00 alleles per locus on average (Table 1) . The H O ranged between 0.58 (RM188) and 0.84 (T193). The highest PIC was observed for the loci T156,T26,BM888,RM188,RT13, and T193; the values reached approximately 80% in the dataset of all individuals. The mean value of the major allele frequency across the loci corresponded to 0.276 in the dataset (Table 1). The number of alleles across all the samples classified into seven groups (Table 2) varied from 10.9 (Doupov red deer) to 4.27 (the Lány and Žlutice sika deer). Except for the RT13 (0.192) and RM188 (0.162) loci in the red deer, together with the BM1818 (0.174) locus in the sika deer, all values of the null allele estimated frequencies were smaller than 0.1 (null allele frequencies of ≥ 0.2 were
Animals 2021,11, 1701 6 of 16 considered large). The mean null allele frequency values equaled 0.07 and 0.04 across the 11 loci in the red and sika deer populations, respectively (Table 3). Table 1. The summary statistics of the microsatellite loci across all samples (n = 123). Marker MAF Genotype No. No. of Observations Allele No. Availability HEHOPIC F OarFCB5 0.4024 34 123 13 1.0000 0.7781 0.6992 0.7561 0.1055 T156 0.2125 55 120 17 0.9756 0.8996 0.8167 0.8920 0.0964 T26 0.1595 48 116 16 0.9431 0.9052 0.7586 0.8975 0.1661 BM888 0.2500 50 120 22 0.9756 0.8739 0.6917 0.8628 0.2125 RM188 0.2479 31 117 17 0.9512 0.8650 0.5812 0.8521 0.3319 RT1 0.2764 31 123 13 1.0000 0.8211 0.6423 0.7998 0.2217 T501 0.3696 34 115 14 0.9350 0.8050 0.6000 0.7868 0.2587 RT13 0.2059 43 119 17 0.9675 0.8737 0.5966 0.8612 0.3209 T193 0.2292 49 120 22 0.9756 0.9034 0.8417 0.8971 0.0725 BM1818 0.3824 30 119 10 0.9675 0.7943 0.6471 0.7745 0.1894 ETH225 0.2967 39 123 15 1.0000 0.8162 0.6585 0.7941 0.1971 Mean/Overall 0.2757 40.3636 119.5455 16 0.9719 0.8487 0.6849 0.8340 0.1971 MAF = major allele frequency; H E = expected heterozygosity; H O = observed heterozygosity; PIC = polymorphic information content; F= within-sampling site inbreeding coefficient. Table 2. The population summary statistics across the red and sika deer sampling sites. Population Code MAF Genotype No. nNo. of Observations Allele No. Availability HEHOPIC F Red deer Doupovské hory RD-D 0.2897 22.2727 38 36.8182 10.9091 0.9689 0.8276 0.7454 0.8096 0.1129 Red deer Slavkovskýles RD-S 0.2866 16.9091 25 24.1818 10.3636 0.9673 0.8242 0.7078 0.8057 0.1619 Sika deer Doupovské hory SD-D 0.4742 6.3636 10 9.7273 5.7273 0.9727 0.6621 0.6364 0.6240 0.0930 Sika deer Slavkovskýles SD-S 0.5323 7.3636 20 19.2727 5.2727 0.9636 0.6088 0.5564 0.5612 0.1126 Sika deer KladskáSD-K 0.4380 5.1818 9 8.8182 4.3636 0.9798 0.6806 0.8398 0.6254 − 0.1766 Sika deer Lány SD-L 0.5584 5.7273 13 12.7273 4.2727 0.9790 0.5850 0.5757 0.5345 0.0574 Sika deer Žlutice SD-Z 0.4659 4.9091 8 8.0000 4.2727 1.0000 0.6584 0.7159 0.6052 − 0.0208 MAF = major allele frequency; H E = expected heterozygosity; H O = observed heterozygosity; PIC = polymorphic information content; F= within-sampling-site inbreeding coefficient. Table 3. The locus-specific diversity measures estimated for the total sample sizes of the red (n = 63) and sika deer (n = 60). Marker AR HEHOPIC FNull Red Deer Sika Deer Red Deer Sika Deer Red Deer Sika Deer Red Deer Sika Deer Red Deer Sika Deer OarFCB5 11.690 6.982 0.8644 0.5986 0.8254 0.5667 0.8416 0.5658 0.0300 0.0357 T156 12.484 9.931 0.8846 0.8352 0.8065 0.8276 0.8650 0.8107 0.0493 0.0000 T26 12.789 10.724 0.8921 0.8505 0.8103 0.7069 0.8732 0.8240 0.0412 0.0832 BM888 16.653 6.991 0.8417 0.6915 0.7333 0.6500 0.8182 0.6514 0.0767 0.0120 RM188 9.675 8.000 0.7873 0.6731 0.5079 0.6667 0.7508 0.6397 0.1624 0.0179 RT1 11.390 5.899 0.6989 0.6423 0.6984 0.5833 0.6668 0.5980 0.0266 0.0474 T501 10.000 9.492 0.8111 0.5363 0.6607 0.5424 0.7859 0.5038 0.1126 0.0000 RT13 15.721 6.931 0.8824 0.7652 0.5410 0.6552 0.8637 0.7242 0.1917 0.0630 T193 14.771 6.900 0.9090 0.7182 0.9000 0.7833 0.8929 0.6767 0.0000 0.0000 BM1818 8.000 4.000 0.8291 0.6628 0.7742 0.5088 0.8037 0.5933 0.0086 0.1740 ETH225 13.810 4.791 0.8555 0.5167 0.7778 0.5333 0.8337 0.4120 0.0626 0.0000 Mean/Overall 12.453 7.331 0.8415 0.6809 0.7305 0.6386 0.8178 0.6363 0.0692 0.0394 AR = allelic richness; H E = expected heterozygosity; H O = observed heterozygosity; PIC = polymorphic information content; F Null = estimated frequencies of null alleles for each microsatellite locus.
Animals 2021,11, 1701 7 of 16 With regard to the average observed heterozygosity (Ho = 0.73) in the red deer, the values ranged lower than the overall average expected heterozygosity (H E = 0.84); consequently, the estimated inbreeding coefficient indicated a certain level of heterozygote deficiency. The total locus-specific gene diversity in the red deer (0.70–0.91) exceeded that established in the sika deer (0.52–0.85); the observed heterozygosities and PIC exhibited similar patterns (Table 3). The Wilcoxon signed-rank test suggested that the locus-specific diversity measures estimated for the total samples reached higher values in the red deer, attaining significant levels (H E and PIC at p< 0.01, p= 0.00169 and Ho at p< 0.05, p= 0.03754). The mean allelic diversity (AR with/out hybrids) ranged markedly lower in sikas (7.3/6.4) than in the red deer (12.5/11.9) across the loci (the Wilcoxon signed-rank test significant at p< 0.01). Populations (Cervus nippon nippon) with superior genetic qualities (higher Ho, H E , PIC, and lower f) were found in the subregions of Kladská, Žlutice, and Doupovské hory. The genetic variability at the Lány deer park and Slavkovskýles proved to be lower (Table 2). The presence of hybrid individuals in our dataset led to a slightly overestimated genetic diversity, mainly in the red deer populations (Table S1). Figure 2shows the NeighborNet dendrogram constructed from the Nei 0 s D A distances between the seven red and sika deer sampled populations. These groups tended to cluster together according to the geography. A similar pattern appeared when we analyzed the relationships between the sampling sites, also via a NeighborNet visualization of the Reynolds’ DRgenetic distances in terms of short-term evolution (Figure S1). Animals2021,11,x7of17 AR=allelicrichness;HE=expectedheterozygosity;HO=observedheterozygosity;PIC=polymorphicinformationcontent; FNull=estimatedfrequenciesofnullallelesforeachmicrosatellitelocus. Withregardtotheaverageobservedheterozygosity(Ho=0.73)inthereddeer,the valuesrangedlowerthantheoverallaverageexpectedheterozygosity(HE=0.84);conse‐ quently,theestimatedinbreedingcoefficientindicatedacertainlevelofheterozygotede‐ ficiency.Thetotallocus‐specificgenediversityinthereddeer(0.70–0.91)exceededthat establishedinthesikadeer(0.52–0.85);theobservedheterozygositiesandPICexhibited similarpatterns(Table3).TheWilcoxonsigned‐ranktestsuggestedthatthelocus‐specific diversitymeasuresestimatedforthetotalsamplesreachedhighervaluesinthereddeer, attainingsignificantlevels(HEandPICatp<0.01,p=0.00169andHoatp<0.05,p= 0.03754).Themeanallelicdiversity(ARwith/outhybrids)rangedmarkedlylowerinsikas (7.3/6.4)thaninthereddeer(12.5/11.9)acrosstheloci(theWilcoxonsigned‐ranktestsig‐ nificantatp<0.01).Populations(Cervusnipponnippon)withsuperiorgeneticqualities (higherHo,HE,PIC,andlowerf)werefoundinthesubregionsofKladská,Žlutice,and Doupovskéhory.ThegeneticvariabilityattheLánydeerparkandSlavkovskýlesproved tobelower(Table2).Thepresenceofhybridindividualsinourdatasetledtoaslightly overestimatedgeneticdiversity,mainlyinthereddeerpopulations(TableS1). Figure2showstheNeighborNetdendrogramconstructedfromtheNei′sDAdis‐ tancesbetweenthesevenredandsikadeersampledpopulations.Thesegroupstendedto clustertogetheraccordingtothegeography.Asimilarpatternappearedwhenweana‐ lyzedtherelationshipsbetweenthesamplingsites,alsoviaaNeighborNetvisualization oftheReynolds’DRgeneticdistancesintermsofshort‐termevolution(FigureS1). Figure2.TheNeighborNetdendrogramconstructedfromtheNei´sDAdistancesbetweenthesevenredandsikadeer samplingpopulations(n=123).Populationcodes:RD‐D(ReddeerDoupovskéhory),RD‐S(ReddeerSlavkovskýles),SD‐ D(SikadeerDoupovskéhory),SD‐S(SikadeerSlavkovskýles),SD‐K(SikadeerKladská),SD‐L(SikadeerLány),andSD‐ Z(SikadeerŽlutice). Asignificantgenotypiclinkagedisequilibriumwasdetectedinonlyonepairofloci (T501×RT13)inbothspecies;thisfactthensupportsthehypothesisthatthegivenloci segregateindependentlyineachspecies’genome(TableS2). WithregardtotheWright’sF‐statistics,theoverallFIS,FIT,andFSTvaluesequaled 0.103,0.273,and0.189(p<0.05),respectively.Concerningonlythenon‐hybridanimals, thevalueswereasfollows:FIS=0.098,FIT=0.285,andFST=0.207(p<0.05).However,FST mightproveinappropriateforcomparinglociwithsubstantiallydifferentlevelsofvaria‐ tion,couldbemisleading,andmayyieldwrongresultsinrecentlyisolatedpopulations (butstillmaycontainsomesimilarityduetocommonancestralpopulation).Forthisrea‐ son,wecalculatedthegeneticdistance,whichshowedmuchhighervalues(TableS3). Figure 2. The NeighborNet dendrogram constructed from the Nei ´ s DA distances between the seven red and sika deer sampling populations (n = 123). Population codes: RD-D (Red deer Doupovskéhory), RD-S (Red deer Slavkovskýles), SD-D (Sika deer Doupovskéhory), SD-S (Sika deer Slavkovskýles), SD-K (Sika deer Kladská), SD-L (Sika deer Lány), and SD-Z (Sika deer Žlutice). A significant genotypic linkage disequilibrium was detected in only one pair of loci (T501 × RT13) in both species; this fact then supports the hypothesis that the given loci segregate independently in each species’ genome (Table S2). With regard to the Wright’s F-statistics, the overall F IS ,F IT , and F ST values equaled 0.103, 0.273, and 0.189 (p< 0.05), respectively. Concerning only the non-hybrid animals, the values were as follows: F IS = 0.098, F IT = 0.285, and F ST = 0.207 (p< 0.05). However, F ST might prove inappropriate for comparing loci with substantially different levels of variation, could be misleading, and may yield wrong results in recently isolated populations (but still may contain some similarity due to common ancestral population). For this reason, we calculated the genetic distance, which showed much higher values (Table S3). 3.2. Clustering Analysis and Interspecific Hybridization The genetic structures of the investigated samples were evaluated by using Bayesian model-based clustering in the STRUCTURE software. The STRUCTURE analysis of the
Animals 2021,11, 1701 8 of 16 microsatellite genotypes in both species separated according to the sampling sites provided the strongest support for the grouping of the genetic variation into two clusters (K= 2) based on ∆ K= 1393.14 (Figure 3a). The potential of the sika deer genetic introgression into the red deer and vice versa was revealed; the individuals were assigned to the sika deer cluster (at the membership probability of Q ≤ 0.01) and the red deer cluster (Q ≥ 0.99). Here, a “nuclear hybrid” was defined via the nuclear/microsatellite markers as an individual returning the Qvalue of 0.01 < Q< 0.99 between two taxa, with the confidence intervals calculated at 95%, as described in the Materials and Methods Section. Using the definitions outlined above, we found 56 pure sika deer (i.e., those with a CI overlapping 0.01), achieving an average Qscore of 0.0025 ± 0.0015 ( ± SD); 54 pure red deer (CI overlapping 0.99), obtaining an average Qscore of 0.9968 ± 0.0033; and 13 hybrid deer, where the Q scores ranged between 0.0193 and 0.9896. In the sika deer populations, we identified four hybrid animals at the membership probabilities of 1.93%, 2.78%, 6.78%, and 31.01%; the red deer populations then yielded nine hybrid individuals (at the membership probabilities of 98.96–94.82%) (Figure 3b). The total hybrid score equaled 10.57% (13/123), with 14.29% (9/63) of the hybrids being among red deer–like individuals and 6.67% (4/60) among sika-like ones. In the sikas from the subregions of Kladská, Žlutice, and Lány, no hybrid animals were identified. In Doupovskéhory, the total interspecific hybrid score reached 18.75% (9/48); 15.79% of the hybrids ranged among red deer–like individuals (6/38), and 30% fell within the sika-like ones (3/10). In Slavkovskýles, the total hybrid score amounted to 8.89% (4/45) with 12% of red-like hybrids (3/25) and 5% of sika-like hybrids (1/20). Animals2021,11,x8of17 3.2.ClusteringAnalysisandInterspecificHybridization ThegeneticstructuresoftheinvestigatedsampleswereevaluatedbyusingBayesian model‐basedclusteringintheSTRUCTUREsoftware.TheSTRUCTUREanalysisofthe microsatellitegenotypesinbothspeciesseparatedaccordingtothesamplingsitespro‐ videdthestrongestsupportforthegroupingofthegeneticvariationintotwoclusters(K =2)basedon∆K=1393.14(Figure3a).Thepotentialofthesikadeergeneticintrogression intothereddeerandviceversawasrevealed;theindividualswereassignedtothesika deercluster(atthemembershipprobabilityofQ≤0.01)andthereddeercluster(Q≥0.99). Here,a“nuclearhybrid”wasdefinedviathenuclear/microsatellitemarkersasanindi‐ vidualreturningtheQvalueof0.01<Q<0.99betweentwotaxa,withtheconfidence intervalscalculatedat95%,asdescribedintheMaterialsandMethodssection.Usingthe definitionsoutlinedabove,wefound56puresikadeer(i.e.,thosewithaCIoverlapping 0.01),achievinganaverageQscoreof0.0025±0.0015(±SD);54purereddeer(CIoverlap‐ ping0.99),obtaininganaverageQscoreof0.9968±0.0033;and13hybriddeer,wherethe Qscoresrangedbetween0.0193and0.9896.Inthesikadeerpopulations,weidentified fourhybridanimalsatthemembershipprobabilitiesof1.93%,2.78%,6.78%,and31.01%; thereddeerpopulationsthenyieldedninehybridindividuals(atthemembershipprob‐ abilitiesof98.96–94.82%)(Figure3b). (a) Figure 3. Cont.
Animals 2021,11, 1701 9 of 16 Animals2021,11,x9of17 (b) Figure3.TheBayesianmodel‐basedclusteringrenderedwiththeSTRUCTUREsoftware:(a)theevolutionofthemeanln oflikelihood(lnP(K))accordingtoK,basedonthefiverunsofthe50,000burn–insand150,000MCMCs(standarddevia‐ tionsindicated);(b)STRUCTUREclusteringresultsatK=2(theadmixturemodeland11‐locusdataset).Eachverticalline representsoneindividual.Thethinblacklinesseparateindividualsfromdifferentsamplingsites(groups).Population codes:RD‐D(ReddeerDoupovskéhory),RD‐S(ReddeerSlavkovskýles),SD‐D(SikadeerDoupovskéhory),SD‐S(Sika deerSlavkovskýles),SD‐K(SikadeerKladská),SD‐L(SikadeerLány),andSD‐Z(SikadeerŽlutice). Thetotalhybridscoreequaled10.57%(13/123),with14.29%(9/63)ofthehybridsbe‐ ingamongreddeer–likeindividualsand6.67%(4/60)amongsika‐likeones.Inthesikas fromthesubregionsofKladská,Žlutice,andLány,nohybridanimalswereidentified.In Doupovskéhory,thetotalinterspecifichybridscorereached18.75%(9/48);15.79%ofthe hybridsrangedamongreddeer–likeindividuals(6/38),and30%fellwithinthesika‐like ones(3/10).InSlavkovskýles,thetotalhybridscoreamountedto8.89%(4/45)with12% ofred‐likehybrids(3/25)and5%ofsika‐likehybrids(1/20). FromtheconfusionmatrixoftheBayesianmethodascalculatedbySTRUCTUREfor eachred‐sikadeersamplingsite,weconstructedacircularideogram(Figure4)displaying variationinthegenomestructures.Figure5indicatesthatthesetof123individualsphe‐ notypicallydesignatedasred(n=63)andsikadeer(n=60)throughgeneticidentification comprised54animalshavingtheQ‐valuesof0.99–1.0(“purereddeer”),nineanimalsex‐ hibiting0.75–0.99(hybridofthe“redtype”),oneindividualshowing0.25≤Q≤0.75(an intermediatehybrid),threeindividualswith0.01–0.25(ahybridofthe“sikatype”),and 56animalscharacterizedby0–0.01(“puresika”).Afterremovingallthered‐sikahybrid animals,thedatasetcontaining110animalswasre‐analyzedinSTRUCTUREandclus‐ teredintotwogroups;theassignmenttestclearlyisolatedtheredandsikapopulations withouttheadmixedtypes. Figure 3. The Bayesian model-based clustering rendered with the STRUCTURE software: ( a ) the evolution of the mean ln of likelihood (ln P(K)) according to K, based on the five runs of the 50,000 burn–ins and 150,000 MCMCs (standard deviations indicated); ( b ) STRUCTURE clustering results at K= 2 (the admixture model and 11-locus dataset). Each vertical line represents one individual. The thin black lines separate individuals from different sampling sites (groups). Population codes: RD-D (Red deer Doupovskéhory), RD-S (Red deer Slavkovskýles), SD-D (Sika deer Doupovskéhory), SD-S (Sika deer Slavkovskýles), SD-K (Sika deer Kladská), SD-L (Sika deer Lány), and SD-Z (Sika deer Žlutice). From the confusion matrix of the Bayesian method as calculated by STRUCTURE for each red-sika deer sampling site, we constructed a circular ideogram (Figure 4) displaying variation in the genome structures. Figure 5indicates that the set of 123 individuals phenotypically designated as red (n = 63) and sika deer (n = 60) through genetic identification comprised 54 animals having the Q-values of 0.99–1.0 (“pure red deer”), nine animals exhibiting 0.75–0.99 (hybrid of the “red type”), one individual showing 0.25 ≤ Q ≤ 0.75 (an intermediate hybrid), three individuals with 0.01–0.25 (a hybrid of the “sika type”), and 56 animals characterized by 0–0.01 (“pure sika”). After removing all the red-sika hybrid animals, the dataset containing 110 animals was re-analyzed in STRUCTURE and clustered into two groups; the assignment test clearly isolated the red and sika populations without the admixed types. Animals2021,11,x10of17 Figure4.TheCircos‐likeplotdisplayingtheinterspeciesvariationbasedonthenDNAdata.Popu‐ lationcodes:SD‐K(SikadeerKladská),SD‐L(SikadeerLány),andSD‐Z(SikadeerŽlutice). Figure5.Theestimatedproportionofinterspecifichybridizationancestry(Q)betweenthesikaand reddeer(n=123),andthethresholdvaluesof0.01<Q<0.99usedtodetectthehybridindividuals inthewesternpartsoftheCzechRepublic.Wedetectedmorered‐likehybrids(n=9)exhibitinga recentsikaancestry(0.50<Q<0.99)thansika‐likeones(n=4)havingarecentreddeerancestry (0.01<Q<0.50). 4.Discussion Thehybridizationbetweentheredandsikadeerwasdocumentedbyusingavailable sourcesincludingtheolderphenotypicandcraniological[18,46]andtherecentgenetic [22]analyses;thesereferencesmakeitobviousthatmerelyonegeneticreportconcerning red‐sikadeerhybridizationintheCzechRepublicwascompletedpreviously.Ourpaper Figure 4. The Circos-like plot displaying the interspecies variation based on the nDNA data. Population codes: SD-K (Sika deer Kladská), SD-L (Sika deer Lány), and SD-Z (Sika deer Žlutice).
Animals 2021,11, 1701 16 of 16 46. Bartoš, L.; Hyánek, J.; Žirovnický, J. Hybridization between red and sika deer (Cervus nippon) in Czechoslovakia. Folia Zool. Brno 1981,31, 195–208. 47. Bartoš, L.; Žirovnický, J. Hybridization between red and sika deer, III. Interspecific behaviour. Zool. Anz. 1982,208, 30–36. 48. Powerscourt, V. On the acclimatization of the Japanese deer at Powerscourt. In Proceedings of the Zoological Society of London, London, UK, 1 March 1884; pp. 207–209. [CrossRef] 49. Bishop, M.D.; Kappes, S.M.; Keele, J.W.; Stone, R.T.; Sunden, S.L.; Hawkins, G.A.; Toldo, S.S.; Fries, R.; Grosz, M.D.; Yoo, J.; et al. A genetic linkage map for cattle. Genetics 1994,136, 619–639. [CrossRef] [PubMed] 50. Barendse, W.; Armitage, S.M.; Kossarek, L.M.; Shalom, A.; Kirkpatrick, B.W.; Ryan, A.M.; Clayton, D.; Li, L.; Neibergs, H.L.; Zhang, N.; et al. A genetic linkage map of the bovine genome. Nat. Genet. 1994,6, 227–235. [CrossRef] 51. Steffen, P.; Eggen, A.; Dietz, A.B.; Womack, J.E.; Stranzinger, G.; Fies, R. Isolation and mapping of polymorphic microsatellites in cattle. Anim. Genet. 1993,24, 121–124. [CrossRef] [PubMed] 52. Buchanan, F.C.; Galloway, S.M.; Crawford, A.M. Ovine microsatellites at the OarFCB5, OarFCB19, OarFCB20, OarFCB48, OarFCB129 and OarFCB226 loci. Anim. Genet. 1994,25, 60. [CrossRef] 53. Wilson, G.A.; Strobeck, C.; Wu, L.; Coffin, J.W. Characterization of microsatellite loci in caribou Rangifer tarandus, and their use in other artiodactyls. Mol. Ecol. 1997,6, 697–699. [CrossRef] 54. Jones, K.C.; Levine, K.F.; Banks, J.D. Characterization of 11 polymorphic tetranucleotide microsatellites for forensic applications in California elk (Cervus elaphus canadensis). Mol. Ecol. Resour. 2002,2, 425–427. [CrossRef] 55. Liu, K.; Muse, S.V. PowerMarker: Integrated analysis environment for genetic marker data. Bioinformatics 2005 ,21, 2128–2129. [CrossRef] 56. FSTAT. Available online: http://www2.unil.ch/popgen/softwares/fstat.htm (accessed on 17 March 2021). 57. Rousset, F. Genepop’007: A complete reimplementation of the Genepop software for Windows and Linux. Mol. Ecol. Resour. 2008 , 8, 103–106. [CrossRef] 58. Nei, M.; Tajima, F.; Tateno, Y. Accuracy of estimated phylogenetic trees from molecular data. J. Mol. Evol. 1983 ,19, 153–170. [CrossRef] 59. Reynolds, J.; Weir, B.S.; Cockerham, C.C. Estimation of the Coancestry coefficient: Basic for a short-term genetic distance. Genetics 1983,105, 767–779. [CrossRef] 60. Huson, D.H.; Bryant, D. Application of Phylogenetic Networks in Evolutionary Studies. Mol. Biol. Evol. 2006 ,23, 254–267. [CrossRef] [PubMed] 61. Pritchard, J.K.; Stephens, M.; Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 2000 ,155, 945–959. [CrossRef] 62. Earl, D.A.; von Holdt, B.M. STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 2012,4, 359–361. [CrossRef] 63. Evanno, G.; Regnaut, S.; Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 2005,14, 2611–2620. [CrossRef] [PubMed] 64. Jakobsson, M.; Rosenberg, N.A. CLUMPP: A cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 2007,23, 1801–1806. [CrossRef] 65. Rosenberg, N.A. Distruct: A program for the graphical display of population structure. Mol. Ecol. Notes 2004 ,4, 137–138. [CrossRef] 66. Krzywinski, M.; Schein, J.; Birol, I.; Connors, J.; Gascoyne, R.; Horsman, D.; Jones, S.J.; Marra, M.A. Circos: An information aesthetic for comparative genomics. Genome Res. 2009,19, 1639–1645. [CrossRef] [PubMed] 67. Slate, J.; Coltman, D.W.; Goodman, S.J.; MacLean, I.; Pemberton, J.M.; Williams, J.L. Bovine microsatellite loci are highly conserved in red deer (Cervus elaphus), sika deer (Cervus nippon) and Soay sheep (Ovis aries). Anim. Genet. 1998 ,29, 307–315. [CrossRef] [PubMed] 68. Krojerová-Prokešová, J.; Baranˇceková, M.; Kawata, Y.; Oshida, T.; Igota, H.; Koubek, P. Genetic differentiation between introduced Central European sika and source populations in Japan: Effects of isolation and demographic events. Biol. Invasions 2017 ,19, 2125–2141. [CrossRef] 69. Ziegrosser, P. Sika—Nep˚uvodní—Invazivní—Druh. Myslivost 2017,12, 18. 70. Pipek, P. Seek Sika. Invasion of Sika Deer in the Czech Republic and Europe. Živa 2018,5, 280–281. 71. Beysard, M.; Heckel, G. Structure and dynamics of hybrid zones at different stages of speciation in the common vole (Microtus arvalis). Mol. Ecol. 2014,23, 673–687. [CrossRef] [PubMed] 72. Macháˇcek, Z.; Dvoˇrák, S.; Ježek, M.; Zahradník, D. Impact of interspecific relations between native red deer (Cervus elaphus) and introduced sika deer (Cervus nippon) on their rutting season in the Doupovskéhory Mts. J. For. Sci. 2014 ,60, 272–280. [CrossRef] 73. Wyman, M.; Locatelli, Y.; Charlton, B.; Reby, D. Female Sexual Preferences Toward Conspecific and Hybrid Male Mating Calls in Two Species of Polygynous Deer, Cervus elaphus and C. nippon.Evol. Biol. 2016,43, 227–241. [CrossRef] 74. Li, Z.; Wright, A.G.; Si, H.; Wang, X.; Qian, W.; Zhang, Z.; Li, G. Changes in the rumen microbiome and metabolites reveal the effect of host genetics on hybrid crosses. Environ. Microbiol. Rep. 2016,8, 1016–1023. [CrossRef] [PubMed] 75. Vähä, J.P.; Primmer, C.R. Efficiency of model-based Bayesian methods for detecting hybrid individuals under different hybridization scenarios and with different numbers of loci. Mol. Ecol. 2006,15, 63–72. [CrossRef] [PubMed]