Variability of diatom community composition and structure in mountain streams

Small rivers support high levels of biodiversity, being especially sensitive to the effects of global change. Temporal records of community composition in minimally impaired streams can be used to explore trends in biodiversity in response to climate change and natural temporal variation. We approached the comparison of two time periods (2003–2008 and 2016–2020) to study whether the composition of diatom assemblages changed over time in twenty-three streams of the mountain range of Picos de Europa (Northern Spain). The stream’s water chemistry indicated significant decreases in N_NO3− and P_PO43− content over time. In these minimally disturbed streams, the specific diatom community was dominated by Achnanthidium pyrenaicum, Achnanthidium minutissimum and Cocconeis euglypta. PERMANOVA analyses did not identify significant changes in diatom assemblage composition between periods or river types. Diatom indices (e.g. IPS, NORTIdiat) indicated high or good ecological status and relatively high alpha diversity values were found in these mountain rivers during the studied years. Although diversity and evenness showed a significant decrease over time, the temporal stability of the river-type diatom reference community between the two periods should be considered as an indicator of biodiversity persistence of high importance when monitoring the ecological status following the reference condition approach.


Introduction
One of the most important global challenges of the twenty-first century is to maintain the benefits that aquatic ecosystems bring to humans without affecting aquatic biodiversity and fundamental ecosystem processes (Pawlowski et al., 2018). To achieve this, extensive national and international regulations have been adopted to protect water resources all over the world, such as the European Union Water Framework Directive (WFD-European Commission 2000/60/ Handling editor: Judit Padisák Supplementary Information The online version contains supplementary material available at https://doi.org/10.1007/ s10750-021-04779-4. EC). Since WFD was implemented in Europe, it has provided an economic, social, and environmental approach to the sustainable exploitation of water bodies based on the protection of their biological and chemical quality at local, regional, national, and European levels (Cruz et al., 2010;Mortágua et al., 2019). This regulation requires that all water bodies including estuaries, coastal waters, rivers, lakes, and groundwaters reach the ''good ecological status'' at ambitious deadlines (Hering et al., 2018). In streams and rivers, aquatic organisms used as biological quality elements are benthic macroinvertebrates, macrophytes, phytobenthos and fish (WFD-European Commission 2000/60/EC, Karaouzas et al., 2019). Specifically, phytobenthos and invertebrates are the most widely used in rivers to assess anthropogenic pressures, due to their sensitivity to different types of stressors and their integration of pressures at different scales of time .
Diatoms constitute a dominant and diverse algae group in aquatic ecosystems, playing a key role in primary production, as well as in the silica cycle (Burliga & Kociolek, 2016). Its importance increases especially in small streams, where diatoms constitute the predominant primary producer component (Kireta et al., 2012); therefore, affecting changes in the structure and functioning of diatom communities to entire ecosystems (Peszek et al., 2021). They have been the subject of studies relating their response to anthropogenic stressors, nutrient enrichment, pollutant episodes (e.g. oil spills) and to already observed effects of climate change (Jiang et al., 2014;Polmear et al., 2015). Many are the physico-chemical alterations attributed to global change (e.g. increase in temperature, eutrophication; Jiang et al., 2014;Raimonet et al., 2018), being closely related to long-term changes in diatom communities (Jiang et al., 2014). Results up-to-date evidence that time series on biotic and abiotic data are useful in the analysis of long-term responses of individual species to the temporal variation of environmental conditions (Schlüter et al., 2012). However, not many studies have aimed to understand the natural variability that exists in diatom communities inhabiting minimally disturbed streams (Rott et al., 2006;Stoddard et al., 2006;Novais et al., 2020) or reference streams (Cantonati et al., 2020;Peszek et al., 2021;Rimet et al., 2004).
The reference conditions of WFD are established following specific and objective criteria and rules on the accepted level of human pressures to select a spatial network of minimally disturbed sites (Pardo et al., 2011(Pardo et al., , 2012. In fact, finding common approaches to define criteria for reference conditions was one of the great challenges of the WFD. The ''biological reference condition'' is based on measurements of biological parameters sampled at nonimpacted sites having minimum levels of anthropogenic pressure (Nõges et al., 2009). The reference conditions cannot be considered static but change because of the impacts of climate change on the physical and chemical conditions of water bodies (Nõges et al., 2007). Given that the effects of climate change can cause changes in the ecosystems at reference conditions, Nõges et al. (2009) stated that the WFD should re-evaluate the reference conditions according to the observed changes, using monitoring networks of potential reference sites or modelling tools that allow distinguishing the effects of climate change. Consequently, restoration objectives (i.e. good ecological status) should also be periodically assessed accepting that some changes (e.g. loss or changes in the distribution of some species) have become permanent.
The ''biological reference condition'' is linked to each stream type and type specific reference sites must present the full range of conditions expected to occur naturally in each stream type (Sánchez-Montoya et al., 2012;Stoddard et al., 2006). Streams and rivers of similar environmental and hydromorphological conditions and inhabited by the type of specific biological communities are classified within the same ''type''. Following WFD, river types are first differentiated by abiotic descriptors (e.g. altitude, catchment area, geology) and secondly, stream types are confirmed with their biological communities. Many reference sites correspond to small streams in nature-preserved areas where streams and rivers tend to be in good ecological status. As the water temperature has generally close values to that of groundwater at the source in headwater streams (Caissie, 2006) there are ideal study cases to understand ecological responses to water temperature gradients (Arai et al., 2015;Delgado et al., 2020). The study of headwater's biological communities should allow us to understand how climate change can affect freshwater ecosystems (Nukazawa et al., 2018). Despite that headwater streams represent [ 70% of the total river network length (Kuglerová et al 2017) they have been understudied since they are usually not included in water management due to their small size (Stubbington et al., 2017). Nevertheless, they constitute valued water resources and they are also very important because of their ecological processes, such as sources of migrants, reservoirs of biodiversity or nutrient absorption and retention (Hlúbiková et al., 2014;Peterson et al., 2001). These headwaters are particularly vulnerable to changes in land-use and pollutants from diffuse sources, as they have a greater contribution of watershed area to stream area compared with larger streams (Selby et al., 1985). Different land-uses nearby headwater streams have repercussions on water physical and chemical variables (e.g. variations in conductivity levels and nutrient content; Pérez et al., 2013), but also on the structure of biological communities (e.g. diatoms and invertebrates; Hlúbiková et al., 2014;Moore & Palmer, 2005).
Conservation of mountain streams diversity is believed to be most effective when management focuses on restoring or conserving environmental conditions at various points within relatively small watersheds (Göthe et al., 2014). Although, major efforts on biodiversity conservation are usually directed towards large size species, such as mammals or vertebrates (e.g. Grenyer et al., 2006), few studies dealt with diatom conservation status valued by the risk of species to disappear from an area (Falasco & Bona, 2011;Lange-Bertalot & Steindorf, 1996).
Climate change has caused in Spain in recent decades a 6-32% reduction in rainfall, an increase of 1.5-3.38C in temperature and a 2-54% reduction in water resources (Senent-Aparicio et al., 2017). The frequency of warm summers has increased in the Picos de Europa National Park (Serrano et al., 2011) and the effects of climate change are expected to worsen in the near future . However, in recent years this National Park has experienced a reduction in agricultural practices and an expansion of forest area (García-Llamas et al., 2019). Therefore, we aimed to test whether diatom assemblages changed over time in well-preserved mountain streams. We used data from diatom samples and environmental variables from streams in a mountainous protected area in the Cantabrian Cordillera, Northern Spain. Samples were collected throughout two periods spanning the last 17 years (2003-2008 and 2016-2020). Our working hypotheses are as follows: (i) the low variability of diatom communities over time during the two studied periods, since the mountain streams studied are located in a well-preserved area; (ii) the absence of differences in diatom composition amongst Spanish-established river types; and (iii) the maintenance of high/good ecological status over time according to diatom indices (e.g. diversity, IPS index (Cemagref, 1982) and NORTIdiat index (Pardo et al., 2018)).

Study area
The study area is in the Picos de Europa National Park (N Spain; Fig. 1). This is a high mountainous area up to 2,650 m above sea level (m. a. s. l.) located in the northern slope of the Cantabrian Mountains . The Park has been a UNESCO biosphere reserve since 2003 and it is part of the Natura 2000 network. The climate is humid and temperate with cool summers and no dry season , with averaged rainfall being approximately 1500 mm, with temperatures ranging between -3 and 17°C and with snow covering the territory above the 1500 m for 6 months a year Elicisimo, 1990). The substrate is mostly limestone, with well-developed karst and with glacial debris (Canseco & Heredia, 2003). The information on land-use activities was estimated from CORINE land cover maps (Bossard et al., 2000; Coordination of Information on the Environment, Land Cover 2000). The land-use activities did not experience remarkable changes over time in this protected area, although it seems that agriculture has decline over the last decade (Perpiña-Castillo et al., 2020). All the points showed between 85 and 100% of natural and semi-natural areas and between 0 and 15% of agricultural areas. The most common riparian species in these Cantabrian Mountains are Castanea sativa Mill., Fagus sylvatica L. and Pinus sylvestris L., which are widely distributed in the north of the Iberian Peninsula (José Vidal-Macua et al., 2017;Nunes et al., 2020). The most frequent understorey species were Plantago lanceolata L., Trifolium pratense L. and Trifolium repens L. (Prince et al., 2012).
Biotic groups that include meiofauna and diatoms are poorly studied in headwaters (Stubbington et al., 2017).

River typology
According to the WFD, the ecological assessment of biological communities has to be 'type specific', i.e. water bodies should be grouped according to their physical, biological and morphological attributes. River types can be defined using two different WFD systems: System A is a standardised system with physical descriptors and mandatory limits, whereas System B is more flexible and allows the inclusion of mandatory and optional variables and limits (WFD 2000/60/EC). Spanish rivers are classified according Fig. 1 Map of the study area in the Cantabrian region of Spain locating the 23 sampling sites (n = 46) in the Sella (SE codes), Deva and Cares rivers (DC codes). The Picos de Europa National Park has been highlighted in green to a System B typology based on physico-chemical descriptors, accordingly, within the Picos de Europa National Park; the rivers are divided into six different Spanish types (Table 1; RD 817/2015; ARM/2656/ 2008). On the other hand, the predictive diatom-based model NORTIdiat developed by Pardo et al. (2018) established river types based on their diatom communities in reference sites and then used mandatory and optional descriptors of the WFD typology to environmentally characterise them, corresponding to a typology B of 4 river types that represented existing gradients in calcareous geology, basin size and slope along the whole North of Spain (From Galicia to the Basque Country). The resulting river types were named accordingly: Type 1, calcareous mountain rivers; Type 2, main river axes; Type 3, siliceous rivers; and Type 4, mixed rivers. According to this typology, all streams analysed in this study correspond to Type 1, calcareous mountain rivers (Pardo et al., 2018).

Diatom sampling
A total of 46 diatom samples and corresponding environmental variables were collected at 23 sampling sites established along the Sella, Deva and Cares rivers during the summer (June-September) of 2003 (n = 14), 2004 (n = 3), 2008 (n = 6), 2016 (n = 1), 2017 (n = 1), 2018 (n = 6), 2019 (n = 11) and 2020 (n = 4). 14 out of 23 sampling sites studied were sampled in both periods. Between 2008 and 2016 no samples were collected. The sampling period and Spanish-type correspondence (Spanish typology) of the 46 samples are described in Table 1. These 46 samples belong to the monitoring network established in the implementation of the Water Framework Directive (WFD 2000/60/EC) and they were not collected on a regular basis every year. For this reason, some river types had more representation in one period than in the other (e.g. R-T32 type overrepresented in the 1st period and R-T22 in the 2nd; Table 1).
Water and diatom samples were collected from stones by brushing the top surface of each with a toothbrush, following the European standards (AFNOR, 2003;CEN, 2004;Kelly et al., 1998) and preserved with formaldehyde solution (37%) immediately after the collection. Then, an aliquot of 1 to 3 ml from each sample was treated with a 65% solution of nitric acid (HNO 3 ) and potassium dichromate (K 2 Cr 2 O 7 ) at room temperature for 24-48 h to remove the organic content. The acid residues were removed by centrifugation (1500 rpm), followed by rinsing with distilled water. After the oxidation, permanent slides were prepared with NaphraxÒ. Diatoms were studied under the light microscope (Olympus BX40), where a minimum of 400 diatom valves were counted on each slide and were identified to the lowest possible taxonomic level. The identification was based on Krammer & Lange-Bertalot (1986, 1988, 1991a, Krammer (1997aKrammer ( , b, 2002, Lange-Bertalot (1993;2001), Prygiel & Coste (2000), Trobajo et al. (2013), Wetzel et al. (2015) and Lange-Bertalot et al. (2017). Names were updated according to Algaebase (Guiry & Guiry, 2021).

Data analysis
Diatoms were collected in parallel to environmental variables, except for the years 2003 and 2020, for which no environmental data are available. An environmental variable sample matrix was produced and grouped by periods: the 2004-2008 samples corresponding to the 1st period and the 2016-2019 samples to the 2nd. The environmental matrix consisted of five variables: electrolytic conductivity (lS cm -1 ; n = 27), pH (n = 25), oxygen saturation (%; n = 14), nitrogen nitrate (mg N_NO 3 -L -1 ; n = 28) and phosphorus phosphate (mg P_PO 4 3-L -1 ; n = 24). Although data on water temperature were not available, it is known that the average values in this area present temperatures between 2.5 and 11 8C (Meléndez, et al., 2019). Electrolytic conductivity and oxygen saturation were determined in situ using a Hanna Instruments HI9829 multiparameter probe. N_NO 3 and P_PO 4 3were measured with an Auto-Analyzer 3 (Bran ? Luebbe, Germany). Normality was checked with the Shapiro test and student's t test was used to compare both periods. These analyses were performed using SPSS v.24 statistical software. The level of significance for the test was set as P \ 0.05. The values of the physical and chemical variables measured per site were compared with each variable boundary established for the classification of the ecological status classes of the corresponding river type defined in the Spanish normative (RD 817/2015). The most correlated environmental variables (Spearman's Rank correlation [ 0.6) were eliminated, and variables with several representative data for both periods were selected and fourth root transformed. Diatom abundances were log (x ? 1) transformed before Bray-Curtis similarity matrices were calculated. Then, a Distance-based linear model (DistLM; Legendre and Anderson, 1999;McArdle and Anderson, 2001) in combination with distance-based redundancy analysis (dbRDA) (Anderson et al., 2008) were used to identify the set of environmental variables that best explained diatoms distribution amongst samples. The best-fit model, based on Akaike's Corrected Information Criterion (AICc; Akaike, 1987) was selected and visualised using dbRDA ordination.
A permutation analysis of variance (PERMA-NOVA; McCune et al., 2002) was used on Bray-Curtis similarity matrix with 9999 permutations to test whether diatom assemblages varied across river basins, periods, river types and ecological status classes calculated with NORTIdiat and IPS. Nonmetric multidimensional scaling (nMDS) was used to summarise the variability of the two periods in terms of diatom community structure based on Bray-Curtis dissimilarity distance. All multivariate analyses were conducted in PRIMER7 v.7 software (Clarke & Gorley, 2015) with the PERMANOVA ? add-on package (Anderson et al., 2008). Student's t test was applied to establish the statistical significance of the difference between both periods of the abundances of the species belonging to the reference community. The level of significance for the test was set as P \ 0.05.
Finally, a SIMPER analysis (SIMilarity PERcentage; cut-off 90%) was performed using the PRIMER v.7 software to assess the degree of similarity between the diatom assemblages amongst periods, river types and ecological status calculated with NORTIdiat and with the IPS. The most common diatom indices in Europe (e.g. IBD, CEE, IDAP) together with the Specific Polluosensitivity Index (IPS, Cemagref, 1982), the Shannon Diversity Index, the corresponding evenness and the Trophic Diatom Index (TDI, Kelly & Whitton, 1995) were calculated with the Omnidia software v 5.3 (Lecointe et al., 2003). The scale of most diatom index values is from 0 to 20, except for TDI that is scaled from 1 to 100. The NORTIdiat index was calculated by applying the Bray-Curtis similarity to the matrix of the relative abundances of log-transformed samples based on the diatom reference community of river Type 1 (calcareous mountain rivers) (Pardo et al., 2018). The IPS and the NORTIdiat were used to calculate the ecological status class of each test sample. IPS values above 17 suggest a high ecological status, whilst IPS between 15 and 17 implies a good ecological status (Noga et al., 2014). NORTIdiat values above 0.930 suggest a high ecological status and values between 0.930 and 0.700 imply a good ecological status (Pardo et al., 2014(Pardo et al., , 2018. Data were Box-Cox transformed when required to meet the requirements of normal distribution (Shapiro-Wilk's test). Student's t test was performed to verify significant differences in diatom indices between the two periods studied and the two river basins (Sella and Deva-Cares rivers). The level of significance for the test was set as P \ 0.05.

Physical and chemical variables
The average electrolytic conductivity (EC) values decreased over time, from 262.0 lS cm -1 in the 2004-2008 period to 222.9 lS cm -1 during the 2016-2019 period. The pH values significantly (P = 0.000) increased with time from 7.7 to 8.4. The oxygen saturation was slightly higher during the 2nd period, with averages from 91.7 to 97.8% (Fig. 2).
Since we do not have dissolved oxygen data for all samples, it could not be included in the analyses, but the values range was between 8.44 and 11.05 mg l -1 .
Nitrogen nitrate (N_NO 3 -) and phosphorus phosphate contents (P_PO 4 3-) were significantly (P = 0.000 and P = 0.015, respectively) lower during the 2nd period Fig. 2 Boxplots of the four physical and chemical variables studied in both periods. The median values (central line), 25th and 75th percentile values (box) and the error bars are shown. Test with * significant at the 0.05 level and *** at the 0.001 level (Fig. 2). The DistLM analysis indicated that the bestfit model explaining diatom's distribution was composed of all the environmental variables (multiple R 2 = 0.28). The pH was positively correlated with the EC; thus, it was eliminated from the analysis. First two dbRDA axes explained 79.8% of fitted model variation (Fig. 3). The dbRDA1 explained 57.2% and separated samples from the first period with high levels of conductivity towards its positive part, from samples from the second period with low levels of conductivity on its negative end. The second axis (dbRDA2) explained 22.6% of the total variance and separate high nutrient content streams towards its positive part, from low nutrient content streams on its negative side (Fig. 3).

Analyses of diatom assemblages
The resulting diatom matrix consisted of 46 samples and a total of 132 diatom taxa (Supplementary material 1). PERMANOVA showed that there were no significant differences between diatom assemblages regardless of the river basin, Sella vs Deva-Cares (pseudo-F 1,42 = 1.65, P = 0.055). Although the nMDS ordination showed a slight grouping of the samples corresponding to the two studied periods (Fig. 4), the composition of diatom assemblages between the study periods was not significantly different either (pseudo-F 1,42 = 3.35, P = 0.164).
Significant differences in diatom assemblages amongst river types (according to the official Spanish typology) were not found (PERMANOVA, pseudo-F 4,31 = 1.42, P = 0.23). Moreover, nMDS ordination did not show a clear separation of samples amongst river types, instead, a dispersed distribution was observed in river type's samples (Fig. 6A, B). The SIMPER analysis indicated that the within group average similarity for each type ranged between 43.9 and 49.0% and the average dissimilarities amongst types ranged between 52.7 and 61.6%. The highest dissimilarities were found when comparing the R-T29 type with the R-T26 and R-T21 types. In both cases, the species that contributed the most to this dissimilarity were Achnanthidium atomoides and Cocconeis lineata.
The result of the diatom indices (IPS or NORTIdiat) indicate that all samples belong to good or high ecological status. PERMANOVA showed that there were no significant differences between the diatom assemblages regardless of the ecological status calculated with the NORTIdiat (PERMANOVA, pseudo-F 1,39 = 1.52, P = 0.08) or with the IPS (PERMANOVA, pseudo-F 1,39 = 1.48, P = 0.09). Thus, no segregation of samples was observed (Fig. 6A, B). It should be noted that many of the samples classified as in good ecological status by NORTIdiat (Fig. 6A) were classified in high ecological status by IPS (Fig. 6B). The SIMPER analysis showed that the dissimilarity between ''good'' and ''high'' ecological status was 56.9% for NORTIdiat and 58.1% for IPS.

Discussion
Our comparative analyses of diatom assemblages between two periods spanning 17 years in nonsignificantly impaired streams demonstrate a general persistence of diatom composition in Picos de Europa National Park. The characteristic species of the rivertype community and the high/good ecological status were maintained over time and the observed variations in their proportions should be attributed to the general maintenance or even improvement of the stream's environmental conditions. This study demonstrates that analysis of data collected during different periods can be successfully applied to support conservation and assess management strategies under changing European climates.

Temporal variation in water chemistry
Stream water composition resulted in slight, although important, variations in conductivity and nutrient contents. Electrolytic conductivity records are typical of non-impaired calcareous headwaters in North Spain during low flow periods in summer (Martínez et al., 2015;Pardo et al., 2018). The pH levels increased significantly, although they always remained within the range of calcareous streams (Corman et al., 2016). The results of the dbRDA revealed that period and EC were found to have the greatest explanatory power for diatom species distribution. However, no significant differences were found amongst diatom assemblages between periods. EC has already been found as the most influential environmental variable for diatom species (Pajunen et al., 2020). This study area presented a lower nutrient content than other calcareous geographical regions (Poikane et al., 2019). In areas with higher nutrient content, nutrients were the most important variables to explain the variation in the diatom species composition (González-Paz et al., 2020). Even though evidencing small decreases with time that should be attributed to a general improvement of water quality due to implemented good practices since 2012 within the natural park aiming to reduce the impact from agriculture and livestock activities. Moreover, these mountain streams have well-oxygenated waters. The observed significantly decrease in nutrients between periods clearly points to diffuse reduction in nutrient inputs. Higher N_NO 3 and P_PO 4 3in streams in the 1st period may be due to the agricultural expansion and aggregation detected between 2006 and 2012 in the Cantabrian Mountains (García-Llamas et al., 2019), as agricultural and breeding activities are the main source for N_NO 3 contents (Carayon et al., 2019;Gaglioti et al., 2019). However, the reduction in agricultural activity in the Cantabrian Mountains has caused land abandonment since 2015 (Perpiña-Castillo et al., 2020) seeming to coincide in time with lower N_NO 3 contents during the 2nd period of study. Our information on land-use indicated that all the studied sites have between 85 and 100% natural and semi-natural areas and between 0 and 15% agricultural areas. Recent studies in this area suggest that this trend of reduction of agricultural area will continue in the coming decades, reducing by 14.5% by the year 2049 with respect to the year 1989 (Hernández-Romero, 2020). Furthermore, the area most affected by livestock in the surroundings of the study area (Cabrales municipality) began to apply best practices in 2012.

Temporal changes in diatom assemblages
Diatom assemblages were stable as they persisted over time because we did not detect significant changes in their composition between the two periods. The most abundant species in the two studied basins, Achnanthidium pyrenaicum, Achnanthidium minutissimum and Cocconeis euglypta, are characteristic of calcareous rivers (Beltrami et al., 2012;Leira et al., 2017). The most abundant species in this study area together with Gomphonema pumilum (Grunow) E.Reichardt & Lange-Bertalot and Navicula tripunctata (O.F.Müller) Bory constitute the most abundant species of the specific biological reference community for type 1 of Mountain calcareous rivers, according to the NORTIdiat classification (Pardo et al., 2018). Both Achnanthidium minutissimum and Achnanthidium pyrenaicum have been found as dominant species in other temperate high mountain protected areas (e.g. Peszek et al., 2021;Poland et al., 2016). Cocconeis lineata, Nitzschia fonticola and Navicula cryptotenella significantly reduced their abundance in the 2nd period in the samples classified as good ecological status, whilst Achnanthidium pyrenaicum increased it. Nevertheless, all of them were part of the dominant species throughout the study. Achnanthidium minutissimum has been considered a very cosmopolitan taxon (Taylor et al., 2007) being dominant in a wide range of calcareous (e.g. France; Rimet et al., 2004) and siliceous geologies (e.g. Italy; Bona et al., 2007). Despite this, it is considered indicator of good water quality (Leira et al., 2017) across Atlantic and Mediterranean basins (Almeida et al., 2014;Martín et al., 2010) and of high ecological status in the IPS index (Lecointe et al., 1993). Moreover, it has been found to be of ''reference'' assemblages (Cantonati et al., 2020;Kelly et al., 2008). Nonetheless, Achnanthidium pyrenaicum is more sensitive than Achnanthidium minutissimum, regarding organic contamination, metals and other pollutants, because of its more specific alkaliphilus character (Cantonati et al., 2014). It should be noted that Achnanthidium minutissimum abundance during the first sampling year (2003) may be biased because Achnanthidium atomoides was not described (Monnier et al., 2004) and until then the latter species abundances were attributed to the very similar Achnanthidium minutissimum. However, Achnanthidium atomoides was not very abundant during the first period. A. atomoides has been found together with Cocconeis euglypta, both species related with high phosphate values (Beltrami et al., 2012). This preference could explain the observed reduction in the abundance of Cocconeis euglypta coinciding with P_PO 4 3decrease during the 2nd period. This species, together with Achnanthidium minutissimum, has been found amongst the most abundance species in the Cyprus stream reference community (Cantonati et al., 2020), in Luxembourg headwaters  and in other small mountain streams from the Iberian Peninsula (Gomà et al., 2005;Novais et al., 2020). Nitzschia fonticola, one of the most abundant species during the 1st period, that tends to positively respond to nutrient enrichment (Cochero et al., 2015), decreased with time suggesting a less favourable environmental conditions for the species.
On the other hand, in this study we only found the invasive species Didymosphenia geminata (Lyngbye) M. Schmidt (RD 630/2013) in one sample from one of the streams outside the park (SE012 in 2018) where its abundance represented 1.4% of the total. This invasive species is native to northern parts of Asia, Europe and North America (Anderson et al., 2020). The species is distributed in several rivers in Northern Spain since 2011 (Ladrera et al., 2018), but it was absent from the studied rivers during the 1st period studied in the Cantabrian Mountains. The species can be a threat for the biotic diversity of aquatic ecosystems and has to be controlled in its dispersion. The persistence of diatom assemblages between periods confirms our predictions that headwaters located in well-protected natural areas are appropriate systems to be considered reference ecosystems with type specific biological communities.

Differences amongst river types
Our results show that studied rivers present a great spatial similarity in diatom composition confirming the existence of a single diatom assemblage, thus corresponding to a single river type, the mountainous calcareous river (Pardo et al., 2018). PERMANOVA analyses confirmed the absence of statistically significant differences amongst the diatom assemblages from the different types. In conclusion, the NORTIdiat typology seems to be a more precise method to define river types in the area, as the other Spanish types did not show differences in their species composition, confirming the existence of a homogeneous diatom assemblage in the area. These results reinforce the prediction that river types based on biological reference conditions (i.e. invertebrates or diatom assemblages) are more reliable in Northern Spain than river types derived solely with physiographic descriptors. Therefore, determining the biological communities inhabiting under reference conditions is of great importance to provide more robust ecological status assessments in the future (Stubbington et al., 2017).

Temporal variation in structural metrics
This study was conducted in mountain streams in Northern Spain belonging to a protected area of high ecological value (Barquín et al., 2015) and this was reflected in high values for all community structural metrics that indicate a high/good ecological status, as for other protected mountain areas in Europe (Falasco & Bona, 2011;Falasco et al., 2012). On the other hand, a significant reduction in diversity and evenness was observed. The decrease in species richness and diversity was attributed to a 30% species disappearance that corresponded to species of the genus Nitzschia during the 2nd period. Most of the valves belonging to this genus were found in samples classified as good ecological status and a few were found under high ecological status class. The genus Nitzschia is considered tolerant to pollution (Hill et al., 2001), thus its disappearance should be related to a significant improvement in diatomological indices and thus in water quality. The IPS high values throughout both periods correspond to high ecological status (RD 817/2015) probably related to the reduction of nutrients that reduced the occurrence of Nitzschia species. Slight increases in IDAP index may be related to the decrease of NO 3 -, since it assesses micro-pollutants and organic contamination based on ammonium and nitrogen (Prygiel et al., 1996). TDI values were high compared to other protected areas (Jüttner et al., 2003;Tang et al., 2006). As it has been found that TDI values are positively correlated with phosphate contents (Jüttner et al., 2003), the decrease in this index during the 2nd period could be related to the decrease in P_PO 4 3contents. The NORTIdiat index showed slight non-significant reduction over time. This should be attributed to the natural variation in the diatom community at non-impaired sites, as the method uses the whole diatom assemblage to assess changes, being less sensitive to few species additions or absences.

Conclusion
Although, climate change has caused a rise in temperatures and fewer rainy days throughout the planet, affecting different ecosystems and reducing diversity, protected areas can act as natural reserves recovering degraded habitats and safeguarding intact ecosystems. The diatom assemblages found in these mountain streams persisted over time (2003-2008; 2016-2020) and maintained the high and good ecological status. Best agricultural practices seem to have reduced N_NO 3 and P_PO 4 3content along the years, further improving water quality. Even though there was a significant reduction in diversity and evenness over time, when comparing the diatom assemblages during two periods, it was observed that dominant species for both periods corresponded to the previously established reference diatom community for this area. This should be contemplated as an indicator of biodiversity maintenance, reflecting that those natural-protected areas have proven to be good biodiversity reservoirs and that headwaters maintain WFD requirements to be considered reference sites. Thus, further studies should address whether the relationship between changes in diatom assemblages and climate change is species-specific, rather than resulting in total assemblage change, as our results point out.