Using incomplete floristic monitoring data from habitat mapping programmes to detect species trends

The loss of biodiversity has raised serious concerns about the entailing losses of ecosystem services. Here, we explore the potential of repeated habitat mapping data to identify floristic changes over time. Using one German federal state as a case study, we assessed floristic changes between the 1980s and 2010s. These habitat data have great potential for analysis because of their high spatial coverage while also posing methodological challenges such as incomplete observation data. We developed a modelling approach that accounts for incomplete observations and explored the ability to detect temporal trends.

Over the last three decades, the European Union reporting obligations (EAFRD regulation, Natura 2000) have resulted in nation-wide assessments, which, however, differ between member states. The report from 2007 to 2012 on the Habitats Directive revealed that only 16% of habitats and 23% of species assessments were in "favourable" condition, while 77% of habitats and 60% of species assessments were in an "unfavourable" state (European Environment Agency, 2015). In Germany, the task of nature conservation and environmental planning is largely independently performed by the 16 federal states, while the German Federal Agency of Nature Conservation just aggregates the data (Hünig & Benzler, 2017;Sachteleben & Behrens, 2010). In addition to different monitoring schemes among federal states, trend analyses are hampered by changes over time in field survey assessment methods, resulting in heterogeneous data sets. This situation of heterogeneous biodiversity records also applies to many other regions worldwide, which has resulted in considerable uncertainties in the assessment of global biodiversity trends (Díaz et al., 2019).
To derive biodiversity trends in Germany, we explored the possibilities and limitations of using these heterogeneous data sets from different points in time. One of the available data sources for biodiversity change assessments is inventory data from a class of habitat surveys called "biotope mapping" (in German "Biotopkartierung"). In the following, we use the term "habitat" as synonymous to "biotope." Habitat mapping projects usually do not cover the entire landscape and focus on priority habitats that are legally protected according to the laws of the respective federal state.
Habitat mapping has been carried out in all federal states in Germany, starting in the 1970s (Kaiser, Schlumprecht, Finck, & Riecken, 2013;Sukopp, Trautmann, & Schaller, 1979). However, the main purpose of habitat mapping has not been the monitoring of biodiversity but the development of inventories of all legally protected habitats for conservation planning and protection. Moreover, only those habitats are mapped in the field that are legally protected according to the laws of the respective federal state. There are some states with repeated mapping campaigns, reaching back to the 1970s. Habitat mapping is usually geographically explicit, resulting in polygons with habitat type assignments. In many cases, species lists per polygon have been recorded in the field, but none of the federal states has combined habitat mapping with a full floristic or faunistic inventory.
Thus, assessing floristic changes based on repeated habitat mapping has so far not been considered feasible.
One such survey with extensive, if incomplete, species lists, with a focus on the habitat-typical species, is the habitat mapping scheme of the Federal State of Schleswig-Holstein. Starting as early as 1978, the first comprehensive mapping campaign across the whole state was completed by 1992 (SH1). A second (SH2) and third surveys (SH3) were carried out to map habitats protected by nature conservation law ( §30 BNatSchG of Germany and §21 LNatSchG of Schleswig-Holstein) and habitats with High Nature Value (HNV) farmland, respectively. A second comprehensive survey was started in 2014 (SH4) and has recently been completed.
For the present analysis, we employed data of the first and last habitat survey in Schleswig-Holstein (SH1 and SH4) to develop a new methodological approach for trend analysis with incomplete data. We matched the species lists of the habitat types in the two different surveys and developed tools compensating for varying degrees of completeness in terms of species recorded. This was needed to assess whether species have increased or decreased in frequency. Given the eutrophication at the landscape level, we expected an increase in nitrogen indicator species, which has been reported for Germany (Ewald et al., 2013;Immoor, Zacharias, Müller, & Diekmann, 2017;Wittig, Waldmann, & Diekmann, 2017) and other parts of Europe (Bernhardt-Römermann et al., 2015;Carey et al., 2008;Duprè et al., 2010;Pannek, Duprè, Gowing, Stevens, & Diekmann, 2015;Payne et al., 2013). In addition, we asked whether species losses are linked to certain species traits. We were particularly interested in the change of plant species that provide resources for pollinators, reflected in traits such as flower colour and nectar supply (Pakeman et al., 2017), which are of obvious relevance as potential drivers of change for pollinating insects (Biesmeijer et al., 2006). Thus, we assessed how past floristic changes might have impacted on a key ecosystem function.
Solving the methodological obstacles of such an analysis would allow using widely available yet currently untapped biodiversity information that is also available in other federal states in Germany.
This would have implications for many other countries, for example those under EU reporting needs for Habitats Directive, as most countries also do not have standardized repeated sampling for most taxa. We are, however, fully aware that any methodological mismatch between two surveys may be misinterpreted as unwarranted species turnover. Thus, we minimized the risk of a false detection error of change by taking most conservative evaluation strategy decisions and used modelled probabilities of species occurrences rather than raw data. For trend analyses with incomplete data, we therefore developed methodological approaches that should be transferable also to other habitat mapping data sets. Finally, we clearly point out the challenges and limitations in the discussion.

