Spatio-temporal Variability in Benthic Macroinvertebrate Communities in Headwater Streams in South Korea

Comprehensive research on the structural and functional variability of benthic macroinvertebrate communities within headwater streams is limited, despite the fact that the majority of streams within a watershed are headwater streams that form the primary link between terrestrial and aquatic ecosystems. Therefore, we investigated the structure and function of benthic macroinvertebrate communities in four headwater streams at two different spatial scales autumn) of the year. Community indices, functional feeding guilds and habit trait guilds varied significantly depending on the seasons rather than on sites in two-way ANOVA based on spatial (i.e., sampling sites) and seasonal effects in each headwater stream. Non-metric multidimensional scaling analyses showed the differences between communities according to the considered spatial and temporal scales. At the individual stream scale, the differences between samples followed seasonal variation more than spatial differences. Site differences became more important when performing an ordination within a single season (i.e., spring, summer, and autumn). Continued research and monitoring employing both multidisciplinary and multidimensional approaches are required to maintain macroinvertebrate diversity within headwater streams.


Introduction
Biodiversity has been declining at an increasing rate worldwide [1] as a result of anthropogenic habitat disruption.Although freshwater occupies less than 1% of the Earth's surface area, and rivers and streams represent only 0.006% of all freshwater resources [2], they exhibit high biodiversity, comprising approximately 10% of known species [3,4].
Headwater streams are extremely heterogeneous ecosystems with high spatial and temporal variation [5], comprising a significant proportion (i.e., more than three-quarters) of the total stream channel length within a watershed [6].Headwater streams are main sources of water, sediments, and organic materials that are transported downstream [7][8][9][10], and their small catchments couple terrestrial and aquatic ecosystems such as food web dynamics [11,12] including allochthonous input [13], inputs of terrestrial invertebrates [14], etc. (see Nakano et al. [15] for a detailed explanation).Furthermore, they are essential for sustaining the structure and function of watersheds [7,8,10,16].Headwater streams provide valuable habitats for unique and diverse communities of aquatic flora and fauna [16][17][18].Therefore, it has become increasingly clear that headwater streams are essential for maintaining biodiversity in both terrestrial and aquatic habitats [7,8,10,17,19,20].
Benthic macroinvertebrates perform central ecological roles in stream ecosystems [21], such as processing of detritus, participation in animal-microbial interactions and functioning as primary and secondary consumers through critical trophic interactions [22,23].Headwater streams are characterized by diverse microhabitats (i.e., refugia) that help protect macroinvertebrates from competition, predation and natural disturbances, and therefore support a rich regional biodiversity [20].Research on the environmental and biological parameters that determine the structure and function of macroinvertebrate community in headwater streams is essential for the basic understanding of the ecology, biodiversity, and conservation of these important ecosystems [24,25].
The composition of the macroinvertebrate community can be differentiated by various factors, including latitudinal gradients [26], stream segmentation and microhabitat [27,28].Heino et al. [26] suggested that local filters (e.g., water quality) in headwater streams were relatively weak whereas they showed the clear latitudinal gradients of macroinvertebrate community composition.Ligeiro et al. [27] found that the composition of macroinvertebrate community was differentiated according to stream segments and microhabitats in a tropical headwater catchment, and García-Roger et al. [28] reported that during the dry season, the species richness was decreased especially in the temporary headwater streams due to the reduction of available habitats.The diversity of different guilds (i.e., functional feeding guilds and habit trait guilds) in headwater streams is affected by pH, stream width, moss cover, stream particle size, nitrogen, and water color [19].Moreover, algae-scraping invertebrates represent longitudinal zonation patterns along the river systems whereas within riffles, algal abundance can determine the invertebrates in small-scales [29][30][31].The distributions of leaf-shredding invertebrates often reflect longitudinal and among-stream variability in riparian conditions [32,33] as well as riffle-scale patchiness of leaf detritus on stream bottoms [34,35].Chung et al. [36] reported that the variation in the trophic structure was affected by habitat characteristics in each channel reach, including channel morphology, proportion of habitat type, and benthic organic matter availability.However, there has been little research on aquatic biodiversity in headwater streams considering both seasonal and spatial differences.
Therefore, we examined the diversity of a benthic macroinvertebrate community in four different headwater streams at two different spatial scales (i.e., sampling sites >samples (riffles)) in three different seasons (i.e., spring, summer, autumn).We tested hypothesis that the composition of macroinvertebrate communities would be spatially and temporally heterogeneous at different spatial scales in headwater streams [37][38][39].We considered only headwater streams free of anthropogenic disturbance to exclude interaction effects between anthropogenic and natural factors on macroinvertebrate communities.

