Hidden diversity in forest soils: Characterization and comparison of terrestrial flatworm’s communities in two national parks in Spain

Abstract Terrestrial flatworms (Platyhelminthes, Tricladida, and Geoplanidae) belong to what is known as cryptic soil fauna of humid forests and are animals not easily found or captured in traps. Nonetheless, they have been demonstrated to be good indicators of the conservation status of their habitat as well as a good model to reconstruct the recent and old events affecting biodiversity. This is mainly due to their delicate constitution, their dependence on the integrity of their habitat, and their very low dispersal capacity. At present, little is known about their communities, except for some studies performed in Brazil. In this work, we analyze for the first time in Europe terrestrial flatworm communities. We have selected two protected areas belonging to the Red Española de Parques Nacionales. Our aims include performing a first study of the species richness and community structure for European terrestrial planarian species at regional and local scale. We evaluate the effect of type of forests in the community composition and flatworms’ abundance, but also have into account the phylogenetic framework (never considered in previous studies) analyzed based on molecular data. We find differences in the species composition among parks, with an astonishingly high diversity of endemic species in the Parque Nacional de Picos de Europa and an extremely low diversity of species in the Parque Nacional de Ordesa y Monte Perdido. These divergent patterns cannot be attributed to differences in physical variables, and in addition, the analyses of their phylogenetic relationships and, for a few species, their genetic structure, point to a more probable historical explanation.