| Data preparation
The habitat mapping surveys in Schleswig-Holstein were confined to areas of interest for nature conservation. We used data from the habitat mapping from 1978 to 1992 (SH1, refers to Schleswig-Holstein, first habitat mapping) and from 2014 onwards (SH4; SH2 and SH3 focused only on selected habitats and were not part of this study). In the following, we refer to SH1 and SH4 as first and second surveys, respectively. SH1 and SH4 were conducted using different habitat identification keys (Landesamt für Naturschutz und Landschaftspflege Schleswig-Holstein, 1991; Landesamt für Landwirtschaft & Umwelt und ländliche Räume des Landes Schleswig-Holstein, 2017). As the last survey had not been concluded at the time of analysis, we only used the intersecting areas of both surveys, thus disregarding habitats in the first survey that had not yet been resurveyed.
In general, habitat types were defined more broadly in the first survey compared with the second survey. This resulted in a much higher number of polygons in the second survey, that is 298,185 and 611,572 polygons in the first and second surveys, respectively.
In the following, we only kept those polygons (13,416 and 40,280, respectively) that were linked to species lists. We unified the nomenclature of all species using a German standard list (Jansen & Dengler, 2008;GermanSL v1.4) and the vegdata R package (v. 0.9.2).
Spelling errors were corrected and names without a match in the standard list (n = 496) were removed. Subspecies were combined at the taxonomic level of species, and where applicable, species at the taxonomic level of species aggregates, using the aggregate definitions by Horst (2010). The surveys together contained 1,547 species, with 1,288 and 1,317 species recorded in the first and second surveys, respectively. Using qgis 2.14, we intersected the polygons of both surveys (for an example, see Figure 1). This resulted in 5,503 polygons (51%) of the first survey intersecting with 20,559 polygons from the second survey (41%). This is due to the fact that some habitat types were recorded only in the second but not the first survey (e.g. mesophytic grasslands). On average, polygons were much larger in the first than in the second survey, with mean areas of 1.57 ha and 1.24 ha, respectively, caused by the broader habitat definitions, yet polygon size distribution showed a large overlap between both surveys ( Figure S1). In total, we retained 60 habitat types, defined according to the first survey's mapping key, which comprised a broad range of habitats such as tall sedge swamps, bogs and dunes. This is a conservative decision as this procedure ignores losses of whole habitats. In the following, we based all evaluations on the polygons of the first survey that intersected with the second survey.

| Accounting for sampling intensity
To derive a reliable estimate of species richness per polygon, we had to account for the higher observation effort in the second survey. first survey and joining their species lists resulted in a higher total number of observations in the second survey. Assuming that species detection probabilities depended on sampling intensity, species richness of these matched polygons would be expected to be higher in the second survey. To account for the higher sampling intensity, we applied rarefaction to the second survey and calculated standardized effect sizes (SESs) according to Gotelli and McCabe (2002) (1).
where SR obs_surv1 and SR exp_surv2 are the mean observed and expected species richness per polygon and σ is standard deviation of expected species richness. The observed value of species richness in the first survey was taken as SR obs_surv1 . To obtain SR exp_surv2 , we randomly resampled from all polygons in the second survey 100 times the same number of species x polygon observations that were encountered in the first survey (126,375). We calculated SESs as the observed value of species richness in the first survey minus the mean richness in the 100 random runs, divided by the standard deviation of mean richness of the 100 random runs (Equation 1). Using the threshold of p < .001, SESs were considered significant (i.e. a significant change in richness between the first and second surveys) when they were 3.29 standard deviations smaller or larger than the observed values. As for all further tests below, we employed this low error threshold of p < .001 to reduce the type I error rate and, thus, only make conservative estimates of a potential change.