Study Area
We studied benthic macroinvertebrate communities at the headwater streams in four different regions of the northern part (Gwangreng: GR and Hongcheon: HC) and southern part (Wando: WD and Geumsan: GS) of South Korea (Figure 1 and Table 1).All streams were in forested areas, free of anthropogenic disturbance (Table 1).For instance, GR and WD are in the National Arboretum and people have rarely visited HC and GS due to the accessibility.Acer pseudosieboldianum, Quercus mongolica and Securinega suffruticosa were dominant trees in riparian areas of GR, Sambucus racemosa L. ssp.sieboldiana and Deutzia grandiflora Bunge var.baroniana were dominant in HC.Meanwhile, the riparian vegetation of GS was mainly composed of Pinus densiflora, Styrax obassia, and Phragmites japonica, and Eurya japonica, Camellia japonica and Quercus acuta were mainly observed in the riparian vegetation of WD.There were no houses or farms in the stream catchments of study areas.All sampling sites were in the first or second order streams based on a geographical map (scale: 1:50,000).There were clear gradients of climate (i.e., temperature and precipitation) according to the climate data from the Korea Meteorological Administration (KMA) [40].Annual precipitation in the study areas was higher in the southern study area (WD: 1532.7 mm and GS: 1512.8 mm) than in the northern study area (GR: 1450.5 mm and HC: 1405.4 mm).Due to the monsoon climate, more than 50% of the precipitation was concentrated in summer (especially, June or July to August); whereas other periods (mainly from October to March) were dry [41].Annual average temperature based on the data from 1980 to 2010 from KMA is the lowest in HC (10.8 ˝C) followed by KR (12.7 ˝C), KS (13.4 ˝C) and WD (14.3 ˝C).Monthly temperature range is the highest in HC from ´11.5 ˝C to 30.2 ˝C followed by GS (´5.8 ˝C-30.3˝C), GR (´5.9 ˝C-29.6 ˝C) and WD (´0.4 ˝C-29.2˝C).

Ecological Data
Benthic macroinvertebrates were collected with a Surber sampler (30 ˆ30 cm, 300 um mesh) to a depth of 10 cm at 12 sampling sites in four different streams (Figure 1).Sampling was conducted seasonally in spring, summer, and autumn in 2009 (GS), 2010 (GR), 2011 (HC), and 2014 (WD).Samples could not be collected in winter because the streams were frozen.In each stream, three riffle sites (e.g., GS1, GS2 and GS3 in GS stream) were selected at less than 0.5-km intervals between the adjacent sites.Within each riffle, three to five replicates were sampled on a longitudinal direction within 1-to 3-m distances between the adjacent sampling replicates (see [26,42]).Therefore, a total of 177 samples were collected (four streams ˆthree sites ˆthree-five replicates ˆthree seasons).In the laboratory, macroinvertebrates were sorted and preserved in 70% ethanol.All the individuals were identified mainly to the species level except Chironomidae under a stereo microscope (SMZ800N) at 400ˆbased on literature [43][44][45][46][47][48].

Ecological Data
Benthic macroinvertebrates were collected with a Surber sampler (30 × 30 cm, 300 um mesh) to a depth of 10 cm at 12 sampling sites in four different streams (Figure 1).Sampling was conducted seasonally in spring, summer, and autumn in 2009 (GS), 2010 (GR), 2011 (HC), and 2014 (WD).Samples could not be collected in winter because the streams were frozen.In each stream, three riffle sites (e.g., GS1, GS2 and GS3 in GS stream) were selected at less than 0.5-km intervals between the adjacent sites.Within each riffle, three to five replicates were sampled on a longitudinal direction within 1-to 3-m distances between the adjacent sampling replicates (see [26], [42]).Therefore, a total of 177 samples were collected (four streams × three sites × three-five replicates × three seasons).In the laboratory, macroinvertebrates were sorted and preserved in 70% ethanol.All the individuals were identified mainly to the species level except Chironomidae under a stereo microscope (SMZ800N) at 400× based on literature [43][44][45][46][47][48].All specimens were categorized into both functional feeding guilds (FFGs, predators: PR, scrapers: SC, collector-gatherers: CG, collector-filterers: CF, and shredders: SH) and habit trait guilds (HTG, clinger: CL, burrower: BU, swimmer: SW, sprawler: SP, and climber: CM) based on Merrit and Cummins [34], except Chironomidae, because of the difficulties in taxonomic classification.