| 7387
ÁLVAREZ-PRESAS Et AL. mean annual precipitation, or type of forest cover, can result in variations of abundance and species composition in the soils at different scales, from microsite, site, local to regional levels (Melguizo-Ruiz, Verdeny-Vilalta, Arnedo, & Moya-Laraño, 2012). On the other hand, taxa from soil communities can exhibit also the genetic imprint of ancient climatic and geographic events that may have been lost in other organisms with higher dispersal capacity (Pfenninger & Posada, 2002;Sunnucks et al., 2006). As a consequence, these groups of organisms allow the reconstruction of old events affecting the generation and maintenance of biodiversity and become excellent indicators of the conservation status of forest soils. However, to be used for such aim, an extensive knowledge about the state and functioning of their communities is necessary. Land planarians (Platyhelminthes, Tricladida, Geoplanidae) belong to this group of animals; they are inhabitants of humid forests soils and top predators of other invertebrates. They are simple animals that do not possess mechanisms for water retention; therefore, they are dependent on soil moisture to maintain their water requirements and use vertical migration through soil, litter, and vegetation to keep their humidity (Winsor, Johns, & Yeates, 1998). Land planarians are in general sensible to disturbed habitats, although some are reported to be adapted to inhabit them (Carbayo, Leal-Zanchet, & Vieira, 2002;Oliveira et al., 2014;Álvarez-Presas, Amaral, Carbayo, Leal-Zanchet, & Riutort, 2015). Based on these features, some studies have highlighted the value of this group of organisms as bioindicators in relation to the habitat perturbations caused by human activities (Sluys, 1998).
Studies conducted in Brazil have shown that the community structure for terrestrial flatworms can be influenced by the vegetation type and by the degree of anthropic alteration (Carbayo et al., 2002;Fick et al., 2006;Fonseca et al., 2009). However, studies analyzing the effect of environmental factors have not found any of them as driver of the abundance or species composition of terrestrial planarians communities (Antunes, Leal-Zanchet, & Fonseca, 2012;Baptista & Leal-Zanchet, 2010;Boag, Jones et al., 1998;Boag, Yeated et al., 1998;Fick et al., 2006;Sluys, 1998;Winsor et al., 1998) with the only exception of pH and organic matter (see Baptista & Leal-Zanchet, 2010). At last, none of the community studies conducted so far has taken into account the phylogenetic relationships between the species found, nor whether these relationships or the climatic and geological history of the area can explain communities' composition differences among areas.
Moreover, none has used molecular data in conjunction with phylogenetic inference methods to delineate genetic lineages, which also provides a more accurate delimitation of molecular operational taxonomic units (MOTUs) or species and avoids the problem of identifying morphologically cryptic or pseudocryptic species (common in terrestrial planarians, Álvarez-Presas et al., 2015;Carbayo, Álvarez-Presas, Jones, & Riutort, 2016;Amaral et al., 2018) that may take much time and lengthen the process of community study.
In Europe, although the diversity of species is still much lower than in tropical regions, recent publications have shown that the species richness of this group in this temperate region is higher than previously suspected (Mateos, Sluys, Riutort, & Álvarez-Presas, 2017;Sluys, Mateos, Riutort, & Álvarez-Presas, 2016;and references therein). In this continent, all the native species belong to the subfamily Microplaninae and to a single genus, Microplana; although there are some doubts about whether the genus Rhyncodemus can also be autochthonous (Jones, 1988(Jones, , 1998Ogren & Kawakatsu, 1998). The European species are in general much smaller than the tropical ones and less colorful, which renders them less prone to be found and identified in any soil community study. Moreover, the fact that this group of animals belong to cryptic soil fauna, due to the difficulty in sampling them by the usual methodologies as using traps or soil fauna extractors, makes difficult their inclusion in any study of communities, nonetheless, they can contribute important information.
In Spain, the forests suitable for terrestrial planarians are few and are located mainly in the north of the Peninsula. In that region, three national parks bear forests with the characteristics needed to host terrestrial planarians: in the Pyrenees, the national parks Aigüestortes i Sant Maurici and Ordesa y Monte Perdido, and in the Cantabrian Mountains Picos de Europa. The two latter present a broader extension of these types of forests and are where we have focused our study. The Red Española de Parques Nacionales is an integrated system intended to protect and manage a selection of the best samples of the Spanish Natural Patrimony (http://www. mapama.gob.es/es/red-parques-nacionales/ last visited January 2018). Among its objectives, one of the most important is the protection and management of its biodiversity in order to ensure the proper functioning of ecosystems. It is desirable therefore that they house a biological diversity that is representative of its original biodiversity (Araújo, Lobo, & Moreno, 2007), that includes the levels of genetic diversity needed for the maintenance of its populations and that it is also representative of the diversity in the region. But for this it is necessary to know what was the original situation of the fauna inhabiting them, knowledge that is currently inexistent for the terrestrial planarian communities, as well as for other cryptic forest soil dwellers.
In this work, we performed a study on the diversity of terrestrial planarian communities focusing on two national parks of the Red Española de Parques Nacionales: Ordesa y Monte Perdido (hereafter referred as Ordesa) and Picos de Europa (hereafter referred as Picos).
These parks have been selected because they are located in the area of the Iberian Peninsula with the highest probability of housing terrestrial planarians, as explained above. The broader extent of Picos and its higher diversity of forests may influence the genetic diversity distribution between and within the parks, predicting finding a higher diversity and species richness in Picos. In addition, the two parks are situated in regions that have gone through different ancient climatic events (Hewitt, 1999;Petit, Brewer, Bordács, & Burg, 2002;Petit et al., 2003) also pointing to an expectation of a higher genetic diversity in Picos than in Ordesa, which could be reflected in a higher number of species and/or within species diversity. Thus, the specific aims of the study are as follows: (i) performing a first analysis of the species richness and the community structure for European terrestrial planarian species at a regional scale; (ii) analyze the effect of the forest type in their communities (local scale); and (iii) explore the drivers of the communities composition in the parks under a phylogenetic framework.

| Study area
The study area comprises two national parks in Northern Spain Picos extends across Asturias, León, and Cantabria provinces, in the Cantabrian mountain range. It has a surface area of 67,127 ha, with a pronounced influence of Oceanic climate (Atlantic climate), with cool summers and comparatively warm winters (Felicísimo, 1994). It has the highest limestone formation in Atlantic Europe, with important karstic processes, chasms reaching more than 1,000 m, very clear glacial erosion and presence of lakes. It is characterized by a narrower range of annual temperatures than those encountered in Ordesa at comparable latitudes, lacking, for instance, the extremely dry summers typical of the other park, which is more influenced by Mediterranean climate, and with a mean annual precipitation between 1,109 and 1,968 mm/year. Unlike Ordesa, most rainfall on Picos comes as drizzle, and mountain fogs are very frequent due to the Oceanic influence (Felicísimo, 1994). Forested area occupies 25% of the park's range, being beech forests of Fagus sylvatica the most abundant forest type (18.4% of the park area), followed by mixed and oak forests (<0.5%). Atlantic mixed forests of Picos, relics difficult to find in Spain, appear on the lower part of the mountain and are interspersed with meadows areas. Common and cornish oaks (Quercus robur or Quercus petraea) and hazels (Corylus avellana) are intermingled with birch (Betula celtiberica) as the main tree species (5%), maples (Acer sp.), linden (Tillia sp.), ash (Fraxinus excelsior), chestnut (Castanea sp.), and walnut (Juglans regia) trees; at its feet, an undergrowth of brambles, briars, and thorn. Our collection sites were established at elevations ranging from 115 to 1,353 m a.s.l. distributed in the three main areas composing the park (western, central, and eastern massifs).

| Collection sites and sampling protocol
When working with cryptic soil fauna, one important question is the establishment of a standard sampling unit and the sampling methodology used in order to get comparable data between these sampling F I G U R E 1 Map showing the situation within Spain of the two National Parks, and a detail of the distribution of plots within each park. Picos and Ordesa map at the same scale. Squares, oak forest; triangles, mixed forest; circles, beech forest; asterisks, pine forest units. The authors working on terrestrial flatworm's ecology (see Baptista & Leal-Zanchet, 2010;Baptista et al., 2006;Carbayo et al., 2002;De Castro & Leal-zanchet, 2005;Fick et al., 2006;Leal-Zanchet et al., 2011) used plots of fixed sizes (50 × 2 m, 2 × 2 m, 7 × 7 m) or transects of fixed length (30 to 50 m) as a sampling unit.
On each plot or transect, the sampling methodology used was always direct searching by a fixed number of researchers (between one and five) during a standardized lapse of time (between 15 min and 1 hr). In our study, each sampling unit consisted of a forest plot to 03/7/2014 in Picos. Every sampling day, plots of different forests type were sampled in order to avoid temporary self-replication. The data of the two sampling campaigns were pooled.
In the same day of the collections, the animals were visualized under a stereomicroscope, their morphological external appearance recorded and photographed, and finally fixed. When the specimens were big enough, a small anterior section was fixed in absolute ethanol for DNA extraction, and the rest (including the parts necessary for histological studies) were fixed with Steinmann fluid and stored in 70% alcohol in order to study the copulatory apparatus and other structures that allow the diagnosis of the species.

| Environmental parameters
Pluviometry data were obtained from two meteorological stations located at the east and west ends of each national park. Total annual rainfall as well as accumulated rainfall of the 3 months previous to each sampling campaign was recorded. For calculation purposes, the mean value of the two stations of each park was considered. In Ordesa, the two meteorological stations were located in Bielsa municipality (at 1,100 m asl) and Torla municipality (1,000 m asl). In Picos, the two meteorological stations were located in Soto de Sajambre municipality (1,500 m asl) and Sotres municipality (1,200 m asl). On each sampling campaign soil temperature, pH and water content were measured. On each plot, soil temperature was taken in three points and the mean value was used for calculations. Also on each plot, two samples of soil substrate were taken to measure hydric content and two other samples were taken to measure the pH in the laboratory; the mean values of each variable were used for calculations. Measures of pH were performed the same day with a portable pH meter (PH25, Crison) by diluting the soil substrate samples in 5 volumes of distilled water. Also the same day, the two soil samples for hydric content were weighted to obtain their fresh weight (Fw). Once back in the University laboratory, the samples were dried in a stove at  Genes Cox1 and EF were aligned based on the amino acid sequences using the Clustal W plugin included in the BioEdit software. 7.0.9.0. (Hall, 1999). Ribosomal RNA gene sequences were aligned using the online version of the software Mafft v. 7 (Katoh & Standley, 2013) applying the G-INS-i iterative refinement method and, subsequently, checked the alignments by eye with Bioedit. Positions that could not be unambiguously aligned were subsequently excluded from the analyses by applying GBlocks v 0.91b (Talavera & Castresana, 2007), with half allowed gap positions and a minimum length of a block of 10 nucleotides, to obtain the maximum number of nucleotides. Based on these alignments, we estimated the DNA sequence evolution model that better fits the data by using jModelTest v. 2.1.4. (Darriba, Taboada, Doallo, & Posada, 2012), applying the Akaike information criterion (AIC).

| Species assignation and phylogenetic inference
For the Cox1 gene, a saturation test was run using the software DAMBE6 (Xia, 2017) by plotting observed transitions and transversions vs. gene divergences under the GTR model.

Two datasets have been used for different analyses: (i) Cox1
dataset, included all Cox1 mitochondrial gene sequences to assign individuals from Picos and Ordesa to MOTUs; (ii) concatenated dataset, included the information of the three nuclear genes to infer a general phylogeny. Both alignments included sequences downloaded from GenBank (see Supplementary Information Table S1 for the individuals included in each dataset).
Phylogenies of the two datasets were inferred applying the maximum-likelihood method (ML) with the software RaxML v.8 (Stamatakis, 2014) estimating bootstrap support values from 10,000 replicates. We also used the Bayesian inference method (BI) as implemented in the software MrBayes v3.2.2. (Ronquist et al., 2012) for the concatenated dataset. Two runs were applied producing 5 million generations for each and of these a tree each 1,000 was stored.
A 25% default burn-in was used, likelihood values (log-likelihood) of cold chains were checked to have reached stationarity, and the convergence of the two runs was verified by the average standard deviation of split frequencies being ≪0.01. A consensus tree from the remaining trees was obtained. For the Cox1 dataset, we used BEAST software v2.3.2 (Bouckaert et al., 2014) since with MrBayes, it was impossible to reach the convergence of the two runs. A relaxed log normal clock was applied, and 100 million generations were run storing a tree of every 10,000.
For the MOTUS shared by the two parks and with a higher abundance (M02 and M78), haplotype networks were constructed using the median-joining method (Bandelt, Forster, & Röhl, 1999) as implemented in the Popart program (available at http://popart.otago.ac.nz, last visited November 2017). To perform these haplotype networks, information from the Cox1 dataset was used, adding sequences from a previous analysis of Microplana terrestris in the Iberian peninsula (Álvarez-Presas, Mateos, Vila-Farré, Sluys, & Riutort, 2012) in the case of M02, and adding known sequences from M. aixandrei from the GenBank database, in the case of M78 (see Supporting information Table S1).

| Numerical methods
On each plot, the measure unit in numerical analyses was the terrestrial flatworm abundance, defined as the number of specimens collected by two researchers during one hour in the two sampling campaigns pooled. We also quantified the species richness (number of species or MOTUS) per sampling plot. These two variables were compared among the parks and the forest types by ANOVA or t-Student tests (with Tukey test post hoc comparisons when necessary) after checking the homogeneity of variances by the Levene test (Levene, 1960). A Pearson correlation test between soil water content, temperature, and pH with respect to flatworm abundance and richness per plot was performed. These tests were performed using R-language package (Team, 2003).
From the matrix of abundance data of flatworm species (eliminating species with overall abundance inferior to three specimens), we carried out an analysis of similarity in flatworm species composition between national parks and among forest types. To do this analysis, the abundance data were pooled by forest type and transformed using log(x + 1) procedure. Differences in flatworm species composition between pairs of sampling points were quantified by the Bray-Curtis similarity index. From the similarity matrix, we performed a similarity analysis (ANOSIM, PRIMER-E 2001), which gives a general R-value and allows pairwise comparisons between the areas compared (parks or forest types).
For each national park, we estimated the three basic components of the diversity (sensu Whittaker, 1960). We define α-diversity as the diversity measured on each forest type, β-diversity as the species turnover between the forest types of each national park, and γ-diversity as the diversity value measured on each national park as a whole (pooling the forest types). Estimates of α-and γ-diversity were done using Shannon-Weiner diversity index (H), and species richness. β diversity was estimated using the Sørensen dissimilarity index between pairs of forests type on each national park; this index, defined as β t in Wilson and Shmida (1984), is bound between 0 and 1, where 0 means that the two sites have the same composition (i.e., they share all the species), and 1 means that the two sites do not share any species. Commands "diversity" (with base 2 logarithms) and "betadiver," from R Vegan package, were used for diversity calculations.
We have used species accumulation curves to model the species increase in relation to sample size at two different scale levels. First, analyzing the two national parks as a whole and second analyzing only beech forests of the two parks (as this is the shared forest type in the two parks). Due to the disparity in the sampled area (Picos is 4.3 times higher than Ordesa), in order to make comparable the species accumulation curves between beech forests in the two parks, the second analysis has been performed constructing four different plot data matrices in Picos, each of them with surface area equivalent to that found in Ordesa (see Supporting information Figure   S1): north (12 plots, A + B), south (12 plots, C + D), east (eight plots, B + D), and west (16 plots, A + C). For the two species accumulation curves analyses we used the "method = random" from the command "specaccum" in R Vegan package, calculating the mean species accumulation curves and its standard deviation from random permutations of the sampling units (plots).
In general, in sampling protocols, not all species are detected in any site, and these unseen species also belong to the species pool (Oksanen, 2016). In order to estimate the number of unseen species on each forest type and on each national park as a whole, command "specpool" of R Vegan package was used. Command "specpool" studies a collection of sites and assumes that the number of unseen species is related to the number of rare species, or species seen only once or twice (see Oksanen, 2016 for a detailed discussion).
Command "specpool" implements several nonparametric estimators of which we selected the Chao 1 estimator (Chao, 1987;Chiu, Wang, Walther, & Chao, 2014), the first order jackknife estimator, and the bootstrap estimator (Smith & van Belle, 1984), calculating the mean estimate values, and associate variances. In Ordesa, we have found four MOTUs all of them having a known distribution through Europe or wide in the Iberian Peninsula, not being in any case endemic from the Park or the area (Table 1). On the other hand, in Picos, seven of the 13 MOTUs found are endemic from the Park, three are found only in wet forests in North-Western Spain, and three have a wide distribution through Europe (Table 1).

| Phylogenetic relationships
The analyses of the Cox1 gene phylogeny showed no resolution for the relationships among MOTUs, and the saturation test shows this molecule to be effectively saturated (Supporting information Figure   S2); for this reason, only nuclear genes were used to infer the relationships among MOTUs. Figure 3 and Supporting Information Figure S3 show the phylogenetic trees inferred from the concatenated dataset by maximum likelihood and Bayesian inference, respectively; there are small differences between both trees affecting nodes with low support but some clear relationships appear. Both trees show a clade constituted exclusively by MOTUS coming from Picos that is highly supported in BI but not in ML (0.93 PP and 50% BP) (A in Figure 3). This group is sister to a clade (B) constituted by M. terrestris (M02) and M. cf. aixandrei-2 (M78); the two species present in the two parks and also in other regions around Europe. M35, a species exclusive from Picos, in the BI analyses is sister to clade B, while in the ML tree is sister to clades A and B but with low support in both analyses.

| Abundance, diversity, and community composition at regional scale (parks)
Given the low general abundance of the MOTUs, we have pooled the data from spring and autumn samplings for the following analyses. The proportion of plots without any terrestrial planarian was higher in Picos (34%, 16 of 47 plots) than in Ordesa (11%, five of 46 plots) (Supporting information Table S2). The abundance per MOTU in each plot was in general low with only a few MOTUs reaching values close or over 10 animals in some plots, resulting in mean abundances per type of forest below 1 (Table 1) with a few exceptions (M02, M28, and M78). The distribution of MOTUs is not uniform through parks (Table 1) as demonstrated by the ANOSIM analysis that shows significant differences in species composition between parks (p level < 0.05, Table 2).
Mean abundances and richness per plot for each park were not significantly different (Table 3). Nonetheless, the total number of species and its diversity was different among parks, presenting Picos a higher richness (13 vs. 4 species, Table 4). The species accumulation models performed with all plots of the two parks ( Figure 5) showed that, with around 12 accumulated plots, the predicted species number is higher in Picos than in Ordesa. Also in Ordesa, the in predicted species number in Picos, but no increment in Ordesa (Table 4). The species accumulation models performed only with beech forests plots of the two parks also show that beech forest in Picos host more species than Ordesa irrespective of the area selected (Supporting information Figure S1). At last, Gamma diversity measured with the Shannon diversity index (H) was higher in Picos than in Ordesa (Table 4).

| Abundance, diversity, and community composition by forest type
In the two forest types analyzed in Ordesa, beech, and pine forests, the proportion of plots with terrestrial planarians was high (Table 1).
In both forests, the same four species and in the same proportion were found (Figure 6a), and M02 was the most common species,  (Table 2). In this park, mean terrestrial flatworm's abundance per plot was similar in the two forest types, while mean richness per plot was significantly higher in pine than in beech forest (Table 3). Total species number (S) was identical in the two forest types (four species), and α-diversity values measured with Shannon index (H) were also very similar (Table 4).
The Sørensen dissimilarity index value (β t ) measured between beech and pine forests in Ordesa was 0.00, meaning that the same species were found in both forest types. As a result, the nonparametric total species estimators showed the same values as actual total species number for beech and pine forests (four species, with only a little deviation for bootstrap estimator, Table 4).
In Picos, three forest types were analyzed, beech forest (11 species of terrestrial flatworms), oak forest (six species), and mixed forest (10 species), the two-first having a lower proportion of plots with terrestrial planarians (Table 1). In beech and oak forests, the abundance of different species was similar, but in mixed forest, M02 was by far the most abundant species (Figure 6b). ANOSIM analysis (Table 2) showed no significant difference in species composition between beech and oak forests (p value > 0.05), but significant differences were found between terrestrial flatworm communities of mixed forest with respect to beech and oak forests (p values < 0.05).
Mean terrestrial flatworm's abundance per plot was significantly higher in mixed forest than in beech and oak forests, while mean species richness per plot was not significantly different between these three forest types (Table 3). Shannon diversity values (H) were higher in beech and oak forests than in mixed forest due to a higher abundance of M02 and M28 in the last forest (Table 4, Figure 6b).
The mean Sørensen dissimilarity index value (β t ) measured between the three forest types was 0.38 (β t between beech and mixed forests = 0.24, between beech and oak forests = 0.41, between mixed and oak forests = 0.50). The three nonparametric species number estimators (Chao1, Jack1 and Boot) showed an increment in total species number in the three forest types compared with actual total species number, especially for oak forest and mixed forest (Table 4).

| Environmental parameters
Total annual rainfall and accumulated rainfall of the three previous months to each sampling campaign are shown in Supporting information Table S3. In Ordesa, the mean total annual rainfall for the years of this study was within the normal range of historical data (mean annual rainfall range between 1,129 and 1,690 mm/year).
In Picos, however, rainfall was very low during the sampling period (especially in 2013), far below the mean annual precipitation range normal for the park (between 1,109 and 1,968 mm/year) and, which is more important to soil humidity, the accumulated rainfall during the three previous months to sampling dates was also very low (also especially in 2013).
The mean soil water content, temperature, and pH values measured in the two sampling campaigns on every forest type of each TA B L E 2 Differences in species composition (ANOSIM results) between parks (regional scale) and between forest types (local scale). Data from fall and spring pooled TA B L E 3 Species abundance and richness per plot on each park (regional scale) and forest type within parks (local scale) and Student's t test or ANOVA for the comparisons of abundance and richness national park are shown in Supporting information Table S4. As the same plots were visited during the two sampling campaigns, altitude measures were the same in the two dates (Supporting information   Table S4). In Ordesa, in general, beech forest soils showed higher mean water content and pH values and lower temperature than pine forest soils. Comparing the two sampling campaign' values, the main differences were the higher mean soil water content and the lower mean soil temperature values in 2014 (especially in beech forest). In Picos, generally beech forest soils showed higher mean water content and lower temperature than oak and mixed soil forests, while mean pH value was higher in mixed forest than in beech and oak forests. Comparing the two sampling campaigns values, the main differences were the higher soil water content and temperature, and lower pH values in 2014.
In the two parks, no significant correlation was obtained between flatworm abundance and richness per plot with respect to soil water content and temperature. In Picos, total terrestrial flatworm abundance and richness per plot were significantly correlated with soil pH values (Supporting information Table S2), presenting higher values of abundance and richness those plots with moderate-to-high pH values (Supporting information Figure S4), while in Ordesa, these correlations were not significant.

| D ISCUSS I ON
This study presents the first community analysis of terrestrial planarians in Europe. The analysis, performed in two national parks in northern Spain, has revealed an unexpected diversity for this group of animals. However, this diversity is not equivalent in both parks; community composition and total amount of species significantly differ. Moreover, within parks, only Picos shows some significant differences among forest types. In the following sections, we analyze the reasons to explain such differences. TA B L E 4 Species diversity per parks (regional scale) and per forest type within parks (local scale)

| Differences in flatworm's abundance and richness between forests within parks
F I G U R E 5 Flatworm species accumulation curves per park.
Vertical lines indicate standard deviations and prey abundance (Boag, Jones et al., 1998;Boag, Yeated et al., 1998;Fick et al., 2006;Johns et al., 1998;Winsor et al., 1998), but none of them point to any of these factors as the main driver of the presence/absence of these terrestrial invertebrates. Antunes et al. (2012), working in an undisturbed area of Araucaria forest, found that terrestrial flatworms were not significantly associated with any particular microhabitat condition, which included environmental parameters as soil humidity, but also leaf-litter and other fauna composition. Comparing land flatworm communities in two types of forests in Southern Brazil, Fick et al. (2006) found that pH differences together with thermal amplitudes may explain dissimilarities in composition between the communities due to disparities in pH preferences for different species.
In the present study, flatworm abundance and richness per plot showed a significant correlation only with pH in Picos (Supporting information Figure S4), with higher values of these parameters for neutral to basic soils (pH values between 6.5 and 7.5). This could be the explanation of the significant differences found for one of the forest types in Picos, the mixed forest, that seems to harbor higher density of planarians per square meter (individuals per plot), and shows a significantly different species composition with respect to beech and oak forests.
In Ordesa, although the total abundance per plot was equivalent between forests, we found a higher species richness in pine forests. Both forest types in Ordesa showed mean pH values over 6 (Supporting information Table S4), around what appears to be optimal for Microplana species attending to the results in Picos (Supporting information Figure S3), which may be one of the factors to explain the lack of differences between beech and pine forests in this park.

| Differences in flatworm's abundance and richness at a regional scale
Our results showed the existence of significant differences between parks in terms of species composition and richness ( Table 2) that cannot be explained by differences in the carrying capacity between parks (Table 3), nor by the inequality in the area of the two parks (Supporting information Figure S1). The environmental factors analyzed also fail to explain such differences. However, the historical analysis based on the phylogeny of the species and the genetic structure within the shared species provides some clues about the differences observed.
Of the two shared species between parks, M02, (M. terrestris, Figure 4a), presented a highly structured haplotype network with no shared haplotypes between parks. In a previous work  (Hewitt, 1999;Petit et al., 2002). As could be ex-  Mateos et al., 2017) and in this case, the different regions share at least one common haplotype (Figure 4b). This could be indicating a relatively recent expansion, although we cannot know, based on the available data, what was the origin and direction of this colonization, whether from north to south or from south to north, and also if the origin is in the rest of Europe or in the peninsula. We cannot discard, also, the possibility of human transport, although the characteristics of the species and the habitat it occupies (they are very small, white animals that hide under rocks and rotten trunks in humid forests) make it unlikely.
Nonetheless, in both parks, we find species with high genetic diversity and occupying basal positions in the phylogenetic tree, as M71 in Picos and M01 and M25 in Ordesa. These species may represent old lineages with a wide distribution through Europe before the LGM that may have recovered from different refugia throughout the continent.
In summary, historical factors offer plausible explanations to be considered the drivers of the present differences among regions in the Iberian Peninsula. Paulo and Florianopolis only in Brazil (Sluys, 1999). These figures simply confirm the known rule that species richness is much higher in the tropics than in temperate and cold regions (Brown, 2014), being terrestrial planarians one more example to certify this fact.

| Final remarks
Nevertheless, our results confirm another rule stating that there is a very large bias in species descriptions according to the number of taxonomists dedicated to the group: Since we started working on species identification in Europe, we have seen that the number of species has increased rapidly.
We aimed to analyze the terrestrial planarian community composition in two protected areas to have a view of the original fauna situation in front of what can be found out of those areas and to have a first approximation of the value of protected areas for soil communities. Our results allow us to certify that in both parks, terrestrial planarians fauna is representative of the diversity found in their area, and in the case of Picos, it certainly protects a highly endemic and at the same time genetically diverse and genealogically old group of species, which completely fulfills the objectives of the park. Our results also point to the importance of the study of soil fauna, generally poorly considered when planning protected areas management or new protected areas.

ACK N OWLED G M ENTS
We thank Amparo Mora Cabello de Alba (Parque Nacional Picos de Europa), Elena Villagrasa, and Ramón Castillo (Parque Nacional Ordesa y Monte Perdido) for providing us with the infrastructure and logistical support necessary to carry out the sampling campaigns. We also acknowledge all the guards of both national parks that have accompanied us during our field work days. We thank Laia Leria, Paula Escuer, and Arnau Poch for help with fieldwork. Two reviewers (Ana Maria Leal-Zanchet and Lisandro Negrete) are gratefully acknowledged for their helpful comments. This research was financed by the Ministerio de Agricultura, Alimentación y Medio Ambiente (Spain), within the program "Ayudas a la investigación en Parques Nacionales" (ref. 589, 2012).

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
The three authors contributed to the design of the project. MAP and EM performed the samplings. MAP was in charge of the laboratory work and analysis of molecular data. EM performed numeric and statistical analyses on communities. The three authors contributed to final analyses, discussion of results, and writing of the ms.