| Quantifying change making use of communitylevel information
The central problem of using species lists from habitat mapping surveys is the incompleteness of the lists. This can be illustrated by analysing log 10 richness-log 10 area relationships ( Figure S2), where most polygons show a species richness much below the average numbers expected in the typical Northern German landscape. Thus, we cannot exclude the possibility that species may appear to have become rarer or more frequent between the two surveys simply because of shifts in attention away or towards these species. The different amount of attention given to a particular species may be the result of different habitat mapping keys, where certain species were part of the definition of a habitat or not. In addition, there is also an unknown observer bias, with people varying in species identification skills, previous experience and amounts of time when taking records in a certain area. Using rarefaction, as we did for assessing changes in species richness (see chapter 2.2) is not a solution here, as it would aggravate the problem of incomplete observations and discard important information on the species composition of habitats. Thus, in order to account for incompleteness we took a different approach and calculated the occurrence probability of every species in every polygon, using Beals' index of sociological favourability (Beals, 1984) according to formula (2).
The probability p pi for species i to occur in a polygon p is calculated from joint occurrences M ij with all species j of the total number of species in that polygon N p , divided by the number of plots M j in which the species j is present. As a result, we obtained occurrence probabilities for every species in a polygon, including those that were not observed in a polygon (see Figure 2a for an example) or even not observed at all in one of the two surveys. We based the co-occurrence matrix M ij on all polygons of both surveys (53,696 polygons × 1,547 species), thus assuming that the species co-occurrence matrix is static in time. Thus, we excluded the possibility that a species' occurrence probability in a polygon with the same species composition could change because the species occurred more or less often with each other in the two surveys. Using the same co-occurrence matrix M ij across both surveys is a conservative decision as a non-static co-occurrence matrix would have allowed larger changes in species occurrence probabilities. Figure 2b illustrates the robustness of co-occurrence probabilities towards overlooked or otherwise underrecorded species. Thus, a change in occurrence probability of a species is mainly brought about by a change in species composition of that polygon, that is by changes in the presence of all the species in that polygon. The occurrence probability increases (or decreases) if the species composition of the polygon becomes more (or less) similar to polygons where the species is usually found. This means that Beals' index does not only capture biotic interactions (Breitschwerdt, Jandt, & Bruelheide, 2018) but also overall habitat quality. Additions or disappearances of single species also change the occurrence probabilities of the other species in a community ( Figure 2b).
We applied the formula of Beals' index to the full data set. We then aggregated the occurrence probabilities for each species across all polygons in the second survey, using the highest value if more than one polygon of the second survey intersected with a polygon of the first survey. This is also a conservative evaluation decision, as this procedure generally tends to result in higher occurrence probabilities of species in the second survey. As a result, we obtained occurrence probabilities for every of the total of 1,547 species in 5,503 polygons for each of the two surveys, from which we only retained the 1,377 species that had actually been observed in either of the two surveys. To determine whether a species declined or increased in occurrence probability, we applied a two-sided paired t test (Table S1). We adjusted the significance levels for multiple testing using Holm adjustment.
In addition, we analysed the change of observed frequencies of species in polygons based on raw data. Although this approach is only suitable for unbiased and complete data, we report these results to allow comparisons with other studies. For every species, we calculated the number of polygons of the first survey in which the species was observed exclusively in the first, second or both surveys. In each polygon, a species could disappear (−1), remain (1). To test whether these changes in frequencies were significant across all polygons, we applied a sign test, correcting for multiple testing by applying Holm adjustment (Table S2). To compare the approach based on occurrence probability with that using raw data, we compared the ranking of species with respect to the change in occurrence probability with the change in actually observed frequency using Spearman's rank correlation coefficient r S .