Data Analysis
We conducted two steps of analyses to compare the differences between macroinvertebrate communities according to the spatial and temporal differences.First, variations of community indices (abundance, species richness, Shannon diversity index, Simpson diversity index, and Evenness) and proportions (%) of each class of FFGs and HTGs were analyzed using two-way analysis of variance (two-way ANOVA) to determine spatial and/or seasonal differences in each headwater stream.Second, we analyzed the abundance of macroinvertebrates using non-metric multidimensional scaling (NMS) and the Bray-Curtis distance to identify the relative differences between the sample units over multiple spatial scales and seasons.NMS is an indirect ordination analysis that compares the distribution of the macroinvertebrate community across all the sampling units without including any prior information about how the structure or taxa of macroinvertebrates could be altered or respond to environmental variables [50,51].NMS was applied to the datasets at two different spatial scales: (1) each individual stream (three sites each) and (2) each site.
Prior to NMS analyses and statistical tests, we transformed the abundance of each taxon that showed large variations using the natural logarithm.Before transformation, the number one was added to the variables to avoid the logarithm of zero [52].
Two-way ANOVA were conducted with the package stats in R software [53], and NMS analyses were conducted with PC-ORD version 5 [54].
The seasonal differences in community indices, FFGs and HTGs were mainly observed more frequently than the site differences except GR (Tables 3-5).For instance, their statistical differences (i.e., community indices, FFGs and HTGs) were relatively larger among sites in GR (9 in 15 cases).Only scrapers and shredders showed seasonal differences or spatial differences in all cases (i.e., sites, season and interaction between sites and season).In HC, species richness, Shannon diversity and scrapers showed seasonal differences.Only swimmers showed significant differences among sites.In GS and WD, the frequencies of seasonal differences were also higher (e.g., species richness, collector-gatherers, clingers, burrowers and swimmers in GS) than among sites (e.g., evenness, predators in GS).In the NMS ordination for each stream, the distribution of the sampling units reflected seasonality rather than the differences among sites (first two stress values in GR: 22.4, HC: 14.9, GS: 22.3 and WD: 14.7) (Figure 2).For example, in HC, sampling units were clearly differentiated into three parts, indicating seasonal effects.The sampling units in the spring (green colored symbols in Figure 3) were located at the lower-left of the ordination map, the units in summer (sky-blue colored symbols) were located in the upper part and the units in autumn (plum colored symbols) were ordinated towards the lower right.In GS, seasonal effects in sampling units were shown according to axis 1.The units in autumn were mainly located in the right part of the NMS, the units in spring were in the middle and lastly, the units in summer were located in left part in the NMS.In the NMS ordination for each stream over different seasons, the sampling units were ordinated mainly according to site differences, especially in summer (Figure 3).For example, in WD in summer, the sampling units at WD1 were mainly located in the upper parts of the ordination, the units at WD2 were in the left part and the units at WD3 were in the right part.In GR, based on the axis 2, the units in GR1 were located in lower parts whereas the units in GR2 and 3 were ordinated in upper parts.In addition, based on axis 1, the units in GR2 were in the left parts whereas the units in GR3 were in the right parts of the NMS.
In the NMS ordination for each stream, the distribution of the sampling units reflected seasonality rather than the differences among sites (first two stress values in GR: 22.4, HC: 14.9, GS: 22.3 and WD: 14.7) (Figure 2).For example, in HC, sampling units were clearly differentiated into three parts, indicating seasonal effects.The sampling units in the spring (green colored symbols in Figure 3) were located at the lower-left of the ordination map, the units in summer (sky-blue colored symbols) were located in the upper part and the units in autumn (plum colored symbols) were ordinated towards the lower right.In GS, seasonal effects in sampling units were shown according to axis 1.The units in autumn were mainly located in the right part of the NMS, the units in spring were in the middle and lastly, the units in summer were located in left part in the NMS.In the NMS ordination for each stream over different seasons, the sampling units were ordinated mainly according to site differences, especially in summer (Figure 3).For example, in WD in summer, the sampling units at WD1 were mainly located in the upper parts of the ordination, the units at WD2 were in the left part and the units at WD3 were in the right part.In GR, based on the axis 2, the units in GR1 were located in lower parts whereas the units in GR2 and 3 were ordinated in upper parts.In addition, based on axis 1, the units in GR2 were in the left parts whereas the units in GR3 were in the right parts of the NMS.

