Genetic insights into an Apennine population of the Italian red deer

The red deer Cervus elephus has been a common species in Italy until the Middle Ages and the Renaissance, when its distribution range started to considerably decrease, due to gradual deforestation and hunting pressure. Afterwards, the red deer has been reintroduced to many regions of the world, including Italy. In the Italian Apennines, the Acquerino-Cantagallo Natural Reserve (ACQUERINO) hosts one of the largest peninsular red deer populations, originated from a series of successful reintroductions. In this study, we meant to detect the level of genetic variability of Acquerino-Cantagallo Natural Reserve deer population and to investigate the genetic relationships with the other Italian and European populations. We identified five mitochondrial DNA control region (D-loop) haplotypes, four falling in lineage A and one falling in lineage C, derived from at least two maternal lineages, confirming that ACQUERINO population should be the result of multiple reintroductions. Haplotype diversity (H = 0.50) and nucleotide (π = 0.004) diversity were low, but included into the deer range values. ACQUERINO population showed low levels of genetic diversity when compared to other European and Mediterranean populations, confirming that this expanding population may have been generated from a low number of founders.


Introduction
Ungulates (Artiodactyla) are a group of mammals which importance in terrestrial zoocenoses is paramount (Putman 1996;Reimoser 2003;Murray et al. 2013). They are often the largest species in size, thus arousing a strong hunting interest and becoming a substantial management commitment, especially in the last decades (Hobbs 1996;Reimoser 2003). Among ungulates, the European red deer Cervus elaphus Linnaeus, 1758 is not an endangered species, but its populations have suffered from increased isolation processes and frequent reintroductions by allochthonous stocks (Hartl et al. 2003;Frantz et al. 2017). Accordingly, introgression of genes of non-native or captive origin red deer into native populations is relatively frequent throughout all of Europe, as a consequence of translocations (Iacolina et al. 2019).
Until the tenth/eleventh centuries, the red deer was distributed over most of Europe, and it was still a common species in almost all the Italian mainland (Mattioli 1990). During the late Middle Ages and the Renaissance, its distribution range decreased considerably, due to gradual deforestation and hunting pressure (Lovari 1993). With the further improvement of firearms, between the eighteenth and nineteenth centuries, red deer populations suffered a rapid decrease, so that, at the beginning of the 1900s, the red deer was nearly extinct over almost the entire Italian territory, except for the Mesola Natural Reserve (Northern-Central Italy, Ferrara District) and Sardinia (Mattioli 1990;Mattioli et al. 2003;Zachos et al. 2009Zachos et al. , 2014Gippoliti et al. 2020). In peninsular Italy, the only native population is currently represented by the Mesola population, on the eastern coastline. At present, in Italy, several populations are the result of a series of reintroductions between 1950 and 2003, followed by a natural expansion in the Alps from neighbouring European areas (Switzerland, Austria, and Slovenia: Mattioli et al. 2001;Carnevali et al. 2009) (Fig. S1). Other small nuclei in the Apennines are indeed the results of local unmanaged escapes (e.g. Monte Penna in Castell'Azzara-Grosseto and Merse river valley in the province of Siena). In 2010, the Italian territory included over 70,000 individuals divided as follows: Central-Western Alps 20,000, Eastern Alps 32,000, Northern Apennines 8,000, Central-Southern Apennines 2,000, Islands 8,000 (cf. Apollonio et al. 2010). Since then, no reliable estimation of red deer population has been conducted in Italy. For a species with a potential annual increase of up to 30%, current numbers could be quite different compared to those collected 12 years ago, especially for populations of Apennine protected areas (see Chirichella et al. 2017). In the Italian Alps, the deer shows a favourable conservation status and has occupied a large part of the potential range again, so much as in some geographical regions, the populations must be controlled in order to avoid excessive forest damages. Red deer populations in Northern Apennines are currently distributed in three main continuous areas (Becciolini et al. 2016;Gippoliti et al. 2020): in the north-west, with a population estimated to be nearly 3,100 individuals; in the south-eastern area, with a population of 2,500-3,000 deer; in the central area, with about 3,000-3,500 individuals. The progressive expansion of distribution area of the deer population presents in the Tosco-Emilian Apennines led in 1999 to the creation of a protocol of agreement for the management of the species at the inter-regional level, named ACATER (Areale del Cervo dell'Appennino Tosco-Emiliano e Romagnolo), which aim is to guarantee a correct management of the red deer in the long term, taking into account the biological and behavioural characteristics of the species. In Central and Southern Apennines, the perspective for a natural expansion of populations seems to be highly probable, given the extent of protected areas in the Abruzzo and Lazio regions. In addition, several areas of the Southern Apennines have environmental conditions suitable for hosting this species and may be interested by future reintroductions.
Mitochondrial markers (for a review see Zachos and Hartl 2011) identified in European red deer modern populations three major haplogroups including the several subspecies described for this taxon (Skog et al. 2009 in Italy. C. e. italicus (formerly in haplogroup C) (Ludt et al. 2004;Skog et al. 2009) and haplogroup E occur in Armenia, Crimea, and Southern Russia (Doan et al. 2022). The above distinct mitochondrial DNA (mtDNA) lineages would reflect the main refugia during the Last Glacial Maximum (LGM). From Iberia, lineage A colonized the British Isles and parts of Central Europe; in the Balkan region, lineage C migrated to South-East Central Europe and the Carpathians. The Mediterranean (lineage B) could be originated from the Italian refugial lineage, and lineage D has been proposed for the Middle East and South-Eastern Europe (Zachos and Hartl 2011;Meiri et al. 2013;Doan et al. 2017;Schnitzler et al. 2018;Queiros et al. 2019). These hypotheses have been recently reconsidered and new data suggested that the species survived the LGM not only in the well-known southern European refugia, but also in more northern areas of Western and Eastern Europe and in the Urals (Niedziałkowska et al. 2021).
The genetic structure of the red deer in Europe has been widely investigated (Ludt et al. 2004;Skog et al. 2009;Perez-Espona et al. 2009;Niedziałkowska et al. 2011;Fernández-García et al. 2014;Zachos et al. 2016;Rey-Iglesia et al. 2017;Schnitzler et al. 2018;Doan et al. 2022) even though a small percentage of data involved peninsular Italian populations (Lorenzini et al. 2005;Hmwe et al. 2006A, BA;Zachos et al. 2016;Doan et al. 2017Doan et al. , 2022. Doan et al. (2022) also highlighted Italy as one of the countries with the lowest genetic amount of data as to red deer populations. The ACQUERINO population is one of the largest populations in Central Italy, and its morphometry is well studied (Becciolini et al. 2016). Even if originated by reintroductions, we considered of interest a genetic characterization of this population that has experienced strong demographic expansion. The ACQUERINO population was originated by seven red deer individuals reintroduced in 1958 and 1965 from the State Forest of Tarvisio, North-Eastern Italy (Mazzarone and Mattioli 1996) in the Acquerino-Cantagallo Natural Reserve followed by some not documented introductions.
In this paper, we aimed at evaluating the level of genetic variability of Acquerino-Cantagallo Natural Reserve deer population and to investigate the genetic relationships with the other Italian and European populations.

Study area
We conducted our study in the Acquerino-Cantagallo Natural Reserve (provinces of Prato and Pistoia, Central Italy, 44.003°N-11.015°E: 100-1100 m a.s.l., hereafter ACQUERINO), an historical reproductive area of the Northern Apennine red deer population (Becciolini et al. 2016). The local climate shows sub-montane characteristics. Average monthly temperatures range from 0 to 6 °C in the coldest months (October-March) and from 17 to 24 °C in the warmest ones (April-September). Annual precipitations are generally over 1,000 mm, with seasonal maximum in autumn and spring, and the minimum in summer (Becciolini et al. 2016). Forests cover the 84% of the study area, with broad-leaved forests being the most abundant type (73%), followed by mixed broad-leaved and conifer (7%) and conifer woodlands (4%). Farmlands and olive groves represent 3% and 2% of the territory, respectively (Becciolini et al. 2016).

Sampling and DNA extraction
We collected 26 blood red deer samples (13 males and 13 females) in 2008-2009 as part of a selective hunting program. Samples were stored in EDTA 0.25 M and kept in a freezer at − 80 °C. DNA was extracted using the Qiagen Qiamp Tissue and Blood kit (Qiagen, Inc, Tokyo, Japan). The analysis of genetic variability was carried out using a portion (738 bp) of the mitochondrial DNA control region (mt CR, D-loop). PCR amplification was obtained using two pairs of primers already available in the literature: CE-CR-FOR (5′-CAA TAC ACT GGT CTT GTA ACC-3′) and CE-CR-REV (5′-TAA TAG GAA GGG GGGCC-3′; McDevitt et al. 2009). The PCR was carried out with an Applied Biosystems 2700 machine in 25 μl comprising buffer 10 × , 12 mM MgCl 2 , 200 μM dNTPs, 0 2 μM of each primer, and 1 unit Taq polymerase (Life Technologies). The amplification conditions included initial denaturing at 94 °C for 2 min, followed by 35 cycles of 94 °C for 45″, 55 °C (T annealing) for 30″, 72 °C for 1 min, and a final extension to 72 °C for 10 min PCR products were run by electrophoresis on 15% agarose gel, containing 05 mg/ml of SYBR green gel staining, purified (Sure Clean method, Bioline), and sequenced with the ABI Big Dye™ Terminator Cycle Sequencing version 20-ABI PRISM (Applied Biosystems, Foster City, USA) kit. The sequence analyses were carried out using an automatic ABI 310 sequencer (Applied Biosystems). Electropherograms have been displayed with the CHROMAS 145 program (http:// www. techn elysi um. com/ au). The sequences have been manually corrected and analyzed with the MEGA 6 program (Tamura et al. 2013).

Genetic and statistical analyses
We used DNA-sp version 5 (Librado and Rozas 2009) to calculate the number of haplotypes (Nh), haplotypic diversity (H), and nucleotide diversity (π). We used DNA-sp also to investigate the demographic history, calculating deviations from neutrality with Tajima's D (Tajima 1989), Fu and Li's F (Fu and Li 1993), and Fu's Fs (Fu 1997) tests. Significant D, F, and Fs negative values indicated an excess of low-frequency polymorphisms relative to expectation, indicating population size expansion (following a bottleneck or a selective sweep) or purifying selection. Fu's Fs is considered a more powerful tool than Tajima's D. Its power comes from the ability to differentiate between population growth and genetic hitchhiking and in rejecting the hypothesis of neutral mutations. JModelTEST 304 (Posada and Crandall 1998) was used to test the most accurate model of substitution using the Bayesian information criterion (BIC; Schwarz 1978) and Akaike's information criterion (AIC; Akaike 1974), corrected for the heterogeneity between sites (gamma [G]). Red deer sequences obtained in this study were aligned with previously published D-loop sequences available on GenBank (http:// www. ncbi. nlm. nih. gov/) (datasets 1 and 2, Table 1). A first dataset (dataset 1) consisted of longer sequences of about 700 bp and a second dataset (dataset 2) consisted of shorter sequences (about 320 bp) but included more European samples (Table 1). A median-joining network was constructed to determine the genealogical relationships and haplotype relative frequencies for ACQ population (Bandelt et al. 1999) by software PopART (Leigh and Bryant 2015). We applied a Bayesian approach performed by MrBayes v 3.2.6. (http:// mrbay es. sourc eforge. net/ downl oad. php). We ran four Markov chains in two independent analyses for 10 million generations, sampling every 1,000 generations. The first 25% of samples were discarded as burn-in. The sequences of C. canadensis and C. nippon (Table 1) were used as outgroups. Support node values were given as posterior probability values (pp).
A likelihood-ratio test (LRT) was performed to test the molecular clock hypothesis based on D-loop gene sequences. The likelihoods with the molecular clock assumption (L0) and without the molecular clock assumption (L1) were calculated using MEGA software. Since the clock-like behaviour of the data was rejected for the complete dataset (the equal evolutionary rate throughout the tree was rejected at a 5% significance level (P = 1.181E-030)), a relaxed clock (uncorrelated lognormal) was applied, as implemented in BEAST 1.5.3 using a Bayesian MCMC approach (Drummond and Rambaut 2007). We used a CR gene substitution rate of 3,849 × 10 -7 (per site per year, Queirós et al. 2019). Two separate MCMC analyses were run for 20 million generations for each dataset with parameters sampled every 1,000 generations. Independent runs were combined using LogCombiner, and the first 20% of the generations from each run were discarded as "burnin". The convergence of the chains was checked using TRACER. The searches achieved adequate mixing as assessed by the high effective sampling size (ESS) values for all parameters of 100 or greater. Node ages and upper and lower bounds of the 95% highest posterior density interval for divergence times were calculated using TreeAnnotator and visualized using FigTree 1.31 (Rambaut 2009).

Results
We detected five new haplotypes among the 26 examined samples: HAP1_ACQUERINO, HAP2_ACQUERINO, HAP3_ACQUERINO, HAP4_ACQUERINO, HAP5_ ACQUERINO (GenBank Acc. No: MZ997342-MZ997346): the HAP1_ACQUERINO is the most frequent one (18 individuals), followed by HAP2_ACQUERINO (5 individuals). The other three haplotypes were found in one individual each (Fig. 1). In the median-joining network, the HAP5_ACQUERINO haplotype is included within Both trees (obtained with dataset 1 and dataset 2; tree with dataset 2 not shown) gave very similar results (Fig. 2): all ACQUERINO haplotypes belong to the  Skog et al. 2009). The red deer haplogroup A, including most ACQUERINO haplotypes, was well supported statistically ( Fig. 2, P = 0.98), unlike haplogroup C that showed a low support value (P = 0.62). All neutrality tests showed a demographic expansion with negative and significant values (Tajima D: − 243, P-value < 0.001; Fu and Li's D: − 401, P-value < 0.002; Fu's Fs: − 413 P-value < 0.002). Times of divergence between the three haplogroups reported in the chronograms (Fig. S2) showed as the European red deer population diverged during the last 3-1 million years.

Discussion
The low level of genetic variability found in ACQUERINO population is related to the small number of source individuals (Tarvisio Forest, NE Italy), along with the unknown multiple reintroductions in the 1980s by some individuals from farms located in the Bologna District (Mazzarone and Mattioli 1996). The ACQUERINO population showed low levels of genetic diversity when compared to those of various European and Mediterranean populations (Doan et al. 2017), whereas being similar to values obtained by Queiros et al.
. Besides, also active migrations may have occurred from close Alpine/Apennine populations since the 1960s throughout the Italian peninsula (Apollonio et al. 2010). However, the increased population size of ACQUERINO population, along with a constant increasing trend of Italian red deer population, is in agreement with the demographic tests militating for a population expansion. Another outcome of mitochondrial analysis, although preliminary, suggested a mixed composition of the ACQUERINO population with the most part of haplotypes grouped in lineage A, identified by Skog et al. (2009). The only one exception is haplotype 5 that grouped in lineage C including Tarvisio, Bulgaria, Hungary, and Poland, the last one representing the closest sequence to HAP5_ACQUERINO. Therefore, although the original population in ACQUERINO derived from free-ranging individuals captured in North-Eastern Italy (Mazzarone and Mattioli 1996), several unofficial releases from local caged areas may have occurred in the last 50 years, thus contributing to this genetic diversity. Genetic studies on the red deer have involved mitochondrial DNA sequences (Ludt et al. 2004;Skog et al. 2009;Perez-Espona et al. 2009;Niedziałkowska et al. 2011;Fernández-García et al. 2014;Borowski et al. 2016;Schnitzler et al. 2018;Queiròs et al. 2019;Doan et al. 2022) andmicrosatellites (Hwme et al. 2006B;Zachos et al. 2016). These authors have shown that the large-scale genetic structure of European red deer has been shaped by the Late Pleistocene and Holocene glacial-interglacial cycles. The role of Italy as a possible refugium and the origin of North African and Sardinian red deer is likely. During the Quaternary ice ages, most of Northern Europe was covered by ice. Conversely, Southern Europe, including Italy, provided ice-free glacial refugia during the Late Pleistocene for many temperate species (Meiri et al. 2013;Carranza et al. 2016). Main European glacial refugia corresponded to the southern peninsulas of Iberia, Italy, and the Balkans. These three main refugia acted as bases for the postglacial recolonization of the continent for different taxa, carrying distinct genetic lineages as an inheritance of their long isolation from each other (Meiri et al. 2013;Carranza et al. 2016). Niedziałkowska et al. (2021) have recently suggested that the species survived the Last Glacial Maximum not only in the abovementioned southern European refugia, but also in more northern areas of Western and Eastern Europe and in the Urals. The history of haplogroups took place in a period of 3-1 mya, i.e. dating older than in previous data (Doan et al. 2022), but in agreement with other studies using the same molecular marker (Randi et al. 2001).
Even if our results are based on only 26 samples and cannot yet be considered conclusive, our data represent the first outcomes regarding the genetics of a peninsular Italian red deer population, except restriction fragment analysis (Lorenzini et al. 2005) and the Mesola Forest Reserve population (Zachos et al. 2014). The phylogenetic analysis provided information on the composition and genetic variability of the ACQUERINO deer population, confirming that the current distribution of red deer in Italy would be the result of a series of successful (re) introductions. Accordingly, it should be emphasized that carrying out initial genetic monitoring of the ACQUERINO deer represented a significant opportunity to predict any negative effects on the population viability (Borowski et al. 2016). Future phylogeographic studies involving the application of nuclear markers, as microsatellites, should be conducted to test whether this population is subjected to inbreeding depression. The results obtained in this paper integrated with eco-ethological data will be essential for the management of the species and to avoid choices that may have repercussions on the well-being of the population. In other words, the best way to avoid any future inbreeding depression is to promote the range expansion of this population (e.g. Kenney et al. 2014). Based on these results, it would be useful to examine a greater number of samples and other molecular markers, to understand what are the gene flows between the various Apennine nuclei and the actual genetic structure of the species in this region.
Acknowledgements Professor Frank Zachos and Stefano Mattioli kindly provided us with useful comments on an early draft. Bruno Esattore and Andrea Viviano helped us with subspecies identification. Dr. E. Basset (University of Birmingham) revised our manuscript and improved the English grammar and syntax. Two anonymous reviewers and the Editor Magdalena Niedzialkowska kindly improved our MS with their useful comments.
Author contribution FD conceived the idea; FG collected data; MB conducted the genetic analyses and organized the dataset; all authors wrote the MS. All authors approved the final draft before submission.
Funding Tuscany Region, Hunting and Fishing Department.
Data availability All data are available by the corresponding author. Tissue samples are stored at the Institute of Biosciences and Bioresources -CNR, Sesto Fiorentino (FI) I-50019.
Code availability Not applicable.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.