| Simulation test for assessing uncertainty of detecting change
To test the robustness of filling gaps by applying Beals' index, we created simulated data sets. To make this test as reliable as possible, we produced the simulated data with the actually observed records of the second habitat survey in Schleswig-Holstein, using the species information of all polygons from that survey that intersected with the first one (133,797 species × polygon combinations). We produced a second data set in which average species richness per polygon was kept constant but some species increased while others decreased. This was achieved by producing an identical copy of the first survey data set and then deleting the same number of occurrences in both sets, deleting occurrences of one set of species in one, and another set of species occurrences in the other data set. To create trends in species frequencies of varying strengths, the two sets of occurrences to be deleted were created by ranking the species by frequency, and then deleting 30 occurrences in every second species, or all occurrences if the species had a frequency ≤30 occurrences. We applied this to the species with uneven and even ranks in frequency in the first and second simulated surveys, respectively.
In this way, we produced two simulated surveys in which the original co-occurrence structure was kept and the observation effort was identical (133,797 and 133,791 species × polygon combinations, respectively), but introduced a clear trend, with the same number of species x polygon combinations showing an increase and decrease (8,281 and 8,275). In this way, we produced a reference pair of two surveys with known changes ("true change").
To test how the lower sampling intensity in the real data of our first survey may affect the detectability of changes, we introduced decreasing levels of observation completeness in the first survey of the simulated reference data set by randomly deleting 10%, 20%, … 90% of all records of the first survey, irrespective of species or polygon. Thus, we assumed that detection probability was not related to species' trends. We then compared the trends in species frequencies in each of the nine resulting data sets of decreasing degree of completeness with the "true" reference. We did this both using raw frequencies of species and occurrence probabilities. The latter were obtained from applying Beals index, calculating occurrence probabilities across both the first survey (thinned out to different degrees) and the second survey (not thinned out), following the assumption that the species co-occurrence matrix is static in time (see above). We calculated the mean change per species in all scenarios and compared the changes with the "true" reference.

F I G U R E 2
Robustness of Beals' occurrence probabilities: example for a randomly selected polygon, where the presence of species is indicated by light grey bars and the probability of occurrence based on Beals' index by dark grey bar. The occurrence probabilities were calculated using the full co-occurrence matrix M ij of 1,547 × 1,547 species. Here, we show only the values for the 15 species with highest occurrence probability. (a) Occurrence probabilities based on seven species. There are eight species that did not occur in the plot but had a higher occurrence probability in this plot than the species with the lowest value (Epilobium hirsutum). (b) As (a), but one species removed (Juncus effusus), indicated by the arrow, thus basing the occurrence probabilities on six species only. As a result, occurrence probabilities change slightly, while the missing species Juncus effusus still ranks among the four species with highest occurrence probability. Note that changes are much smaller when more species are present in the plot. This example is based on seven species only, while the average number of species recorded was much higher (see Section 3)

| Species characteristics
We tested for a systematic bias for species that were explicitly listed in the habitat definitions in the second survey, to which surveyors had to pay attention because they were listed as indicators for environmental conditions or those that were commonly known to be characteristic of Natura 2000 habitats. We compared both the change in the probability of occurrence and raw frequency of species that were either listed or not listed in one of these three lists. As also more attention might have been paid to threatened species, we also compared the amount of change with respect to the species' Red List status (Landesamt für Natur und Umwelt des Landes Schleswig-Holstein, 2006), using the three categories of conservation concern: (1) critically endangered, (2) endangered and (3) vulnerable.

| Polygon-weighted characteristics
To compare changes in the species composition of polygons, we calculated polygon-weighted (=community-weighted) means or medians.
For all vascular plant species, we retrieved Ellenberg indicator values (IV) for nitrogen (N), moisture (F), soil pH (R), light (L) and temperature (T) from the Biolflor database (Ellenberg et al., 1992;Klotz, Kühn, & Durka, 2002). For the first and second surveys, we calculated the weighted median for IV for every polygon, weighting the IV of all species in a polygon by the occurrence probability in that polygon. For a comparison with the raw data, we also calculated unweighted medians across all species that were actually observed. Finally, we analysed the relationship of the observed change with respect to species characteristics. As important traits for pollinators, we selected flower colour and nectar supply, retrieving these traits from Biolflor (Klotz et al., 2002).
We converted the colour names (blue, brown, yellow, green, orange, purple, red, pink, violet and white) into RGB colour code values, using the col2rgb command in R. We then calculated weighted means for each of the ten colours above and the R, G and B values as well as for the information of whether the species supplies nectar to pollinators, either based on the probability of occurrence or on unweighted means based on frequency for every polygon in the first and second surveys.