Discussion
Headwater streams are highly heterogeneous environments [9,10,26,55], supporting unique faunas that can differ from those in larger downstream areas [11].Further, spatial and seasonal variations of various environmental factors create complex habitat conditions [56].Upstream diversity influences the diversity of species found downstream and thus is important for the reestablishment of populations following local extinction events [57].Despite the importance of headwater ecosystems for the resilience of species diversity upstream and downstream, little attention has been given to scale-dependent or multi-scale dependent variability in macroinvertebrate communities in headwater streams [58,59].
Our results showed that community indices were significantly different between seasons and sites that were closely located geographically (<500 m).The differences in species richness at the local scale could be caused by local processes such as habitat heterogeneity [60], biotic interactions [61], and biogeographical processes [62].Moreover, because all the riparian zones were predominately forested, with no anthropogenic disturbance, the main factors differentiating the community composition at the stream and site scales likely relate to the natural variability of physical habitats and seasonal changes (e.g., canopy cover and the degree of autumn-shed leaves) [63].For example, the differences in riparian vegetation, latitude, discharge rate and substrate composition prevailing among riffles and/or sites in each stream sections influence the distributions of macroinvertebrates.The amount, magnitude, and intensity of precipitation could also differ between headwater streams, reflecting regional differences (i.e., southern and northern regions in Korea) [64].In addition, each season can harbor unique habitats with interactions among differential environmental factors and organisms.Periphyton biomass can be limited by light in autumn and summer but not in spring,

Discussion
Headwater streams are highly heterogeneous environments [9,10,26,55], supporting unique faunas that can differ from those in larger downstream areas [11].Further, spatial and seasonal variations of various environmental factors create complex habitat conditions [56].Upstream diversity influences the diversity of species found downstream and thus is important for the re-establishment of populations following local extinction events [57].Despite the importance of headwater ecosystems for the resilience of species diversity upstream and downstream, little attention has been given to scale-dependent or multi-scale dependent variability in macroinvertebrate communities in headwater streams [58,59].
Our results showed that community indices were significantly different between seasons and sites that were closely located geographically (<500 m).The differences in species richness at the local scale could be caused by local processes such as habitat heterogeneity [60], biotic interactions [61], and biogeographical processes [62].Moreover, because all the riparian zones were predominately forested, with no anthropogenic disturbance, the main factors differentiating the community composition at the stream and site scales likely relate to the natural variability of physical habitats and seasonal changes (e.g., canopy cover and the degree of autumn-shed leaves) [63].For example, the differences in riparian vegetation, latitude, discharge rate and substrate composition prevailing among riffles and/or sites in each stream sections influence the distributions of macroinvertebrates.The amount, magnitude, and intensity of precipitation could also differ between headwater streams, reflecting regional differences (i.e., southern and northern regions in Korea) [64].In addition, each season can harbor unique habitats with interactions among differential environmental factors and organisms.Periphyton biomass can be limited by light in autumn and summer but not in spring, while nutrients can limit periphyton when light availability is higher [30].Furthermore, seasonality in hydrology can be influential to structure macroinvertebrate composition [65].During spring, snow-melting can be the main source of surface water supply as well as groundwater recharge.Particularly, in Korea, sequential floods (i.e., summer) and droughts (i.e., autumn) are main natural disturbances in headwater streams that affect the composition of benthic macroinvertebrate composition [66].
Differences within FFGs and HTGs were also observed among sites and seasons.Taxa associated with a particular habit category (i.e., HTGs) exhibit certain morphological, physiological and behavioral adaptations to various microhabitats in freshwater ecosystems [67].They can exist at low discharge rates compared with areas downstream because headwater streams are generally supplied by small catchment areas [64].Clingers have morphological adaptations (e.g., curved tarsal claws, dorsoventral flattening, ventral gills arranged as a sucker, suction discs, and use of silk to construct attached retreats) that allow them to cling to substrate surfaces [68].Therefore, in this study, the differences between hydrological variables as well as substrate compositions may have caused the significant differences in the abundance of clingers among streams.Furthermore, scrapers showed differences among streams and sites over time compared to other FFGs in this study.This was likely due to the differences in stream width and canopy cover.For example, the distribution of grazing invertebrates is directly influenced by the distribution of benthic algae, and therefore indirectly influenced by canopy cover [69,70].Many researchers have suggested that scraper abundance tends to exhibit small-scale patchiness, resulting in localized variations depending on their algal food resources [29,30].
In our study, in NMS, samples were differentiated by seasons more than by spatial differences in each headwater stream.Within each season, the longitudinal differences in benthic macroinvertebrate communities were reflected in the NMS ordination.The units were clearly differentiated according to site differences even though the ordination patterns in each season were dissimilar.This indicated that in spite of their short distances between the adjacent sites in each stream (i.e., less than 500 m) without anthropogenic disturbances, they have their own habitat characteristics among sites, which have different resilience and resistance in comparison to seasonal effects, reflecting complicated interactions among spatial and temporal cues.

