Conservation genetics of American crocodile, Crocodylus acutus, populations in Pacific Costa Rica

Maintaining genetic diversity is crucial for the survival and management of threatened and endangered species. In this study, we analyzed genetic diversity and population genetic structure at neutral loci in American crocodiles, Crocodylus acutus, from several areas (Parque Nacional Marino Las Baulas, Parque Nacional Santa Rosa, Parque Nacional Palo Verde, Rio Tarcoles, and Osa Conservation Area) in Pacific Costa Rica. We genotyped 184 individuals at nine microsatellite loci to describe the genetic diversity and conservation genetics between and among populations. No population was at Hardy-Weinberg Equilibrium (HWE) over all loci tested and a small to moderate amount of inbreeding was present. Populations along the Pacific coast had an average heterozygosity of 0.572 across all loci. All populations were significantly differentiated from each other with both FST and RST measures of population differentiation with a greater degree of molecular variance (81%) found within populations. Our results suggest C. acutus populations in Pacific Costa Rica were not panmictic with moderate levels of genetic diversity. An effective management plan that maintains the connectivity between clusters is critical to the success of C. acutus in Pacific Costa Rica.


introduction
Threatened and endangered species face many challenges including habitat fragmentation and destruction, human population growth, and loss of genetic variability. Maintenance of genetic diversity is of increasing importance in the preservation of threatened and endangered species (Lacy 1997;Haig 1998;Reed and Frankham 2003;Reed et al. 2007). Lack of genetic diversity can lead to inbreeding depression (Frankham 1995), decreased immunity (O'Brien et al. 1985), decreased reproductive performance (O'Brien et al. 1985;Parker et al. 1991) and eventual extinction (Frankham 2005). Effective management strategies for threatened and endangered species require integration of all aspects of the species' biology, including both demography and genetics (Lande 1988) The American crocodile (Crocodylus acutus) is mainly a coastal species ranging from the extreme southern tip of Florida, through the Caribbean, and Central and northern South America (Mazzotti 1999;Thorbjarnarson 2010). Populations rangewide are threatened by habitat destruction and fragmentation, poaching, and past overexploitation (Ross 1998;Thorbjarnarson et al. 2006;Mazzotti et al. 2007). Current threats to C. acutus in Costa Rica include habitat loss (particularly nesting habitat), deliberate killing (Machkour-M' Rabet et al. 2009), and pollution (Rainwater et al. 2007;Rainwater et al. 2011). Crocodylus acutus is a wide ranging and ecologically plastic species (Thorbjarnarson 2010) with substantial genetic differentiation among populations (Menzies and Kushlan 1991;Rodriguez 2007;Porras et al. 2008;Thorbjarnarson 2010). Determining the status and ecology of C. acutus in Costa Rica is a priority project of the IUCN Crocodile Specialist Group Action Plan (Ross 1998). Genetic evaluations of C. acutus range-wide were further named as a moderate priority project in the 2010 Action Plan (Thorbjarnarson 2010). Therefore, it is important to describe genetic diversity and differentiation of this species throughout its range, including Costa Rica.
In the present study, we investigated the genetic structure of C. acutus populations in several areas of Pacific Costa Rica (Fig. 1) and the degree of effective migration occurring between populations. A series of estuaries provide pockets of suitable crocodile habitat along the Pacific coast which made it an optimal area for studying gene flow between potential metapopulations. Crocodiles are known to migrate long distances (Kay 2004;Campos et al. 2006;Read et al. 2007;Campbell et al. 2013). There have been accounts of C. acutus migrating over 35 km for nesting (Cherkiss et al. 2006) and movements over 388 km (Cherkiss et al. 2014) in southern Florida. The ability of crocodilians to migrate and disperse long distances increases the amount of potential gene flow between neighboring or distant populations. It is possible that crocodiles are dispersing between habitat patches in Pacific Costa Rica; thus, facilitating gene flow along the coast. We used microsatellites to test the hypothesis that C. acutus populations do not exist as a continuous, panmictic population in Pacific Costa Rica.

Study area
We sampled 5 localities in Pacific Costa Rica for C. acutus (Fig. 1). Site LB (Parque Nacional Marino Las Baulas) was in the Tempisque Conservation Area (ACT); site PV (Parque Nacional Palo Verde) was in the Arenal-Tempisque Conservation Area (ACA-T); site SR (Parque Nacional Santa Rosa) was in the Guanacaste Conservation Area (ACG); site RT (Rio Tarcoles) was in Central Pacific Conservation Area (ACOPAC); sites RS (Rio Sierpe), T (Terraba Delta), PL (Pejeperro Lagoon), PTL (Pejeperrito Lagoon), RE (Rio Esquinas), RC (Rio Coto) and PB (Parrot Bay Lodge, Puerto Jimenez) were in the Osa Conservation Area (ACOSA). Crocodiles were sampled from seven areas in ACOSA; however, they have been combined as one population due to low sample numbers. Localities ranged from large river systems (PV, RT, and ACOSA), to estuaries (LB and SR) and coastal lagoons (SR and ACOSA). (See Mauger et al. 2012 for additional study location information.)