| Temporal trends and their reliability
Mean species richness, based on the raw data of the polygons in the first survey and the polygons of the second survey matched to those of the first, was 23.0 and 32.2, respectively (significantly different according to a paired t test, p < .001). However, when applying rarefaction and using the same sampling intensity as in the first survey, mean richness of the second survey was only 18.5 species, that is lower by 4.5 species than in the first survey (SES = 124.2, p < .001).
Our simulations showed that the ability to detect trends decreased with decreasing degree of completeness ( Figure 3). When based on raw frequencies, the correlation coefficient between the species changes in the "true" (albeit simulated) reference data set and thinned-out data sets decreased by 0.081 per every 10% of data loss in one of the surveys (Figure 3a). In other words, the raw data reflected the true changes less well with decreasing completeness of observations. In contrast, the correlation coefficients were much larger when the analyses used occurrence probabilities based on Beals' index. Moreover, "true" trends were increasingly underestimated with decreasing completeness when based on raw frequencies, while the direction of changes was more or less equally well detected when based on occurrence probabilities ( Figure 3b).
As supported by the simulation results, we present the observed results of the real data for temporal trends based on the probability of occurrence (Beals' index), and for all results based on observed frequencies (i.e. raw data) refer to the Supporting Information. Out of the total of 1,377 species, 499 and 878 declined and increased in occurrence probability, respectively. A total of 502 of these changes were significant at p < .05 after adjusting for multiple testing, with 172 and 330 species showing either a significant decline or a increase in occurrence probability (Table S1). The absolute change in occurrence probability between both surveys was lower for declining (−0.00153) than increasing species (0.00173). Similarly, in the analysis of raw data the mean change in observed frequencies was −0.567 compared to 0.573, with 103 and 157 species losing or gaining in observed frequency, respectively (Table S2). The two approaches of using probabilities of occurrence or observed frequencies resulted only in a poor correlation of the species' trends (Spearman r S = .54, p < .001 across all n = 1,377 species, and r S = .45, p < .001 across all species with a significant change in occurrence probability). For example, some nitrophilous species such as Urtica dioica and Galeopsis tetrahit showed a decrease in occurrence probability but an increase in observed frequency. Other cases, such as the moisture indicator species Peucedanum palustre and Cirsium oleraceum decreased both in occurrence probability and observed frequency.

| Species characteristics
We did not encounter a systematic bias in favour of species that were either explicitly listed in the second survey in habitat defini- There were strong differences between habitat types with respect to changes in species occurrence probability and observed frequencies (Table S3, Figure S3 and S4). The strongest mean decline of species' occurrence probabilities for habitats with high number of species observations (>10,000 species x polygons observations) was found in hedgerows ("Knicks"), headwaters, tall sedge swamps, dry excavation areas, fens, wet grassland and mesophilous grassland. In contrast, lagoons, salt marshes of the North Sea, coastal dunes and bogs showed the strongest increases in species occurrence probabilities. The changes per habitat type based on occurrence probabilities and observed frequencies were strongly correlated (r S = .79, p < .001).
Among the significantly increasing species were many neophytes,  (Figure 4). The decline in RL category 3 species was significant both with respect to probability of occurrence ( Figure 4) and observed frequency ( Figure S5).  (Table 1). When these colours were analysed for their RGB components, there was a decrease in reddish and bluish colour values and an increase in green values, both with respect to occurrence probability ( Figure 6) and raw frequency ( Figure S7, all significant at p < .001

| Polygon-weighted species characteristics
according to paired t tests). This shift towards green colour values was mainly brought about by a proportional increase in the colours green and yellow, and a decrease in the colour pink. Simultaneously, plants with nectar supply decreased (Figures 7 and S8).

| D ISCUSS I ON
We were able to match the two habitat surveys in the Federal State

| Reliability of the analyses
While the habitat mapping of our study region was carried out by experts, there remains much uncertainty in the conclusions we can draw from these data due to differing methodologies, both in spatial allocation of recording sites and the incompleteness of species lists. The former could be addressed by very conservative comparisons, using only spatially overlapping units. The incomplete species lists, however, are a severe problem that can be alleviated by applying Beals' index. This was tested and confirmed by our simulation approach, which showed that more reliable trend estimates are obtained from occurrence probabilities than from raw species frequencies.
A key feature of using occurrence probabilities instead of raw frequencies is the correction for sampling intensity. While the raw frequency of a species may simply increase by adding species lists of several intersecting polygons in the second survey, a species' Beals' index in that polygon only increases if the overall species composition suggests that this species would be expected to grow there. However, the ability to detect a given trend diminishes with decreasing degree of completeness. With an observed mean species richness of 32.2 species per polygon, our observations probably range at the lower end of the range of incompleteness studied in our simulated experiment. However, the results suggest that there is certainly a limit of applying Beals' index to account for incompleteness with very gappy data.
Our finding revealed similar patterns of changes in probability of occurrence and changes in raw frequencies. This shows that the data were still sufficiently complete to detect trends also by raw frequencies. However, more trends were found in occurrence probabilities. There were almost twice as much species with significant changes in occurrence probabilities (502) as compared to observed raw frequencies (260). The conservative decisions in the evaluation procedure make us confident that the observed trends reflect true changes. Because of these conservative decisions, we are aware that we probably underestimate the severity of plant diversity changes.
Given that the habitat mapping was never intended to monitor biodiversity change, we cannot fully rule out a systematic bias in the recorded species. For example, we suspected that in the second survey specific species received more attention than in the first one, as a consequence of the more detailed mapping instructions that encouraged surveyors to record certain species mentioned in the habitat identification keys. This may explain why a common species such as Urtica dioica that was listed as characteristic for different types of nitrophytic habitats (species list of ruderal species; Landesamt für Landwirtschaft & Umwelt und ländliche Räume des Landes Schleswig-Holstein, 2017) increased in recorded raw frequency. This trend was reversed when being based on occurrence probabilities, indicating that the species may have been recorded more often in habitats that had not the typical species composition of nytrophyitic communities. Overall, our finding that most of the indicator species listed in the mapping instructions of the second survey did not increase but significantly decreased in occurrence probability also suggests that no systematic bias was introduced by mapping instructions. Another systematic bias may have been introduced in the joint co-occurrence matrix by the more narrowly defined habitat types in the second survey. As a consequence, the occurrence probability of some species may have changed simply because they are no longer associated with some species from the first survey. For example, if nitrophilous species were associated with species of more oligotrophic habitats in the first but not in the second survey, they may decrease in occurrence probability. This would explain some discrepancies between species that increased in observed frequency but decreased in occurrence probability.
Another problem in the data set is the difference in sampling intensity between both surveys, resulting from more F I G U R E 4 Change in probability of occurrence by status in the Red List of vascular plants of Schleswig-Holstein (Landesamt für Natur und Umwelt des Landes Schleswig-Holstein, 2006). 1, critically endangered (81 species); 2, endangered (101 species); 3, vulnerable (116 species). The change in probability of occurrence was significantly higher than 0 for category 1 (p = .013), not significantly different from 0 for category 2 (p = .446) and significantly lower than 0 for category 3 (p < .001). The plots show mean occurrence probability of all species as horizontally jittered points, the median as horizontal bar, the 95% frequentist confidence interval as box and smoothed density as coloured "bean" (Kampstra, 2008). For raw frequencies see Figure S5 differentiated habitat types and a larger number of polygons mapped in the second survey. One consequence of this bias was that there were many more species that increased in occurrence probability than declined, which is contradictory to steadily increasing Red Lists of endangered species in Central Europe

| Temporal trends
An unexpected result was the absence of a trend for highly threat-  (Schaap et al., 2018). One key issue here is that we can only compare two points in time, which precludes statistical analysis of the underlying drivers of biodiversity (Pakeman et al., 2017). However, there is plenty of evidence from plot-based studies that community-weighted mean values in nitrogen indicator values have increased in northern Germany (Immoor et al., 2017;Wittig et al., 2017).  , 1978-1992) and second (SH4, 2014 onwards) habitat mapping surveys. N nutrient supply, F soil moisture, R soil pH, L light and T temperature. To make the differences visible, the crosslines show means rather than the overall median. Differences between the two surveys were tested for a statistically significant differences by a paired Wilcoxon test. ***p < .001, **p < .01, n.s., not significant. Arrows indicate the direction of change. For median IV based on frequency see Figure S6 in mean nitrogen values. Similarly, the decrease in mean soil moisture value supports resurvey results of wet grasslands near Bremen, which corresponds to extensive drainage and decreasing groundwater tables (Immoor et al., 2017;Wittig et al., 2017).
In addition, the landscape has become greener in terms of flower colour with losses of blueish and reddish flower colours, which suggests a decrease in insect-pollinated plants and an increase of wind-pollinated species that do not need to attract pollinators and most often have greenish flowers, such as grasses (Poaceae) and sedges (Cyperaceae). Consistent with our results, grass species have been reported to proportionally increase in floodplain meadows in northern Germany  as well as in acidic grassland plots in Germany, The Netherlands and United Kingdom (Biesmeijer et al., 2006;Duprè et al., 2010) or Scottish grasslands (Pakeman et al., 2017). This interpretation is corroborated by a concomitant decrease of species providing nectar. A lower availability of resources for pollinators has also been reported from permanent plots of Scottish grasslands (Pakeman et al., 2017). To our knowledge, our study provides the first evidence that these trends are relevant at the landscape scale. There may thus be an important link to the observed insect decline in Germany (Hallmann et al., 2017), given that many insects are bottom-up controlled by resource supply of plants (Scherber et al., 2010). This should raise serious con-