Conclusions
Our study examined the structure and function of the macroinvertebrate community at two different spatial scales during three seasons.Community and functional diversity indices varied significantly within seasons and/or sites as well as by the category of FFGs or HTGs.In NMS, within a single headwater stream, samples were separated by seasonality rather than spatial differences.Within each season, sample ordination reflected site differences, suggesting that macroinvertebrate communities respond to multiple and interacting spatial and temporal cues.Therefore, continuous monitoring and research on the interactions between species diversity and spatio-temporal and physiochemical effects are fundamental to maintain catchment biodiversity and to provide strategies for watershed restoration of macroinvertebrate communities.

Figure 1 .
Figure 1.Locations of the sampling sites in four different headwater streams.

Figure 1 .
Figure 1.Locations of the sampling sites in four different headwater streams.

Figure 2 .
Figure 2. Spatial and/or temporal changes in macroinvertebrate communities using NMS ordination in four different headwaters.Acronyms in NMS units stand for the samples: the first numbers indicate sampling sites (i.e., 1, 2 and 3) in each headwater and the last numbers represent replicates in each sampling site (1, 2, 3, 4 and 5).Each axis was rescaled on the 0-100 range based on the min-max scores of the NMS axes.(The stress values of the first two axes at GR: 22.4, HC: 14.9, GS: 22.3 and WD: 14.7).

Figure 2 .
Figure 2. Spatial and/or temporal changes in macroinvertebrate communities using NMS ordination in four different headwaters.Acronyms in NMS units stand for the samples: the first numbers indicate sampling sites (i.e., 1, 2 and 3) in each headwater and the last numbers represent replicates in each sampling site (1, 2, 3, 4 and 5).Each axis was rescaled on the 0-100 range based on the min-max scores of the NMS axes.(The stress values of the first two axes at GR: 22.4, HC: 14.9, GS: 22.3 and WD: 14.7).

Figure 3 .
Figure 3. Spatial differences in macroinvertebrate communities using NMS ordination in four different headwaters in each season.Acronyms in NMS units stand for the samples: the first numbers indicate sampling sites (i.e., 1, 2 and 3) in each headwater and the last numbers represent replicates in each sampling site (i.e., 1, 2, 3, 4 and 5).Each axis was rescaled on the 0-100 range based on the minmax scores of NMS axes.(The stress values of the first two axes at GR: spring 14.8, summer 14.8, and autumn 10.3, at HC: spring 13.4, summer 18.4, and autumn 12.3, at GS: spring 10.5, summer 16.0, and autumn 16.4, and at WD: spring 16.8, summer 16.4, and autumn 104).

Figure 3 .
Figure 3. Spatial differences in macroinvertebrate communities using NMS ordination in four different headwaters in each season.Acronyms in NMS units stand for the samples: the first numbers indicate sampling sites (i.e., 1, 2 and 3) in each headwater and the last numbers represent replicates in each sampling site (i.e., 1, 2, 3, 4 and 5).Each axis was rescaled on the 0-100 range based on the min-max scores of NMS axes.(The stress values of the first two axes at GR: spring 14.8, summer 14.8, and autumn 10.3, at HC: spring 13.4, summer 18.4, and autumn 12.3, at GS: spring 10.5, summer 16.0, and autumn 16.4, and at WD: spring 16.8, summer 16.4, and autumn 104).

Table 1 .
Average (standard deviation) of physico-chemical characteristics of headwater streams.
Values in parentheses for each headwater stream indicate the sampling year.

Table 3 .
Summary of two-way ANOVAs for community indices at different sites and seasons.
Df: degree of freedom and MS: mean square.

Table 4 .
Summary of two-way ANOVAs for functional feeding groups at different sites and seasons.
Df: degree of freedom and MS: mean square.

Table 5 .
Summary of two-way ANOVAs for habit trait groups at different sites and seasons.
Df: degree of freedom and MS: mean square.