Sample collection
We collected blood and tissue samples at the beginning of the rainy season in SR (2007) and PV (2005, 2008 and 2009), throughout the year in LB (2007 -2009), during the rainy season in RT (2005RT ( -2006 and during the end of the dry season in ACOSA (2006, 2008 and 2009). We captured crocodiles mainly during spotlight surveys using the break-away snare method (Hutton et al. 1987;Hutton and Woodhouse 1989), snake tongs or by hand (see Mauger et al. 2012 for additional information on sample collection). Blood and/or tissue was collected from 184 individuals (see Table 1 for size class distribution of samples). In sites where a large number of hatchlings were captured, a random number selector was used to randomly select four to six hatchlings for genetic analysis. Tissue was collected from the caudal scutes during marking and blood was collected from the caudal vein or the dorsal sinus. Tissue samples were preserved in 95-100% ethanol. Blood was preserved on Whatman FTA Cards for DNA Preservation® (GE Life Sciences).

DNA isolation and microsatellite amplification
We isolated DNA from tissue samples using the DNeasy Blood and Tissue Kit™ (Qiagen) and purified from blood cards with two five-minute washes with FTA Purification Reagent (Whatman) and two five-minute washes with Tris-EDTA (TE; 10 mM Tris-Cl, pH 7.5, 1 mM EDTA) buffer. Each wash consisted of 50 µl of solution.

Genetic diversity
Data files were converted to formats supported by various genetic programs in CREATE 1.0 (Coombs et al. 2007). Probability of Identity (PI) was estimated in GENALEX 6 (Peakall and Smouse 2006) to determine the minimum number of microsatellites needed to identify individuals. Allelic richness (A R ) and the number of private alleles (A Priv ) were estimated in HP-RARE (Kalinowski 2005). We estimated numbers of alleles, allele frequencies and gene diversities in FSTAT 2.9.3.2 (Goudet 1995).
Observed versus expected number of heterozygotes were estimated in Genepop on the Web (Raymond and Rousset 1995;Rousset 2008). Departure from Hardy-Weinberg Equilibrium (HWE) and linkage disequilibrium (LD) were estimated in Genepop on the Web (Raymond and Rousset 1995;Rousset 2008). Departure from HWE was tested using an exact test (Guo and Thompson 1992) and a chi-square goodness of fit test with a dememorization number of 10,000, and 1,000 batches of 10,000 iterations each. Linkage disequilibrium was tested to determine if small effective population sizes within the different localities caused nonrandom association of alleles at different loci. Linkage disequilibrium was tested for all pairs of loci used by the log likelihood ratio statistic under the same parameters as HWE. All p-values were adjusted to allow for multiple comparisons. Weir and Cockerham's (1984) inbreeding coefficient, F IS , was estimated for each population in FSTAT 2.9.3.2 (Goudet 1995) with and without randomly selected hatchlings. Homozygote excess at each locus in each locality was estimated by MICROCHECKER (Van Oosterhout et al. 2004). MICROCHECKER was also used to identify if null alleles were present at each locus in each locality. Null alleles were suggested for loci with a general excess of homozygotes for most allele size classes.
We identified whether rare alleles had been lost due to previous genetic bottlenecks under the infinite alleles (IAM) and stepwise mutation (SMM) models in BOTTLE-NECK version 1.2.02 (Piry et al 1999). Both the standardized differences and Wilcoxon tests were run.

Population genetic structure
We used an analysis of molecular variance (AMOVA) to estimate the percentage of variance within and among populations with GENALEX 6 (Peakall and Smouse 2006). Population differentiation was estimated for all population pairs using several methods. We estimated F ST and R ST for each population pair in FSTAT 2.9.3.2 (Goudet 1995) and Arlequin ver. 3.11 (Excoffier et al. 2005) respectively. All p-values were adjusted to allow for multiple comparisons.
We used Mantel's test to determine the relationship between geographic and genetic distance. Isolation by Distance Web Service (IBDWS; Jensen et al. 2005) was used to test for the presence of isolation by distance (IBD) between population pairs. Sites RT and T were excluded from the analysis because GPS coordinates of crocodile captures were not available. Each Mantel test was performed with 30,000 randomizations. Distances between populations were estimated using an oceanic/coastline route. Rousset's genetic distance (F/(1-F)) was calculated using genetic differentiation (F ST ).

Genetic diversity
The nine microsatellites chosen for this study had an average probability of identity (PI) of 4.96 -6 across all five Costa Rican populations (SR=5.4 -6 ; LB=7.9 -6 ; ACOSA=8.7 -8 ; PV=1.5 -7 ; RT=1.1 -5 ). This indicated that there was a low probability that two individuals chosen at random would have the same genotype. These microsatellites were sufficient for this study.
We identified 88 alleles in five C. acutus populations sampled in Pacific Costa Rica across all nine microsatellite loci. Average A R and A Priv over all loci were estimated using a corrected sample size of 34 alleles. The A R ranged between 4.22 and 5.64 and A Priv ranged between 0.27 and 1.36 in the sampled locations (Table 2). Allele frequencies for each microsatellite locus ranged from 0 (in localities where the allele was not genotyped) to 0.882 (Appendix 1: Allele Frequencies). Allele frequencies were also calculated with the hatchlings removed. There were no substantial differences in allele frequencies when hatchlings were removed. We tested for genetic bottlenecks under the IAM and SMM models for all samples combined and each site separately. Bottlenecks were not detected under IAM (p>0.05), but were detected in all populations under SMM (p<0.009) with the standardized differences test.

Population genetic structure
An analysis of molecular variance estimated that 19% of the variation occurred between populations, while 81% of molecular variance occurred within individual populations. This suggested that individual populations were genetically diverse. Population differentiation was measured between all population pairs using F ST and R ST . All population pairs were significantly differentiated (p=0.05) using both measures of population differentiation (Table 4). We observed the least amount of differentiation between ACOSA with LB and SR and the highest level of differentiation between LB with SR and RT.
Isolation by distance (IBD) was estimated between each sampled crocodile in all localities in Pacific Costa Rica. No IBD was observed (p=0.92; Fig. 2). Isolation between populations did not restrict gene flow.

Discussion
The nine microsatellites chosen in this study provided data on the genetic structure of C. acutus populations along the Pacific coast of Costa Rica. Average heterozygosity of crocodiles along the Pacific coast of Costa Rica was slightly higher than or comparable to that in other crocodilian populations (Glenn et al. 1998;Davis et al. 2001;Dever et al. 2002;Ryberg et al. 2002;Verdade et al 2002;de Thoisy et al. 2006;Rodriguez et al. 2008). However, in this study, several individual loci did have lower heterozygosity values, possibly due to inbreeding observed in all populations. Crocodiles within site SR had higher F IS values than other sites. This may be because we observed and captured few individuals that had exceeded minimum breeding size. This site also represented the smallest crocodile habitat with lower encounter rates (Mauger et al. 2012). Fewer breeding crocodiles and lower encounter rates could explain the higher inbreeding levels since presumably the gene pool is limited to fewer individuals. In all surveyed localities, size class distributions estimated during the study showed a higher percentage of juveniles (0.5 -1.25 m) than adults (>2.25 m) (Mauger et al. 2012). The higher percentage of juveniles in these localities could indicate population recovering from past bottlenecks (Ouboter and Nanhoe 1989) and explain the higher inbreeding levels observed in this study.
We identified previous genetic bottlenecks in all sampled localities. Past bottlenecks were confirmed under the SMM model by BOTTLENECK (Piry et al. 1999). Based on these analyses, we concluded that all populations underwent a previous reduction in population size. Population bottlenecks lead to the rapid loss of rare alleles and results in the loss of the total number of alleles at a faster rate than a loss in overall heterozygosity (Ortego et al. 2010). The low number of private alleles observed in this study (Table 2) could be an artifact of past genetic bottlenecks. Null alleles were also detected for at least one locus in each locality (Table 3). The absence of these alleles in the analysis could help explain the loss of genetic variation observed in these localities. Crocodylus acutus populations declined range-wide through the mid-20 th century as a result of hunting and illegal poaching (Ross 1998;Thorbjarnarson et al. 2006;Thorbjarnarson 2010). These human activities have not been documented in Costa Rica; however, it is possible that populations here experienced similar pressures and population declines. We observed some poaching and killing of large individuals in LB during the study and have received reports of similar activities at other sites in Pacific Costa Rica (personal observation and communication with local people). Anecdotal data suggests that crocodile numbers are also increasing in LB (F. Paladino, personal communication) and PV (Bolaños-Montero 2012). Recent introduction of tilapia into the Tempisque River Basin (site PV) has provided a continuous food source for crocodiles (Sandlund et al. 2010), potentially contributing to the recent population growth. As a result, crocodile numbers in this region have increased precipitously in recent years, causing adverse interactions between human and crocodile populations (Bolaños-Montero 2012). Additional factors are most likely at play, which are contributing to population growth.
No population was in Hardy-Weinberg Equilibrium for all microsatellite loci. This could be due to the heterozygote deficiency observed at loci that were not in HWE, to inbreeding documented in several populations, or to effective migration. Our results suggest that crocodile populations in Pacific Costa Rica were differentiated from each other, as supported by studies on C. acutus and other crocodilian species (Farias et al. 2004;Porras et al. 2008;Machkour-M'Rabet et al. 2009;Thorbjarnarson 2010). All population pairs exhibited significant genetic differentiation, with all but two population pairs (LB and ACOSA, and SR and ACOSA) showing moderate differentiation from each other (R ST <0.05; Table 4). These values supported the hypothesis that migration occurred between populations; however, some population pairs were more differentiated from each other and thus had lower historical migration rates between patchily distributed habitats in Costa Rica. The level of subdivision observed suggested that crocodiles in Pacific Costa Rica were not panmictic; however, genetic connections did exist. Additionally, the majority of molecular variance was observed within populations. This could explain why moderate differentiation was observed between most population pairs. However, the highly mobile nature of crocodiles (Kay 2004;Campos et al. 2006;Read et al. 2007;Campbell et al. 2013) could be facilitating gene flow between populations along the Pacific coast. A recent study on the spatial ecology of C. acutus in Panama, suggests that males have a larger home ranges, but females have larger average movement distances (Balaguera-Reina et al. 2016). Balaguera-Reina et al. (2016) also noted dispersal differences between age classes and dry and wet seasons. Dispersal differences between age classes, i.e. subadults dispersing to find mating territories, could also explain the departure from Hardy-Weinberg Equilibrium and moderate differentiation levels. GPS-based tracking studies of Costa Rican C. acutus would contribute important information on contemporary crocodile dispersal abilities and maximum home ranges in patchily distributed habitats.

Conclusions
The data presented here supported moderate differentiation and an absence of isolation by distance in Pacific Costa Rica. Our results suggested the loss of genetic variation through a lack of connectivity between some localities and previous population bottlenecks. The moderate heterozygosity values and genetic differentiation described here emphasized the need to protect all potential crocodile habitat, to write management plans across conservation areas and national parks in Costa Rica, and the need for conservation and management units to extend over the entire span of a species' range.

Abstract
Overuse of natural resources by humans is a major threat to biodiversity. Overuse often involves species of economic or esthetic value, and fish are a typical example for a group that is exploited both for economic reasons (for human consumption) and for esthetic reasons (e.g. by aquarists). Pseudorhodeus tanago (Tanaka, 1909) (formerly known as Tanakia tanago) is a small colorful but legally protected (fishing, keeping and transfer are banned) bitterling fish distributed around Tokyo, Japan. Whereas it is critically endangered and more and more habitat loss has occurred, at least four stocks have been newly found during the last decade. To explore whether emergence of these newly found habitats is a consequence of incomplete survey, we genotyped mitochondrial cytochrome b sequence of P. tanago from 17 localities and an illegal home aquarium. Populations known by the past extensive survey (13 localities) showed geographically structured population genetic characteristics. Population-specific haplotypes were common indicating past divergence and bottleneck events. Four (north, {center + west}, south_1, south_2) or five (north, center, west, south_1, south_2) geographic groups were detectable as for these known localities.
On the other hand, newly found stocks were polymorphic and showed identical haplotypes from distant known localities. If we assume historical basis of distribution and genetic characteristics of these newly found stocks, it must be a series of unlikely geological events and haplotype sorting. We discuss potential issues posed by these questionable stocks.

Keywords
Bottleneck, bucket biology, conservation genetics, fish dumping, home aquarium, phylogeography, poaching, satoyama, Tanakia tanago introduction Overuse of organisms by hunting or fishing for trade or esthetic purposes is one of the biggest threats to biodiversity. Controlling these activities is fundamental for ensuring the persistence of endangered organisms, particularly those that have traits attractive to humans. Alerts to poaching are thus necessary on those organisms legally protected. Not only tusks of elephants and rhinos but beautifully colored bodies of fish as well have attracted violators. The red list of Japanese brackish and freshwater fish (Ministry of the Environment, Japan 2015) acknowledges 39 species out of 168 of endangered and vulnerable categories have threats of overfishing in connection with home aquarium. Inversely, unauthorized fish release also threatens the integrity of natural fish population structures. A big source of invasion of alien freshwater fish is fish dumping from home aquaria (Lintermans 2004, Gertzen et al. 2008, Fuller et al. 2013, Ishikawa and Tachihara 2014. Poaching and dumping of endangered and protected fish, if any, disturb conservation programs and policy making and ultimately threat that species. Pseudorhodeus tanago (Tanaka, 1909) (formerly known as Tanakia tanago transferred to the new genus, Chang et al. 2014) is a colorful small bitterling fish endemic to Japan with a limited geographic range around Tokyo. Habitats of P. tanago are in small water bodies such as ditches or ponds with spring water in hills (Nakamura 1969, Mochizuki 1997, Maita 2002, Ishinabe 2014. These habitats link with traditional agricultural landscape known as 'satoyama'. Traditional farming activities have maintained ditches, ponds and vegetation around waters in farmlands in hills (Kobori and Primack 2003), and thus maintained habitats of P. tanago. However, because of heavily populated and highly developed areas around its habitats, P. tanago is critically endangered. Urbanization accompanying change of farming style gives rise to habitat loss (Mochizuki 1997, Ishinabe 2014). Agency for Cultural Affairs, and Ministry of the Environment, Japan have endorsed this bitterling as a legally protected species (Ministry of the Environment, Japan 2015). Fishing, keeping and transfer of P. tanago are banned. Since around 1970, conservation measures of P. tanago have been taken under local bases including extensive search for habitats, development of ex situ breeding techniques, protected area enclosure, mitigation of civil engineering activities (Tochigi Prefectural Fisheries Experimental Station 1973, Akiyama et al. 1994, Mochizuki 1997, Kubota et al. 2010, Ishinabe 2014. A national breeding program for the protection of the species began in 1995, which included ex situ breeding and habitat restoration (Agency of Environment, Japan et al. 1995).
Records of 45 localities of P. tanago habitats were established as part of an extensive survey (Tanaka 1909, Nakamura 1969, Maita 2002, Ishinabe 2014 report) (Fig. 1a). Some of the localities have been kept secret by authorities because of risks of poaching. By around 2000, most known habitats were lost, and further local extinction has occurred in these days in spite of conservation measures (Fig. 1b). Many populations have undergone habitat degradation to decline to ca. 1/100 in a few years, and in an extreme case the population is vanishing in the final remnant habitat of a short (< 100 m) stretch of a small (ca. 60 cm wide) stream (Mochizuki 1997).
From the habitat characteristics of P. tanago, confined and scattered among headwaters of fine branching dales in agricultural landscape on hill terrains, we assume geographically structured population genetic characteristics of this species. Kubota et al. (2010) outlined geographic population structure of this species, but lack of localities from southern part of its geographic range and usage of questionable specimens obscured clear-cut geographic structure. Analysis of geographic population structure of this critically endangered fish with habitat loss based on specimens from reliable sources is necessary for setting up conservation programs (Waples and Gaggiotti 2006).
Whereas more and more local populations have been lost recently, P. tanago were newly found in a few localities in the last decade ('newly found' locality or stock hereafter) (Fig. 1b, hatched areas). Are these newly found stocks are simply because of incomplete survey of distribution? The aims of this study were to delineate geographical genetic structure of P. tanago and to identify symptom of disturbance, if any, in this structure by an analysis of genotypes of 80 individuals from 18 populations or stocks. Four of them are from newly found localities, one from an illegal home aquarium, and the others are from localities previously known by the past extensive survey ('known' localities or populations hereafter).

Materials and methods
Pseudorhodeus tanago from 13 known localities collected from 1993 through 2013 (55 individuals), four newly found stocks collected from 2010 through 2014 (22 individuals), and three individuals seized from an illegal home aquarium were materials of this research (Table 1, Fig. 1b). Fishes from two out of the 13 known localities were of ex situ preserved stock extinct in the wild in early 1990s (#8) and late 1970s (#10). Habitat characteristics of the collecting localities were identical to descriptions in literatures (Nakamura 1969, Mochizuki 1997, Maita 2002) except for two newly found localities (#14, 17) where civil engineering activities including straight cut of channels with concrete enforcement have suffered their habitats.
Extraction of DNA was done from fin clips with QuickGene DNA Tissue kit on QuickGene-810 (Kurabo, Neyagawa, Japan). PCR primers were L14695 on the L-strand (AATTYTTGCTCRGACTCTAACC) and H15910 on the H-strand   Table 1. Asterisks indicate extinction in the wild and specimens came from ex situ preserved stocks. Habitats #14-17 in hatched areas (exact places are kept secret by authorities) are newly found. Hatched areas thus do not represent exact geographic ranges.
(GATCTTCGGATTACAAGACCGAT) which worked for amplifying a 1219 bp mitochondrial DNA fragment encompassing the whole cytochrome b gene and flanking tRNA partial sequences. PCR reaction mixture of 12.5 µL contained 1µL of template DNA, 0.96 µL of dNTP mix (2.5 nmol each), 1.2 µL of 10x ExTaq buffer, 0.06 µL (0.3 U) of ExTaq (Takara, Shiga, Japan), 1 µL of primers (5 pmol each), 7.28 µL of Milli-Q grade water. PCR reaction was of touchdown profile (Don et al. 1991) in which annealing temperatures dropped from 59°C down to 53°C in the initial 7 cycles and was constant at 55°C in the remaining 28 cycles (35 cycles in total). The PCR reaction started with 3 min at 94°C followed by 35 cycles of 30 sec at 94°C, 30 sec at the above touchdown annealing temperatures, 120 sec at 72°C with final extension at 72°C for 5 min.  et al. 2002). The net nucleotide divergence was calculated by π xy -(π x +π y )/2 where π xy is average number of nucleotide differences between populations x and y, and π x and π y stand for this value between individuals within populations x and y. Parsimonious haplotype network was drawn with TCS v.1.2.1 (Clement et al. 2000).

Haplotype grouping
Sequencing Pseudorhodeus tanago mitochondrial cytochrome b revealed 10 haplotypes (Table 2, Hap01-10). There were 19 variable sites, and 14 nucleotide substitutions observed between the most distant haplotypes (Hap01, 10) ( Table 3). All the observed nucleotide substitutions were transitions, two of which were non-synonymous (positions 15100 and 15415). Most specimens from known localities (11/13) were monotypic, whereas newly found and seized stocks except for #17 with a single individual examined were polymorphic.
A haplotype network coincides with geographically structured population genetic characteristics of this species (Fig. 2). Regarding 13 known localities, haplotypes from each locality frequently were specific to that locality (9/13). Haplotypes placed on the upper-left side of the figure (Hap09, 10) were from northern part of the geographic range (Fig. 1b, #12, 13). Likewise, haplotypes arranged on the right side (Hap01-05) were all from the southern part (Fig. 1b, #1-6), and Hap08 on the lower-left side of the network was from the western part (Fig. 1b, #10) of the range. Placement of the other two haplotypes (Hap06, 07) was on the center of the network, and their geo- graphic placement was roughly the center of the range (Fig. 1b, #7-9) except for a single exception from locality #6 in the southern part (Table 2). This haplotype grouping was delineated by gaps with three or more missing haplotypes, whereas one or no gap lay between haplotypes within these groups. Haplotypes of the south group, however, contained distant types (four nucleotide differences at most), and there might be subgrouping among them.
(ca. 1-16 km, Table 5) divided the southern localities into two geographic groups, whereas moderate differences (0 to four base changes) among localities of the center and west geographic groups distant from each other (ca. 60-100 km except between #7 and #8) could be pooled. Geographic grouping of known localities was largely coincident with river connectivity (Fig. 3). The North and the West geographic groups appeared each in a single river basin respectively. The South geographic group covered two river basins close to each other, one of which is kept secret by the authority (not shown on the map). Genetic (π xy , below diagonal) and geographic (km, above diagonal) distance among collecting localities.  Potential connection between these two basins through deltas is possible. The Center geographic group spanned a wider range on the periphery of a lowland plain which roughly corresponds with the area encircled by the gray solid line of this group plus coastal area on the east (Fig. 3). Localities of this geographic group presented in three river basins potentially connected with each other through flat low overflow terrain of 12 m high (Fig. 3, arrow a) or overflow plus switching dales (Kagose 1979) (flat hill top at 90m high and opposing headwaters with imbalanced slopes at 60 m high, arrow b).
On the other hand, placement of specimens from newly found localities onto the haplotype network obscured the geographic population genetic structure, though SAMOVA analysis apportioned these localities to either of the South or Center geographic groups according to their haplotype composition. Coverage of variation among geographic groups reduced from 74.43-80.5% to 72.76-76.18% (Table 4). Newly found localities except for #14 contained haplotypes of only the South group in spite of potential river connectivity with localities with the Center or the North geographic groups (Fig. 3). Locality #14 showed a mixed composition of both the Center and the South haplotype groups.

Genetic architecture of known populations
The low sequence diversity of Pseudorhodeus tanago represented as monotypy in most known localities (11/13) conformed to their habitat characteristics (Table 2). Specific haplotypes from many of the known localities (9/13) indicated that populations of P. tanago rarely exchange with each other. Confinement in small isolated water bodies and subsequent bottleneck events might have brought about the present genetic structure of populations.
Isolation on a wider and geological scale is also responsible for divergence among geographic groups. The North geographic group contained haplotypes distant from others (Hap01 and Hap10), and ancestors of the North and other geographic groups thus diverged first. The central part of the entire geographic range of P. tanago encircled the flat plain that underwent repeated marine transgressions (Fig. 3). Marine transgression at > 400 ka (Sugai et al. 2013) and/or change of river flow at 500 ka (Fig. 3, blue arrow) (Koike et al. 1985, Kubota et al. 2010) might be responsible for the isolation of these geographic groups. Repeated marine transgressions afterward have facilitated further isolation among geographic groups.

Implication of genetic characteristics of newly found stocks
Polymorphism and haplotype composition among newly found stocks, on the other hand, posed a question about the characteristics of genetic and geographic population structure of P. tanago. If we assume historical basis of distribution and genetic characteristics of these newly found stocks, it must be a series of unlikely geological events and haplotype sorting.
Including these newly found stocks, P. tanago as a whole showed a bipartite population structure in which some localities show genetic variation while others rarely do so (Table 2). We should assume population admixture and rampant migration among these newly found stocks or one way migration from northwest part of the map including known populations (Fig. 1a), whereas known populations underwent isolation and bottleneck events. Marine transductions were, however, more frequent near the east coast since geological ancient until historic medieval times where newly found stocks inhabit (Kubo 2007, Sugai et al. 2013, and thus habitats of newly found stocks might be ephemeral prone to bottleneck. River connectivity implies that the hypothetical migration was from central or north geographic group (Fig. 3), but newly found localities (#15-17) consisted solely of south haplotype group. Admixture between populations around south group of known and newly found localities upon ace age marine regression is unlikely because of absence of continental shelf at the south-east coast (Sugai et al. 2013). Even if it was possible, river connection should involve localities of center geographic group (#7, 8) (Fig. 3), the main haplotype of which is absent from the newly found localities. Among newly found stocks, #14 consisted of haplotype of center group and was close to populations of center geographic group (#7, 8), but its natural distribution is questionable because of disturbance of the habitat by civil engineering. Straight cut channelization with concrete enforcement both on the banks and bottom has eliminated mud-sandy substrates for bitterlings' spawning host mussels.
There is no information such as confession or witness of violators about unauthorized releasing of P. tanago into the newly found localities. Genetic polymorphism of newly found localities may arrange stocks of these localities as conservation Natural distribution in these newly found localities is, however, questionable because of unlikely association between haplotype/population relationships and river connectivity. Because bitterling fishes have attractive coloration to humans, there have been reports on unauthorized intentional releasing activities possibly implicated in home aquaria (Miyake et al. 2011, Kitazima et al. 2015, Saitoh et al. 2016. Under an assumption of unauthorized releasing or dumping of P. tanago, we confront difficult issues. Genetic polymorphism of the newly found and seized stocks (Table 2) implied a number of poaching activities. Because violators do not assess stock size of their target fish in the field, poaching would sometimes be devastating. Unauthorized releasing or dumping, if any, also poses problems. Unauthorized releasing of P. tanago is illegal, because unauthorized transfer of this legally protected fish is prohibited. Released P. tanago individuals themselves are, however, legal target of conservation (Ministry of the Environment, Japan 2015). Presence of questionable stocks thus disturbs prioritization of conservation targets under limited funds and human resources. Releasing or dumping onto habitats of natural populations is a direct threat to the integrity of natural properties of P. tanago. Our report has cautionary implication to potential violators.

Acknowledgments
We dedicate this paper to Toshihiro Ishinabe (deceased) (Kannonzaki Nature Museum) and his contributions to conservation of Pseudorhodeus tanago over two decades carrying out habitat monitoring and ex situ breeding. KS, NS and KI conceived this research. KS analyzed data and drafted the manuscript. NS, MO, KI, TS, TM, TT and MT conducted field survey. All authors read, improved and approved the manuscript. This research was financially supported by a Protective Breeding Program of Ministry of the Environment, Japan. Permission for research use of P. tanago was provided by Agency for Cultural Affairs and Ministry of the Environment, Japan to NS and KI. We have no financial conflict of interests about contents of this publication. Specimens were collected under a condition in which the authorities expect non-disclosure of exact places of habitats except for a protected and publicly opened locality (#13). NS and KI are employees of two of these authorities that do not disclose some localities.

temporal changes in vegetation of a virgin beech woodland remnant: stand-scale stability with intensive fine-scale dynamics governed by stand dynamic events
The very intensive fine-scale dynamics of the studied beech forest are profoundly determined by natural stand dynamics. Extinction and colonisation episodes even out at the stand-scale, implying an overall compositional stability of the herbaceous vegetation at the given spatial and temporal scale. We argue that fine-scale gap dynamics, driven by natural processes or applied as a management method, can warrant the survival of many closed forest specialist species in the long-run.
Nomenclature: Flora Europaea (Tutin et al. 2010) for vascular plants; Soó 1968Soó -1980  The significance of the herbaceous layer in forest biodiversity and ecosystem functioning has been widely appreciated (Whigham 2004, Gilliam 2014, generating an interest in its dynamics. During the last decade, remarkable changes in the herbaceous layer have been shown by several studies (e.g. reviewed in Verheyen et al. 2012). Resurvey studies investigated the understorey to detect potential long-term changes and to assess the role of different environmental variables as driving forces in deciduous forests. The observed changes were attributed to different background causes, such as combined effect of temperature increase and canopy opening (De Frenne et al. 2013, abandonment of former forest management practices (Hédl et al. 2010, Heinrichs et al. 2014), lack of disturbance (Brewer 1980, Taverna et al. 2005, deer herbivory ( (Brewer 1980) and herbivory (Rooney and Dress 1997). A loss of rarer native plants (Hédl 2004) and the homogenization of the herbaceous layer was detected, and was attributed to high browsing pressure by deer herbivory (Wiegmann and Waller 2006), acidification (Durak 2010) and light deficit (Davison andForman 1982, Heinrichs et al. 2014). Fine-scale herbaceous layer dynamics is also often linked to stand-scale dynamics. Gap formation significantly increases solar radiation and soil moisture (Collins et al. 1985, Moore and Vankat 1986 leading to intensive alteration of composition and cover of the herbaceous layer , Kelemen et al. 2012) in large (Degen et al. 2005) and small gaps (Mountford et al. 2006). Stand dynamic events were matched with certain vegetation changes in different Slovakian primeval forests (Ujházy et al. 2005, Martináková and Martinák 2012. The herbaceous layer was best developed at the beginning of the decaying stage, and its cover sharply decreased when the tree seedlings and saplings overgrew them at the end of the decaying stage and in the growing stage. Herbaceous layer vegetation change can have a cyclic pattern, especially in forest stands that do not experience the strong successional effect of abandonment or other marked environmental change. Ujházy et al. (2007) detected low species turnover and cyclic changes in a Slovakian primeval forest at the stand-scale. Changing abundance or dominance of herb species characterized the cyclic changes, which was interpreted as internal community dynamics (Van Der Maarel 1996, Ujházy et al. 2005. Carøe et al. (2000) could not detect a successional trend either, only changes in species abundance in a managed Danish beech forest. Neither could Martináková and Martinák (2012) in a Slovakian natural fir-beech forest; however, the time interval was short in both investigations.
These short-and long-term changes in the herbaceous layer are most often detected by resurveys of phytosociological relevés or smaller quadrats after a few years (Ujházy et al. 2007, Martináková andMartinák 2012), one or a few decades (Davison and Forman 1982, Taverna et al. 2005, Łysik 2008, Vanhellemont et al. 2014 or even after a half-century interval (Brewer 1980, Hédl 2004, Durak 2010, Hédl et al. 2010, Šebesta et al. 2011. In some cases the relocation can be precise, when the time span is short or permanent corner markers were used, but more often former plots are relocated only approximately. The error resulting from relocation uncertainty can be high, especially in the case of small plots . However, if relocation is done carefully, the resurvey can be robust to localization uncertainty, and can provide evidence of long-term changes in the herbaceous layer (Chytrý et al. 2014, Kopecký and. Ross et al. (2010) also demonstrated that if temporal change in the vegetation is greater than recent spatial heterogeneity, results of the resurveys can be interpreted with some confidence.
In our study, we were especially interested in the dynamics of the herbaceous layer at multiple spatial scales in relation to simple stand dynamic events, such as the opening and closure of smaller and larger gaps of a primeval beech forest remnant in Hungary. The exact location of the original sampling plots surveyed in 1996 was unknown, but the sampling plots could be relocated with 1 m accuracy after 17 years.
We hypothesized that in this virgin woodland remnant: 1) fine-scale changes in the herbaceous layer are evened out at the stand-scale, where no significant changes can be detected in the species pool (low species turnover) and total cover of the herbaceous layer at this temporal and spatial scale; 2) fine-scale changes in the herbaceous layer are governed by dynamic events in the tree canopy, such as opening and closing of the canopy gaps; 3) small gap stand dynamics warrant the survival of closed forest specialist herbs at the stand-scale.

Study area
The study was carried out in Kékes Forest Reserve (63 ha), which is one of last and best remnant of central European montane beech woods in Hungary. Kékes is the highest point in Hungary (1014 m) and is situated in the Mátra Mountains, northern Hungary. The climate is relatively continental. Mean annual precipitation is around 840 mm, of which 480 mm fall during the growing season. Mean annual temperature is 5.7°C, with cold winter (-4.7°C in January) and mild summer temperatures (15.5°C in July). The bedrock is andesite and andesitic tuff. The topography is extremely steep; scree slopes are characteristic. The shallow brown forest soils are mainly covered by montane beech forest (Aconito-Fagetum Soó (1930Soó ( ) 1960. Mixed maple-ash-lime forest (Phyllitidi-Aceretum Moor, 1952subcarpaticum (Dostál 1933) Soó 1957 occurs in the most humid and rocky patches of the reserve on ranker type soils (Kovács 1968, Soó 1980. According to historical records (Czájlik 2009) the area has never been used for timber production. As a result, most parts of the reserve show the characteristic fine-scale mosaic of forest developmental stages -sensu Korpeľ (1995) -of central European montane beech forests. For the purpose of the original survey we selected a roughly 1.5 ha plot in the reserve that contained a closed old beech stand, a small gap in the beech stand, and a larger collapse on a rockier site.

Data collection
In 1996, vegetation of the herbaceous layer was systematically sampled in a 120 m × 120 m patch (approximately 1.5 ha), which was divided into 16 30 m × 30 m squares with wooden sticks as field marks in the four corners. On the 120 m × 120 m patch a grid with 5 m intervals was laid out. Altogether 576 plots with 0.25 m 2 were set out on the grid (see Suppl. material 1). All grasses, sedges, herbs as well as trees and shrubs lower than 0.5 m in height were included in the sample (Gilliam 2007). Total vegetation cover and cover for each vascular plant species were estimated in percentage.
To make future assessment of stand dynamics possible, in 1997 a tree stand position map (each tree individual recorded by X and Y coordinates) was created for the 1.5 ha plot (see Suppl. material 2). The map was also used to assist the approximate relocation of the 0.25 m² plots, as they were not permanently marked in 1996. Partial resurvey of the herbaceous vegetation was carried out in 2013. The aim was to compare 17-year vegetation change within and between the patches with different histories of stand dynamics. To select the appropriate patches, we used the stand position map and a series of archive aerial photographs to detect the location and time of certain events in stand development. Using the same methods, we resampled 306 plots representing five different stand developmental situations (see Suppl. materials 1 and 2). However, because of the lack of permanent marking of the original sampling locations, reloca-tion was only possible with a 1 m accuracy. To overcome this problem we used a cluster of four 0.25 m² plots at each of the 306 points to represent the original locations ( Figure 1).
The five stand developmental situations (SDS) were identified as follows (Table 1 The 306 plots resurveyed in 2013 were selected to include both the centres and the peripheries of the five SDSs (see Suppl. material 1). The results are presented for three different spatial scales as follows (Table 2): Stand-scale: Data from all the 306 plots are lumped together to characterise the 1 ha patch.
Stand developmental situation (SDS) scale: Each SDS is represented by the joined data of 25 plots located in the centre of an individual SDS. Thus, an approximately 400 m² area is sampled (which is comparable to the size of standard phytosociological relevé) for each SDS.    Fine-scale: All the 306 plots are treated separately so changes at the 0.25 m² scale can be described. Table 2 summarises the scales and numbers of sampling units in this study.

Data analyses
For our analyses we used the following variables: Cover of the herbaceous layer: For each plot we used the sum of estimated cover values of individual species, so > 100% values can occur. Mean cover and maximal cover were calculated from the plot data for the whole stand and for the five SDSs. For the 2013 data the overall means were used: 1224 plots (306×4) for the stand-scale, 500 plots (4×125) for the SDS-scale. Changes in the abundance of herbaceous vegetation were calculated as change in total cover. To get a better understanding of the dynamics of individual species at the stand-and SDS-scales, both net and absolute cover change were calculated.
Species richness of the herbaceous layer: The number of species was calculated for all the three spatial scales and for the two sample years. The average species number per plot was also compared.
To quantify changes that occurred during the 17 years we used the mean values of the four subsamples of 2013 data (mean species richness, mean cover of the herbaceous layer of each plot). In this way we used an equal number of plots as in 1996 for the comparisons (306 plots for stand-and fine-scale, and 125 plots for SDS-scale analyses). The relationships between stand developmental situations (SDS) and magnitude of change in species richness and total herbaceous cover were studied by using Kruskal-Wallis test (H statistics), the non-parametric analogue of classical analysis of variance (Conover 1980). Wilcoxon matched pairs test was used to test the differences between the two surveys regarding temporal changes of species richness and total herbaceous cover at SDS-and fine-scale.
In order to quantify the changes in a species pool the following variables were calculated, where in the case of 2013 data the average values of the four subsamples were used for all calculations: Number of colonisation events: Individual colonisation events were detected at the fine-scale (0.25 m²). The number of colonisation events (appearance of a species) was expressed as sum (all new occurrences of all species) and mean (average number of new occurrences per plot) for each spatial scale studied.
Number of extinction events: Individual extinction events were detected at the finescale (0.25 m²). The number of extinction events (disappearance of a species) was expressed as sum (all disappearances of all species) and mean for each spatial scale studied.
Absolute species turnover in the herbaceous layer: Absolute turnover in species composition between successive sampling years was calculated as (E + C)/2, where E and C are the number of species extinctions and colonisations, respectively (Williamson 1978).
Relative species turnover in the herbaceous layer: Relative turnover was calculated as (E + C)/(S1 + S2) × 100%, where E and C are as above and S1 and S2 are the number of recorded species present in the two years (Diamond 1969).
Colonization and extinction events and the species turnover in the herbaceous layer were also assessed in a qualitative way. The behaviour of individual species was analysed at the stand-and SDS-scales.
When the effect of relocation was analysed, we compared the size of time effect (i.e. the difference between the 1996 value and the average value of the four 2013 subsamples) and the size of relocation effect (i.e. average of paired differences between the four subsamples) by Wilcoxon matched pair test. A one-sided alternative hypothesis was applied, i.e. higher time than relocation effect. Significant differences imply that the observed differences between two sampling times cannot be caused by the imprecise re-allocation of the subsamples only.

Cover of the herbaceous layer
At the stand-scale we observed a general decrease in the abundance of herbaceous vegetation between 1996 and 2013 ( Figure 2). The mean total cover/plot changed from 20.2% in 1996 to 7% in 2013.
The observed changes in mean cover were obtained as the balance of positive and negative changes in individual species. Most changes were attributed to only a few species. As Table 4 shows, the sum of absolute cover changes was the highest for Galium odoratum, followed by Urtica dioica and the two fern species. A large difference between absolute and net change within an SDS indicates intensive fine-scale dynam-ics (increase and decrease in cover occur simultaneously in individual plots belonging to an SDS). This behaviour is characteristic of forest specialist species such as Galium odoratum (3YG, 1YG), Cardamine bulbifera (C, OBG) and Mercurialis perennis (3YG, 1YG, C, OBG)

Species richness of the herbaceous layer
In 1996 and 2013 we recorded 42 and 48 species, respectively, though in the latter case a four times larger area was sampled because of the four subsamples used (306 plots×4). At the stand-scale our full sample (306 locations, two surveys) contained data from 54 species (45 herbs, 6 trees, 3 shrubs), of which 50 species occurred in the 125 plots representing the five SDSs.
In contrast to the relative stability in species richness at the stand-scale, a much more variable picture was obtained when the SDS-scale was studied. In 1996 species richness varied between 7 and 23 (Table 5) in the five SDSs (25 plots each), whereas the same number were in the range of 14-31 in 2013. As Table 5 and Figure 3 show, the largest change in SDS-scale species richness took place in the three-year-old gap (3YG), where it grew from 7 to 31. Species richness was the most stable in the control, where it changed from 15 to 14. When mean species richness/plot was considered, the largest change was found in OCO, where it decreased to almost one-tenth of the original 3.7 in 1996.

Number of extinction and colonization events and species turnover in the herbaceous layer
Quantitative Species extinctions and colonisation occurred at all spatial scales studied during the 17 years between the two samplings. At the stand-scale there were only six species that dis-appeared, but all of them were present at very low frequencies and abundances in 1996. Half of the recorded 12 newly occurring species were actually present in the 576 plots sampled in 1996, but were not included in the sample (306 plots) used for the resurvey. At the SDS-scale there were considerable differences among the individual SDS. As Table 6 shows, the number of extinctions was the highest in OCO and OBG, whereas -as was expected -colonisation was the most intensive in 3YG. Mean absolute turnover was the highest in OCO and OBG, mostly because of the high number of extinctions. If we also consider the size of the species pool, relative turnover was the highest in OCO, followed by 1YG, where the initial species richness was very low. Fine-scale (0.25 m²) colonisation was in the range of 0 and 4.75 in the 306 plots, with a mean value of 0.717. The same values for extinctions were 0, 10 and 1.4, respectively.

Qualitative
Quantitative results indicate relative stability at the stand-scale and rather intensive dynamics at the SDS-and fine-scales. However, to get a better understanding of the processes, it is worth looking at what species become extinct or colonise in different situations. Special interest is devoted to closed forest specialist species (Schmidt et al. 2011), which are marked in bold in the text below.
At the stand-scale, the six species that were not included in the 2013 sample were all forest species (Chrysosplenium alternifolium, Epilobium montanum, Hordelymus europaeus, Melica uniflora, Monotropa hypopitys, Viola odorata) that occurred at very low abundances in 1996. Half of the newly colonising 12 species were forest species (Actaea spicata, Campanula rapunculoides, Epipactis helleborine, Hieracium acuminatum, Poa nemoralis and Polygonatum verticillatum) that were recorded as new because of similar chance effects since, according to our field records, they were all present in 1996 within the studied 1.5 ha stand. Among the colonisers there were two ruderal species (Tussilago farfara, Cirsium arvense) and one newly appearing adventive species (Impatiens parviflora).
At the SDS-scale it is worth comparing the list of disappearing and newly appearing species. In the OBG, most extinctions were of closed forest species (e.g. Galium odoratum, table 6. Number of extinction and colonisation events (overall mean of the four subsamples at each plot), absolute and relative species turnover in the five stand developmental situations (SDS) based on resurvey data collected in 1996 and 2013. For explanations of the stand developmental situations, see methods.

Mercurialis perennis, Viola reichenbachiana, Cardamine bulbifera)
and there were only a few extinctions of gap or disturbance indicating species (Chelidonium majus, Geranium robertianum, Rubus idaeus, Urtica dioica). Colonisation events were almost exclusively due to closed forest specialist species. At the rocky OCO site the observed extinctions were mainly attributable to those species (e.g. Dryopteris filix-mas, Impatiens nolitangere, Solanum dulcamara, Urtica, dioica) that occurred with high frequencies in 1996. The low number of colonisation was due to similar species because site conditions are not favourable for most closed forest specialist herbs. On the contrary, in the control plot the species pool did not change much, despite the relatively intensive dynamics (77 extinctions, 20 colonisations). Extinctions and colonisations were both due to real closed forest specialist herbs (e.g. Cardamine bulbifera, Euphorbia amygdaloides, Galium odoratum, Mercurialis perennis). In 3YG and 1YG only a low number of extinctions were recorded, mostly due to the disappearance of forest species. The large number of colonisations in these gaps was attributable to both typical gap species and to those forest species that are able to react relatively fast to changing light conditions. Table 7 shows how successfully closed forest species (Schmidt et al. 2011) survived in the individual SDS. Only a few species with low original frequencies disappeared, indicating the overall successful survival of this species group.

Effect of relocation error
Our results on the effects of relocation error showed that time effect was significantly higher than relocation effect both in the case of species richness (306 pairs of plots, Wilcoxon matched pair test results, T = 8475, Z = 7.769055, p < 0.00001) and in the case of total cover of the herbaceous layer (306 pairs of plots, Wilcoxon matched pair test results, T = 9932, Z = 6.983762, p < 0.00001). Consequently, differences between data from 1996 and average data of the four subsamples of 2013 cannot be merely attributable to imprecise relocation of subsamples.

Discussion
After 17 years we found low species turnover and a general decrease in herbaceous cover at the stand-scale. Therefore, our first hypothesis is only partly supported. We found an overall stability of the species pool in the herbaceous layer at the stand-scale, whereas the abundance of vegetation (measured as total herb cover) changed considerably.
Several analyses pointed out that a decrease in herbaceous cover in a temperate forest is attributed to light deficit caused by a denser canopy. This can be a consequence of less intensive forest management or even abandonment, as well as the lack of natural disturbance including the activities of grazers (Brewer 1980, Davison and Forman 1982, Vera 2000, Verheyen et al. 2012, Kopecky et al. 2013. A drastic decrease in herbaceous species cover was found by Łysik (2008) in a primeval beech forest in Poland and by Ujházy et al. (2007) in a primeval fir-beech forest in Slovakia. Both studies connected this sharp decrease with the massive recruitment of woody regeneration. The observed decrease in our case can also be related to the partial closure of the canopy, and to the infilling of the two main former canopy gaps with young woody species (see below). The majority of the reduction in the herb cover was due to a few forest species (Galium odoratum, Mercurialis perennis) and to species (Urtica dioica, Dryopteris filix-mas, Athyrium filix-femina) that were abundant in gaps in 1996. Hence, we argue that stability of herbaceous cover might be realized on a larger study site, where all stand dynamic situations are represented in natural proportions. We also assume that the extremely dry years of 2012 and 2013 as a potential climate factor could have contributed to the observed decrease.
The majority of the temperate forest resurvey studies document high species turnover, as a consequence of environmental change. Denser canopies effect not only the cover of the herbaceous layer, but also trigger trait-based reactions as the decrease in light-demanding species and increase in species with good abilities for living in shade (Brewer 1980, Hédl et al. 2010, Heinrichs et al. 2014. A direct trend in species trait shift can be observed in spite of the formerly assumed or postulated stability of a site (Brewer 1980, Taverna et al. 2005, Łysik 2008). Despite the expected stability, results of Taverna et al. (2005), for example, showed systematic changes in herbaceous species richness in a climax hardwood forest in North-Carolina, as a potential result of the elimination of ground fires and widespread grazing typical at the beginning of the twentieth century. Invasive species can significantly change the composition of the herbaceous layer even in primeval forests (Łysik 2008). One of the few examples, where the ground vegetation showed an overall stability, is reported from a Danish beech forest by Carøe et al. (2000), but the temporal scale was very short (5 years). Similar stability of species pool was found by Woods et al. (2012) in unmanaged northern hardwood forests over a three-decade period. These results highlight the rarity of the observed low species turnover of Kékes Forest Reserve after almost two decades, where all extinct species were represented with very low abundances in 1996; similarly, all colonization happened with low abundance and many of the colonizing species were present in the species pool of the close surroundings.
Our results support the concept that the experienced intensive fine-scale dynamics is profoundly governed by the stand dynamic events of the studied forest stand. Significant differences were found between the individual SDSs and between years in herbaceous cover, species richness and species turnover at the SDS-scale. Highest species richness, mean species number/plot and highest herb cover characterized the two main gaps (OCO and OBG) in 1996. An expressed decrease in the herbaceous layer cover was observed after 17 years in these two, gradually closing gaps, where the saplings have overgrown and started to cast shadow on the herbaceous plants. Lateral canopy expansion of bordering trees has also contributed to gap closure. A very high number of extinctions with high absolute species turnover were detected in these situations. In the old collapse (OCO) the cover of Dryopteris filix-mas, Athyrium filix-femina, and Urtica dioica was drastically reduced. In the old beech gap (OBG) pronounced recession of Galium odoratum, Dryopteris filix-mas, Fagus sylvatica (young beech individuals have grown out of the herbaceous layer) and Mercurialis perennis was detected. Competition between the herbaceous plants and the saplings, especially with the very competitive beech saplings, and the resulting drop in herbaceous plant cover, was observed by several authors (Davis et al. 1998, Łysik 2008, Ujházy et al. 2007. On the other hand, the two younger gaps, created in 2010 and 2012 (3YG, 1YG) had low species richness, low mean species/plot and low herbaceous cover at the beginning, as the canopy was closed in 1996. The resurvey detected a sharp increase in both species richness and cover in the case of the three-year-old gap (3YG), where the massive colonization process was performed by gap specialists (Urtica dioica, Solanum dulcamara), by responsive forest species (Dryopteris filix-mas, Galium odoratum) and by tree seedlings (Fagus sylvatica, Acer pseudoplatanus). All these changes were more moderate in the case of 1YG, due to the shorter time-span. Despite the high absolute species turnover, the cover of the herbaceous layer and species richness changed only slightly in the control SDS. These results support our second hypothesis that stand dynamic events govern the herb cover, and drive the herbaceous layer changes, at least during this nearly two decades in our study site.
The Slovakian fir-beech forest investigated by Ujházy et al. (2005Ujházy et al. ( , 2007) had a very similar species pool in the herbaceous layer as our studied site. They found relatively low species turnover -compared to commercial forests -in the three developmental stages; rather the abundance or dominance values differentiated the growing, optimum and decay stages (Ujházy et al. 2005). The investigated short-term (4 year) changes were attributable to stand dynamics; they recorded a rapid change in the herb layer during the four-year period established in young gaps (Ujházy et al. 2007). However, they did not analyse colonisation and extinction patterns for all species.
Although the four subsamples of each plot were relocated only with a 1 m accuracy in 2013, statistical tests proved that the time effect (change in vegetation) was significantly higher than the relocation effect. With this we managed to assure that the interpreted temporal changes were not artefacts resulting from imperfect relocation.

Conclusion
We found that intensive local-scale extinction and colonization episodes were balanced at the stand-scale, resulting in overall stability of the species pool in the herbaceous layer vegetation. These results have important implications for both forest management and conservation.
The observed stand dynamics driven changes of vegetation indicate that the use of small regeneration areas in forest management can prevent competitive ruderal species (e.g. Calamagrostis epigejos) from causing serious problems in regeneration , Kelemen et al. 2012.
From a conservation viewpoint, one of the most important species-related implications is that the long-term survival of closed forest specialists, including ancient forest indicators Game 1984, Hermy et al. 1999), could be guaranteed by management techniques mimicking natural small gap dynamics. This is important for several of those species that are seriously dispersal limited, meaning that once they disappear it is extremely hard, or it takes decades to centuries, to recolonize secondary forests (Kelemen et al. 2014). We assume that the observed stability in the species pool was secured by the primeval character of Kékes Forest Reserve. By primeval we refer to both structurally rich old-growth character and long forest continuity. The former is illustrated by the presence of all important forest developmental phases produced by natural stand dynamics (natman report 2002) and the high amount and diversity (size and decay phase) of deadwood (Christensen et al. 2005) and also of forest specialist organisms (e.g. Ódor 2000, Ódor and Standovár 2001, Heilmann-Clausen et al. 2014). Long continuity is convincingly documented by Czájlik (2009). Similar relationships between the richness in ancient forest indicator plants and that of other organism groups were shown, for example, by Hofmeister et al. (2014) for macrofungi.
Several studies documented that habitat specialist species occur in the highest numbers in forested landscapes with long continuity of forest structures and habitats (for a review see e.g. Nordén et al. 2014). This emphasizes the specific conservation role of even small remnants of virgin or ancient forest, as they secure the long-term survival of habitat specialists, and they might also serve as a source for dispersal from these remnants into the surrounding landscape.
In areas, where nature conservation is not the sole purpose of management, managers have started to apply retention forestry as a means of integrating conservation concerns into forest management. A thorough review (Fedrowitz et al. 2014) showed that retention cuts supported greater abundance and higher species richness of forest species than clear-cuts. This positive effect increases with the proportion of retained trees and time since harvest. However, they also showed that retention cuts cannot be a substitute for uncut forest because it had negative impacts on mostly closed forest specialist species.
All these suggest that for successful conservation of forest biodiversity we still need to preserve existing remnants of woodlands with high conservation value and also to apply a mix of management approaches in their immediate surrounding that support natural processes and the creation of important habitats.

Keywords
Gopherus flavomarginatus, burrow, plant cover, habitat, temperature, microclimate introduction Research on the ecology of ectothermic organisms has established the importance of vegetation structure for their microhabitat selection (Hertz et al. 1994, Vitt et al. 1997, Litzgus and Brooks 2000, Bryant et al. 2002. Changes in vegetation produce variations in other microhabitat attributes, like light intensity, wind speed, air and soil temperatures (Pringle et al. 2003). Variation in these features influences thermoregulatory behaviors and activity levels in ectotherms (Adams and Decarvalho 1984, Huey 1982, 1991, Huey and Kingsolver 1989, Webb et al. 2005, Turbill et al. 2011, resulting in greater impact on species that are thermally sensitive to changes in habitat structure (Walther et al. 2002, Pringle et al. 2003. Population ecology theory predicts that in a changing environment, a population can adapt to new conditions, migrate to a place that favors its survival, or become extinct (Pease et al. 1989, Pringle et al. 2003 if the species presents a capacity of dispersion and limited evolutionary responses (Allendorf and Luikart 2007). Long term studies have established relationships between changes in vegetation density and animal movements and extirpations of populations with small distributions (Fitch 1999, Pringle et al. 2003. For example, abundances of forest birds in New Hampshire decreased considerably over a period of 30 years causing local extinction of four species; the most important local factor affecting bird abundance was temporal change in forest vegetation structure (Holmes and Sherry 2001). Likewise, it has been reported that for Gopherus polyphemus in sites invaded by an introduced weed, tortoises avoided areas where weeds had formed a dense monoculture, suggesting that habitat selection increases isolating effect of habitat fragmentation on this tortoise (McCoy et al. 2013).
The Bolson tortoise, Gopherus flavomarginatus (Figure 1), is North America's least studied tortoise; it is considered as Vulnerable by IUCN Red List (2015), and has a geographical distribution restricted to the Mapimí Basin in the Mexican Chihuahuan Desert (Aguirre et al. 1984). This restricted distribution is likely due to specific habitat requirements (Aguirre et al. 1997), including constant temperatures and humidity levels provided by their burrows throughout the year, as G. flavomarginatus seems to have a limited thermoregulation capacity (Adest et al. 1989); Adult individuals of this species have a high fidelity to their burrow, spending about 95% of their life hibernating or aestivating within this structure, and remain only 5% outside of them during the summer season (Adest et al. 1989, Lovich and Daniels 2000, Daren-Riedle et al. 2008, and adult tortoises are unlikely to be naturally depredated (Treviño et al. 1995).
Therefore, if Bolson tortoise requires specific microclimatic conditions to inhabit burrows and survive, variations in microhabitat are expected to influence either their use or abandonment. An analysis of microhabitat variation is shown here in relation to the occupation of burrows of G. flavomarginatus. Our objectives included: 1) charac- terization of the environmental factors of air temperature, relative humidity, substrate temperature and pH; physical factors of width and height of burrows and 2) determine how these factors are related to plant cover and occupancy of burrows. This information can increase understanding of this species' response to variation among microhabitats, and support conservation efforts for this species.

Methods
The 100 hectare study site, Tortugas, is located in the south-central portion of the Mapimí Biosphere Reserve, in Mexico (26°00', 26°10'N and 104°10', 103°20'W; CONANP 2006) and within the region known as the Mapimí Basin (Figure 2  At Tortugas, we followed the monitoring protocol established by CONANP to find adult tortoise burrows (CONANP 2006). Burrow monitoring was performed during two consecutive days in summer (September, 2010). There was no rainfall before or after those two days. We classified burrow occupancy (active, inactive, or abandoned) based on measuring external characteristics according to Auffenberg and Franz (1982) and Cox et al. (1987). Accordingly, an active burrow shows foot or plastron prints on the access tunnel and the surrounding mound; the soil is loose, with little compaction. In an inactive burrow, no tortoise tracks are seen, and soil at the burrow's entrance and the mound looks compacted. Finally, an abandoned burrow entrance shows an accumulation of debris, such as branches, grass, cobwebs, and the soil of the mound is clearly compacted.
In every burrow, was measured microhabitat structure considering the variables width (W) and height (H) of the entrance and the substrate's pH 30 cm inside. Dataloggers (Datalogger USB-WK057, accuracy: ± 1.0) were used to measure environmental factors continuously, including air temperature (T a ) and relative humidity (RH) inside (30 cm depth) and outside (30 cm above surface) the burrow, except pH, all environmental data were recorded each hour for 24 hours; substrate temperature inside the burrow (T s ) was also recorded using dataloggers (in contact with the surface). Also, we measured plant cover (PC) using an ellipse area formula (π × a × b/4, where a = major axis and b = minor axis), within three meters of each burrow.
Discriminant analysis was used to determine which habitat and environmental factors differentiate burrows categorized by their occupancy status. Normality was not achieved (Kolmogorov-Smirnov tests; P ≤ 0.05) and we transformed the continuous data (W, H, pH, T a , RH, T s ) with the logarithmic formula (X´ = LOG10(X + 1)), and PC with the arcsine formula (X´= Arcsin√X), according to Zar (1999). A Post Hoc test (LSD) was performed to identify differences among the averages of the three status groups. Lastly, correlation and linear regression analyses were performed to quantify the relation between significant environmental factors and PC in the sampled burrows and was plotted temporal variation of temperature. All statistical analyses were made using STATIS-TICA 10.0 (StatSoft 2011) software and considered statistically significant at P ≤ 0.05.

List of abbreviations W
burrow width H burrow height T ai air temperature inside the burrow T ao air temperature outside the burrow RH i relative humidity inside the burrow RH o relative humidity outside the burrow T si substrate temperature inside the burrow PC plant cover LSD least significant difference d.f. degrees of freedom SD standard deviation

Results
We located and measured a total of 61 burrows at the Tortugas study site. There was significant difference in the T si among the three types of burrows (F = 32.40, d.f. = 2, 58, P < 0.001; Table 1). Post hoc analysis (LSD) showed that abandoned burrows had higher T si (x -= 31.1°C, SD = 5.24) than active (x -= 28.0°C, SD = 4.7) and inactive (x -= 27.0°C, SD= 3.8) ones (Table 1). Results of discriminant analysis were as follows: the first function was statistically significant (ᴧ = 0.241, x²= 76.74, d.f. = 18, P < 0.001; n = 61), while the second function was not (ᴧ = 0.942, x²= 3.25, d.f. = 8, P < 0.917; n = 61). The first function's autovalue analysis indicates that this function explains 97.9% of the variation in burrow activity status, where T si showed the higher scores (Table 2 and Figure 3). table 2. Discriminant canonical function 1 scores with relation to burrow entrance width (W), height (H), air temperature inside the burrow (T ai ), air temperature outside (T ao ), relative humidity inside (RH i ), relative humidity outside (RH o ), substrate temperature inside (T si ), plant cover (PC), and substrate pH.  An inverse relationship was observed between PC and T si (y = -0.2181x + 41.504), indicating that the higher the plant cover around the burrow, the lower the substrate temperature inside it (Figure 4). Correlation and determination coefficients were high (R = 0.98, R 2 = 0.96, respectively); plant cover around the burrows influences 96% the increase of substrate temperature inside the burrows. This relationship was found to be highly significant (F = 1408.949, d.f. = 1,59; P < 0.001), the temporal variation of temperature is shown in Figure 5. Adest et al. (1989) described individuals of this species emerging from their burrows at night as a response to increasing substrate temperatures inside their burrows (>34°C)  because of thermal delay (which it is described as the speed at which the temperature fluctuations penetrate the substrate (Körtner et al. 2008)). Also, they described that around 0700 hr at a distance of 15 cm inside the burrow, the substrate temperature is still above 31°C, while the temperature of an adult individual (7.3 kg) is below 30°C when beginning foraging bouts at the surface. These observations support a hypothesis that the substrate temperature inside G. flavomarginatus burrows influences its occupancy dynamics, increasing the possibility of abandonment when the substrate temperature inside these structures consistently is equal to or greater than 31°C.

Discussion
Our analyses provided evidence that an increase in substrate temperature inside the burrows and their consequent abandonment at our Tortugas study site was correlated with vegetation cover at a scale of 3 m. Aguirre et al. (1984) described that the presence of G. flavomarginatus burrows at the Mapimí Biosphere Reserve seemed to be related to shrub type (Larrea tridentata (Coville 1893), Prosopis glandulosa (Torrey 1827) and grasslands (Hilaria mutica (Bentham 1881)). McCoy et al. (2006) and Waddle et al. (2006) reported that habitat quality reduction was the apparent explanation for the increase of abandoned burrows in G. polyphemus population in Florida, USA. Additionally, Aresco and Guyer (1999) stated that land cover changes around G. polyphemus burrows can result in their abandonment in certain habitats. Similarly, Boglioli et al. (2000), Hermann et al. (2002), Jones and Dorr (2004), Baskaran et al. (2006), and Ashton et al. (2008) all described that for other species of Gopherus the presence of burrows is associated with the vegetation, and that the permanent abandonment of these burrows seems to happen as a response to unfavorable habitat conditions. Moreover, Huey (1982), Hertz et al. (1994), Vitt et al. (1997), and Bryant et al. (2002) mentioned that vegetation structure plays a key role in the activity, feeding, and distribution of some ectothermic organisms. These previous studies support our conclusion that the occupation of G. flavomarginatus' burrows are related to microhabitat structure, with vegetation cover being one of the main environmental factors that can affect habitat selection, this interaction of temperature and microhabitat is key to the species' survival. With predicted increasing temperatures as climate change effects become more pronounced in the deserts of America (Friggens 2012), this interaction will be critical over the coming years. On the other hand, we consider that relative humidity inside the burrows is important for Bolson tortoise. However, this variable showed low scores to discriminate activity status of the burrows and did not present significant differences when comparing between activity status; therefore, it was no possible to determine its influence on burrows occupation.
It is important to note that G. flavomarginatus might not have originated as a desert ecosystems species, they appeared toward the end of the Tertiary, so they could have spent more than 94% of their evolutionary history during the Quaternary (Pleistocene-Holocene) living in non-desert grasslands (Van Devender and Burgess 1985). Therefore, their current restriction to cool microclimates in their summer burrows could be an extension of a physiology geared to a cooler, more mesic climate. Consequently, it is likely that their thermal physiology and even their social behavior reflect more their burrow microhabitat than the surface environment of the Chihuahuan Desert (Adest et al. 1989).