| Conclusions
Our attempt to make use of disparate information on the state of biodiversity in one German federal state has revealed both the possibilities and challenges of using such heterogeneous data. As all federal states in Germany have carried out habitat mapping, and repeated TA B L E 1 Flower colours of species that occurred in a polygon in the first (SH1, 1978(SH1, -1992 and second (SH4, 2014 onwards) habitat mapping surveys Note: n species = the number of species with the respective flower colour; R, G and B = colour values used in the ternate diagram of Figures 6 and S7. Occurrence probabilities in the first and second surveys are mean values of the colour being present in a polygon, weighted by occurrence probability of the respective species. Observed frequencies are the unweighted mean values of the colour being present in a polygon. Differences between the two surveys were tested for statistically significant differences by a paired t test.
habitat mapping data are also available in many other European countries, we see a high potential of applying our approach to other projects and compare species trends across regions. Importantly, our method reveals relevant points for designing future habitat mapping and biodiversity monitoring programs, and highlights the advantages of analysing changes at the species level compared to the community level. However, the uncertainties in our data clearly show that such an analysis cannot replace a structured and well-designed monitoring program in Germany, which calls for a coordinated and standardized scheme across all German federal states. Furthermore, our results point to the importance of paying attention to species characteristic of oligotrophic habitats and wet soil conditions. Thus, indicators in biodiversity monitoring should not just comprise species sensitive to change but also those with important ecosystem functions. In particular, we suggest taking into account flower colour and nectar supply in the design of monitoring programs and the selection of target plant species as these may be crucially linked to changing abundances of insect pollinators and thus indicate changes in species networks and ecosystem service provision.

| DATA ACCE SSAB ILIT Y
The data of both surveys are available in an Oracle database of the State Agency of Agriculture, Environment and Rural Areas of Schleswig-Holstein (LLUR). The observed changes in occurrence probabilities and observed frequencies are listed in Supporting Information Tables S1 and S2.

ACK N OWLED G EM ENTS
We are grateful to field surveyors who mapped biotopes in the field since the 1980s and made this analysis possible. The present study is an outcome of the sMon project (Trend analysis of biodiversity data in Germany) of the German Centre for Integrative Biodiversity  , 1978-1992) and second (SH4, 2014 onwards) habitat mapping surveys. The ternate diagram shows the weighted mean values for red (R), green (G) and blue (B) in the first and second surveys in red and cyan, respectively. The centroids in dark red and dark cyan show the mean colour values across all polygons in the first and second surveys, respectively. For mean colours based on observed frequencies, see Figure S7. For the observed frequencies of the original colours, see Table 1 F I G U R E 7 Mean nectar supply values for all polygons of the first (SH1, 1978(SH1, -1992 and second habitat mapping surveys (SH4, 2014 onwards). The species' nectar supply (yes = 1, no = 0) was weighted by the species' probability of occurrence in the polygons. Differences between the two surveys were tested for statistically significant differences by a paired t test. ***p < .001. For mean nectar supply values based on observed frequencies, see Figure S8