Impacts of streamflow alteration on benthic macroinvertebrates by mini-hydro diversion in Sri Lanka

<jats:title>Abstract</jats:title><jats:p>Our study focused on quantifying the alterations of streamflow at a weir site due to the construction of a mini-hydropower plant in the Gurugoda Oya (Sri Lanka), and evaluating the spatial responses of benthic macroinvertebrates to altered flow regime. The HEC–HMS 3.5 model was applied to the Gurugoda Oya sub-catchment to generate streamflows for the time period 1991–2013. Pre-weir flows were compared to post-weir flows with 32 Indicators of Hydrologic Alteration using the range of variability approach (RVA). Concurrently, six study sites were established upstream and downstream of the weir, and benthic macroinvertebrates were sampled monthly from May to November 2013 (during the wet season). The key water physico-chemical parameters were also determined. RVA analysis showed that environmental flow was not maintained below the weir. The mean rate of non-attainment was ~ 45% suggesting a moderate level of hydrologic alteration. Benthic macroinvertebrate communities significantly differed between the study sites located above and below the weir, with a richness reduction due to water diversion. The spatial distribution of zoobenthic fauna was governed by water depth, dissolved oxygen content and volume flow rate. Our work provides first evidence on the effects of small hydropower on river ecosystem in a largely understudied region. Studies like this are important to setting-up adequate e-flows.</jats:p>

section_content: The alteration of river flow regimes is claimed to be the most serious and continuing threat to the ecological sustainability of riverine environments 1, 2. The duration and seasonal timing of associated low flow conditions, along with the reduced flow variability, strongly influence riverine organisms directly and via changes to habitat 2, 3. It is thus important to provide means of ensuring future developments that are sustainable and able to protect biological richness and ecosystem functions. This has led to a rapid increase of studies aimed to quantitatively understand aquatic ecosystem responses to various degrees of flow alteration 4, 5. 

Small hydropower projects offer one of the most promising energy resources for long-term sustainable development in Sri Lanka. Small hydropower plants have been, and to some extent still are, viewed as an environmentally benign energy source, and are categorized by the Sustainable Energy Authority of Sri Lanka as a green and renewable technology 6. Small hydropower can, however, exert multiple impacts on local environment (e. g. alteration of flow regime with consequent changes of water physico-chemical parameters and habitat structure), and the impacts which are perceived to be of critical importance are ecological, centered on aquatic flora and fauna. The underlying cause could be attributed to the non-maintenance below the weir of sufficient 'environmental flow' (e-flow, i. e. a natural flow paradigm comprising the five components of magnitude, frequency, duration, timing and rate of change), which is recognized as the key to sustaining biodiversity and ecosystem integrity 7. Previous works on the effects of water depletion on benthic macroinvertebrates in river reaches downstream from diversion structures have frequently pointed out significant impairment, particularly in terms of richness reduction related to habitat trivialization 3. The impact is higher in river reaches where hydrologic alteration is more severe, i. e. where the entire flow or a very large proportion of it is diverted [8] [9] [10], and where flow alterations cause relevant habitat changes in terms of substrate heterogeneity reduction, nutrient enrichment and temperature regime alteration [11] [12] [13]. Biological responses to flow alterations thus strictly depend on the context but also vary according to the kind of off-stream diversion scheme. To date, much of the research has been conducted on dams and impoundments [14] [15] [16] [17], mainly of large size, and few studies have assessed the ecological effects of run-of-theriver schemes 18, 19 and minor intakes [20] [21] [22]. 

Benthic macroinvertebrates are an important component of the river biota and are indicators of river health 23, 24. There is a wealth of literature suggesting that macroinvertebrate community composition is tightly linked to instream hydraulic conditions [25] [26] [27]. The knowledge of associations between environmental factors and zoobenthic assemblages is essential in understanding how aquatic communities in a particular geographic area are structured by the physical and chemical make-up of their environment, how they are affected by alterations to those conditions and, as a consequence, how the health of the entire riverine ecosystem is determined. This knowledge is the base to improve managerial decisions in water resource governance. 

A growing need to predict the biological impacts associated with water management activities and to set water management targets that maintain the integrity of riverine ecosystems has created the scientific discipline of 'in-stream flow' modelling and design. The primary application of in-stream flow models has been the design of e-flow regimes to guide sustainable water abstractions. Thus, it is evident that decisions on water regulation projects would benefit if they were informed by quantitative predictions of the ecological effects of varying degrees of streamflow alteration. Richter et al. 7, 28 developed and demonstrated the Range of Variability Approach (RVA), a holistic e-flow methodology, for establishing river management targets by incorporating the concept of natural hydrologic variability. RVA accepts that it is not possible to maintain the full range of natural streamflow variability in regulated or otherwise affected river systems, but supports efforts to manage hydraulic alterations in a manner that minimizes impacts on natural hydrologic variability and advocates conservation of native aquatic biodiversity and protection of natural ecosystem functions. In recent years, this approach or a revised version was successfully applied by many authors [29] [30] [31] [32]. 

Most developing countries still lack the technical and institutional capacity to establish environmental water allocation practices and policies 33. The existing methods of e-flow assessment are either complex and resourceintensive or not tailor-made for the specific conditions of a country, region or basin 34. This is the situation in Sri Lanka. In this country, hydroelectricity is the oldest and main source of electricity generation, with a share of nearly 45% of the total available grid capacity in 2010 6. However, there exists no lawful regulation stipulating the exact quantity of streamflow that should be released as e-flow. In this context, RVA could allow to evaluate the post-diversion circumstances in relation to pre-diversion flows and set management targets to maintain the riverine ecosystems. Since the e-flow concept is relatively new to Sri Lanka, it is also of importance to create awareness among responsible authorities through relevant studies on e-flow assessment. 

Therefore, in our study, a holistic e-flow approach for determining streamflow requirements to sustain benthic macroinvertebrate communities and consequently river health below the weir of a small (approx. 5000 kW) hydropower plant in the Gurugoda Oya (a river in Sri Lanka) was explored. Specific aims of the study were (1) to calibrate and validate Hydrologic Engineering Centre-Hydrologic Modelling System (HEC-HMS) 3. 5 model and generate long-term flow data for the assessment of post-weir hydrologic alterations, (2) to investigate the spatial variation of water physico-chemical parameters and zoobenthic assemblages above and below the weir, and (3) to identify a-posteriori which physico-chemical parameters are useful to predict the macroinvertebrate community composition through the application of Distance-based Linear Model (DistLM). 

We expected that flow abstraction could determine marked changes of key environmental factors which in turn could be associated to biodiversity reduction. We carried out our work during the wet season when the differences due to flow abstraction between impacted and unimpacted sites could be greatest.

section_content: Degree of hydrologic alteration due to mini-hydro diversion. Figure 1a and b show the goodness of fit of simulated flow values against observed flow values for the calibration and validation of HEC-HMS 3. 5 model for the Holombuwa catchment. Calibration for the time period 1991-2001 yielded a best scenario of 75. 8% of residual points falling within ± 1 SD range, and 96. 1% within ± 2 SD range. Moreover, a R 2 value of 0. 66 and a normal distribution of residuals were detected. However, the model slightly over-predicted flows at comparatively lower monthly flow ranges (10-15 m 3 /s). Percentage of residual points within ± 1 SD and ± 2 SD, and R 2 value of validation results for the next 12 years (2002-2013) were 78. 9%, 96. 5% and 0. 67, respectively. Also the comparison of residuals between observed and simulated flows for the entire study period (1991-2013) yielded results (percentage of residual points within ± 1 SD and ± 2 SD, and R 2 value of 79. 3%, 96. 7%, and 0. 66, respectively) above the limits suggested by Mood et al. 35, indicating that the validity of the model held even for long periods of runoff simulations in the Holombuwa catchment. The model could thus be successfully applied to the Gurugoda Oya sub-catchment, and RVA targets and rate of non-attainment for the 32 considered Indicators of Hydrologic Alteration (IHAs) were calculated (Fig. 1c and Table 1 ). Due to the skewness in the distribution of the pre-weir annual values for certain indicators, the mean -1 SD values fell outside (below) the pre-weir low range limits. For these parameters, the pre-weir minima of their range were selected. 

The rate of non-attainment of the group 1 IHAs, which represent the magnitude and timing of flows, was above 30% for the post-weir period. During certain months the rate of non-attainment was 100%. The decrease in the monthly mean flows suggests a drastic drawdown in the water table in the areas downstream of the weir. Group 2 parameters indicated that the daily, weekly, monthly and quarterly minimum flows were negatively influenced by weir regulation. Relevant alterations were also observed in parameters of groups 4 and 5, which represent the timing and frequency, and rate of changes of flow regimes, respectively. Rate of non-attainment of the group 4 IHAs varied between 0 and 100%, and of group 5 IHAs mostly between 67 and 100%. 

In summary, out of the 32 IHAs, 11 parameters scored a rate of non-attainment of 33%, and 12 parameters between 50 and 100%. The calculated mean rate of non-attainment of the flow of the Gurugoda Oya below the weir was around 45% (i. e. moderate alteration). Our results proved that the IHA method encasing RVA targets is an easy and useful tool to quantify the hydrologic alteration in the study area as already reported for many rivers worldwide [36] [37] [38] [39] [40]. However, these results were preliminary since only three years of data (2011-2013) were used for IHA calculation for the post-weir period. Hence, the high level of fluctuation of rate of non-attainment scored within a single IHA category. It is expected that with time, the fluctuations will decrease and a more reasonable outcome be proposed. However, the implementation of the RVA at an early stage of weir operations could set baselines upon which future river management decisions could be taken and weir operations be performed. For example, the fact that certain parameters scored a rate of non-attainment of 100%, suggests a significant level of hydrologic alteration and calls for immediate changes in weir operations. Thus, the application of RVA to the weir site could be seen as a timely approach to be used as the base for future analyses.
section_name: Results and discussion

section_content: A significant difference of the measured physico-chemical parameters (except for pH and five-day Biological Oxygen Demand-BOD 5 ) was recorded between study sites upstream (F-E-D) and downstream (C-B-A) of the weir (one-way ANOVA, p < 0. 05, and Tukey test for pairwise comparison between sites F-E-D and C-B-A, p < 0. 05). The most prominent feature was the drastic reduction in mean flow velocity and volume flow rate values at the sites below the diversion point (Table 2 ). It was also evident that, within the sampling period, the volume flow rate above the weir fluctuated heavily whilst below the weir fluctuated within a very narrow range. This shows that the flow released downstream of the weir was heavily regulated irrespective of the rainfall to the locality. Hydropower exploitation frequently induces reduction in streamflow magnitude, which is a strong predictor of biological integrity 41. However, the impact on aquatic communities is expected to be lower in case of run-of-the-river schemes 18, 19 and minor intakes [20] [21] [22] than in case of reservoirs [14] [15] [16]. In the latter case, more severe reductions of the flow variability and changes of flow timing add to the substantial decrease of the flow magnitude. 

During the sampling period (May-November), 16 benthic macroinvertebrate taxa belonging to three phyla, i. e. Annelida (Oligochaeta), Mollusca (Gastropoda and Bivalvia) and Arthropoda (Malacostraca and Insecta), were found at the six sampling sites (Table 3 ). 

Average taxon richness recorded at sites above the weir (12-13) was higher than that at sites below the weir (7-9). However, the cluster analysis detected three major groups at around 50% of similarity (Fig. 2 ). Sites above the weir differed from sites below the weir, and, among the latter sites, sites A and B differed from site C, i. e. the site closest to the weir (one-way ANOSIM, p<0. 05). Also the NMDS plot (Fig. 3 ) shows how the assemblages collected at site C differed from the assemblages collected at the other sites, with a higher level of intra-site variation between sampling occasions. Moreover, the NMDS plot gives a visual representation of how clustering becomes complex as the similarity level increases. The average abundance of all taxa, except for most of the insect taxa and Melanoides sp. were higher at the upstream sites than at the downstream sites. Gyraulus sp., Lamellidens sp., Paludomus sp., and Heleocoris bengalensis were even absent at the downstream locations. Pila sp. dominated instead all the six study sites except for site C where insect nymphs and larvae displayed similar relative abundances and the lowest average value of total density was observed (Table 3 ). This pattern was detected Table 3. The overall average composition (AA = absolute abundance-n. individuals, RA = relative abundance-%) and distribution of benthic macroinvertebrate communities in the Gurugoda Oya, above (sites F-E-D) and below (sites C-B-A) the weir (see Fig. 5 ). by Shannon-Wiener and Pielou's evenness index, both recording the highest average values at site C (1. 86 and 0. 85 respectively). Average values highly (0. 44 and 0. 20) or slightly (0. 78 and 0. 40) lower than those detected at the three upstream sites (diversity range 1. 03-1. 31, and evenness range 0. 42-0. 51) were instead recorded at sites A and B respectively. The analysis of benthic macroinvertebrate communities thus revealed marked dissimilarities between the sites downstream of the weir and the control sites, upstream of the weir. Mini-hydro diversion was associated to taxa loss and consequent richness reduction; only the assemblage detected at site C showed an increase of diversity and evenness. This variable result is common in studies on benthic macroinvertebrate responses to flow alterations. In different contexts, reduced streamflow is reported to induce both decrease and increase in benthos abundance, richness and diversity 15, 16, 21, 42, 43. DistLM (best solution: AICc = 281. 9, R 2 = 0. 48, RSS = 27844) identified flow depth, dissolved oxygen (DO) concentration, and volume flow rate (in order of decreasing importance) as the environmental variables driving differences in the structure of macroinvertebrate assemblages of the Gurugoda Oya (Fig. 4 ), and depicted how different sites are governed by different factors, i. e. site F by the volume flow rate, site D and E by the water depth and site C by the DO content. The dominance of Pila sp. may be due to the micro-habitat conditions created by varying flows and effects of drought during the dry season 44. The major diversity at site C could be attributed to the low depth of the site and the consequent mixing of the water which makes available sediment-bound nutrients into the water column besides increasing oxygen content and micro-habitat diversity. The habitat changes generated by streamflow reduction due to water abstraction can vary widely, depending on channel morphology and substrate stability, and on the possible related nutrient enrichment and temperature regime alteration [11] [12] [13] 25. How ver, the best fit model showed a cumulative percentage variation below 50% on the two major axis (Fig. 4 ). A major reason for this could be the contribution of other larger-scale environmental or biotic factors on the taxa to observed trends in abundance rather than reduced flow. Realistically, it is an interplay of numerous physico-chemical and biotic factors in the hydro-climatic region. For instance, Jayawardana et al. 45 recently highlighted that local riparian forest cover is important in structuring macroinvertebrate communities in the Uma Oya catchment. Also allochthonous nutrient inputs are reported to influence zoobenthic assemblages of Eswathu Oya and Yan Oya 46. Unfortunately, few other studies on benthic macroinvertebrate communities of Sri Lanka rivers are available 47, 48, and none on the effects of hydropower. Direct comparisons of our data with literature data were thus not feasible, also due to the different resolution level used for taxa identification and sampling methods 49. To our knowledge, our results provide first evidence for a Sri Lanka river of associations between changes in environmental factors and zoobenthic assemblages consequent to water abstraction due to small hydropower. This is the first step towards the acquisition of the information required to setting-up adequate e-flows [50] [51] [52] [53].
section_name: Effects of weir on water physico-chemistry and benthic macroinvertebrates.

section_content: The HEC-HMS 3. 5 model can be endorsed to be reliably used to simulate Gurugoda Oya flows with proper calibration and validation. As the transformation method in the model, 'snyder unit hydrograph method' could be recommended for the Gurugoda Oya basin with the 'initial and constant rate loss' as the loss method. 

The RVA showed that the e-flow was not maintained below the weir throughout 2011-2013. The level of alteration of flows caused by hydropower plant operations was concluded to be moderate (45%). Our results also revealed that the prevailing physico-chemical parameters as well as zoobenthic assemblages varied significantly among the study sites located up-and downstream of the weir. Differences in macroinvertebrate assemblage structure were associated to water depth, dissolved oxygen concentration and volume flow rate. 

Although further research is needed to investigate the impact of mini-hydro diversion at a larger spatial and temporal scale, our study should help to pay attention to this relevant and increasing issue in a largely understudied region such as Sri Lanka. We expect that this finding may be used by Government authorities and other policy makers involved in disciplines related to water governance (e. g. Irrigation Department, Sustainable Energy Authority, Ceylon Electricity Board) to make informed decisions on licensing procedures of mini-hydropower plants, and decisions focused on sectoral water allocations for sustainable development while maintaining the integrity of riverine ecosystems.
section_name: Conclusions

section_content: Study area. The study area is located within Kegalle district in Sri Lanka. It lies entirely within the 'low country wet' agro-ecological zone with an expectancy of annual rainfall of 1900-3200 mm 54. The dry season (rainfall range 50-150 mm) is between November and April, while the wet season (rainfall range 300-500 mm) due to the south-west monsoon between May and October. The mean annual temperature is 27. 8 °C and the 2 ) among the six sampling sites (sites F-E-D are located above the weir and sites C-B-A are located below the weir, see Fig. 5 ). mean relative humidity is approximately 80%. Rubber plantations and home gardens cover most (74%) of the area. Elevation ranges between 20 and 1240 m asl. Our study was carried out at the Hungampola South/Morontota village areas located mid-east in the Gurugoda Oya sub-catchment which was created within the main Holombuwa catchment (Fig. 5 ). The reason for the division of a 'Holombuwa catchment' and a 'Gurugoda Oya sub-catchment' (catchment in relation to the weir site of the mini-hydropower plant) was due to the availability of flow data. Since only the Holombuwa catchment has a flow gauging station, the HEC-HMS 3. 5 model was calibrated and validated to this catchment and applied to the Gurugoda Oya sub-catchment to simulate longterm flows at the weir site. 

Model calibration, validation and application. Daily precipitation and monthly evaporation data (1991-2013) at the Undugoda rain gauging station located within the Holombuwa catchment (Fig. 5a ), along with daily streamflow data (1991-2013) at the Holombuwa gauging station were obtained from the Hydrology Division of the Department of Irrigation. 

For HEC-HMS 3. 5 model simulation, the 'initial and constant rate loss method' was selected as the loss method and the 'Snyder unit hydrograph method' was selected as the transform method based on the model calibration and validation reported by Halwatura and Najim 55. The Holombuwa catchment was used to accomplish the calibration goal. Daily rainfall data for 11 years (1991-2001) were used in the model and the flows simulated from each calibration run were tested statistically against the actual measured flow to produce a best fit model. The calibrated model was then applied for the period 2002-2013 to accomplish the validation goal and the simulated flows were statistically compared with observed flows for the same time period. Concurrently, a simulation run was also performed with the entire data set (1991-2013), for a further validation. The model parameters were adjusted by certain percentages until the statistical evaluation resulted in more than 68% of the residual points (observed value -simulated value) falling within ± 1SD, more than 95% falling within ± 2 SD and a R 2 value nearest to 1. According to Mood et al. 35 a successful model calibration presents the mentioned thresholds. The best parameter distribution scenario was chosen for the subsequent simulation processes. The model was applied to the Gurugoda Oya sub-catchment to generate the daily runoff values for the past 23 years (1991-2013), with the catchment outlet defined as the weir site. The flows which prevailed below the weir (construction period 2006-2010) for the time period 2011-2013 were obtained from data logs of studies conducted in the region and of this study field-sampling campaign.
section_name: Methods

section_content: The holistic e-flow assessment methodology used in this study is based on the RVA developed by Richter et al. 7. In the RVA, the pre-weir streamflow regime (1991-2006) was compared with the post-weir flow regime (2011-2013) using 32 ecologically relevant hydrologic parameters (IHAs) (Table 1 ). 

Measures of central tendency (means) and dispersion [range limits (low and high) and standard deviation] were calculated from the pre-weir annual series for each of the 32 parameters, which produced 64 inter-annual statistics for each annual data series (32 measures of central tendency and 32 measures of dispersion), which were used to characterize inter-annual variations. Values at ±1 SD from the mean were selected as the RVA targets (lower and upper RVA limits) for each of the 32 IHAs as recommended by Richter et al. 7. 

The degree of hydrologic alteration or rate of non-attainment of each hydrologic parameter (values that fall below the lower limit and above the upper limit of calculated RVA targets) after the construction of the weir, was then calculated using the following equation 28 : D is the degree of hydrologic alteration/rate of non-attainment, N o is the observed number of post-impact years for which the parameter value falls within the RVA target range, N e is expected number of post-impact years for which the parameter value falls within the RVA target range. N e can be estimated by P x N T where P is the percentage of pre-impact years for which the parameter value falls within the RVA target range and N T is total number of post-impact years. Values between 0 and 33% represent little or no alteration, 33-67% moderate alteration, and 67-100% high alteration. Finally, the degrees of hydrologic alteration of all 32 parameters were averaged to obtain a single level of alteration of flow regime for the weir site. The IHA software (version 7. 0) developed by the Nature Conservancy was used to calculate the 32 IHAs, the RVA targets and the levels of alterations of the parameters. 

Collection of benthic macroinvertebrates and water physico-chemical parameters. Six sampling sites were selected in the study area to capture the effects of different flow regimes on benthic macroinvertebrates. Sampling sites F-D-E and C-B-A were established separately 300 m away from each other along the Gurugoda Oya, the former three above and the latter three below the weir (Fig. 5b ). The study was carried out from May to November 2013 with monthly intervals between sampling occasions. 

At each sampling site, the cross section of the stream was divided into three equal parts and sampling was carried out in the center of these three sections separately. The three replicates were then pooled into one integrated sample. Macroinvertebrates were collected using a standard D-framed dip net consisting of a D-shaped metal frame (0. 3 m width and 0. 3 m height) holding a conical net (mesh aperture 400 μm). A dip and sweep method was employed where organisms were collected by aggressively disturbing the target habitat. A sweep of 0. 5 m length was made per sampling effort. The net was dipped into the substrate and three such sweeps were performed 56 to collect bottom sediments covering an area of about one square meter 57. River substrate varied from pebbles to boulders in different percentages with a local presence of soft sediment. The depth of each sampling point was measured using a pole and tape in tandem with benthic sampling. Later in the laboratory, the material retained
section_name: Calculation of Indicators of Hydrologic Alteration and rate of non-attainment.

section_name: Data availability

section_content: Main data generated and analyzed during the current study are included in the manuscript, further information is available from the corresponding author on reasonable request.
section_name: Data availability

section_name: Data analysis.

section_content: www. nature. com/scientificreports/ was wet sieved through a mesh (0. 5 mm aperture) 57 and identified to the nearest possible taxonomic category using the naked eye and a binocular microscope following the standard identification keys provided by Mendis and Fernando 58 and Starmühlner 59. At each sampling site, each time a biological sample was taken, the physico-chemical parameters of the water immediately above the bottom were measured. Such parameters were temperature, pH, conductivity, Total Dissolved Solids (TDS), Dissolved Oxygen (DO), five-day Biological Oxygen Demand (BOD 5 ) and water flow velocity. Temperature, pH, conductivity and TDS were measured in situ using the Yellow Springs Instrument (YSI) 55 water quality logger (USA). The DO concentration was measured in the laboratory using the Winkler method 60, after collecting water samples into amber-colored glass bottles (250 mL) and fixing DO in situ using manganous sulfate and Winkler reagents. Additional water samples were collected into amber-colored glass bottles (250 mL), brought to the laboratory, and incubated for five days at room temperature in total darkness. After the five-day incubation period, DO concentration in each bottle was measured using the Winkler method 60 and BOD 5 was determined. The water flow velocity was measured by using a float (a piece of Styrofoam) to drift for a known distance along the water current for a known period of time. Each water parameter was measured in triplicate and the mean value was calculated later on.

section_content: Similarities of macroinvertebrate assemblages were assessed using Bray-Curtis similarity clustering method 61 on square root transformed abundance data. Moreover, Non-metric Multi-Dimensional Scaling (NMDS) 62 was performed to better represent the spatial and temporal clustering of the benthic macroinvertebrate communities between the study sites. Analysis of similarities (ANOSIM) was carried out to detect significant differences of community composition among the six sampling sites. Total density, taxon richness, Shannon-Wiener diversity index and Pielou's evenness index were then calculated for each sample. One-way analysis of variance (ANOVA) and Tukey post-hoc test were used to compare the considered physico-chemical parameters among the six study sites. 

Distance-based Linear Model (DistLM) 63 was applied to perform permutational regression between zoobenthic assemblages and environmental variables assessing the relative contributions of environmental variables structuring macroinvertebrate communities. Prior to DistLM, draftsman plots and correlation matrices were produced to assess the distribution of each variable and to identify co-correlating variables. Environmental variables were square root transformed to normalize their distribution, and since no pairs of variables had a Pearson's correlation coefficient larger than 0. 9, all were included in the analysis. The Bray-Curtis matrix on square root transformed macroinvertebrate abundances was used. DistLM was performed with selection based on the Akaike information criterion (AIC), step-wise selection procedure and 999 permutations. AIC was chosen as the method to create the most parsimonious model, as it adds a 'penalty' for increases in the number of predictor variables 63. Step-wise selection was chosen as it allows for both the addition and removal of a term to the model at each step 63. Distance-based Redundancy Analysis (dbRDA) plot was used to provide the best possible 2-dimensional visualization of DistLM result. 

All the statistical analysis were carried out using PRIMER V 6. 1. 16 (equipped with PERMANOVA+ V 1. 0. 6) statistical software.
section_name: Data analysis.

section_content: D. S. N. M., M. M. M. N. and M. M. M. conceptualized and designed the study, analyzed the data and drafted the manuscript. S. Q. contributed to data interpretation, and manuscript writing and editing. All authors reviewed and approved the manuscript.
section_name: Author contributions

section_content: The authors declare no competing interests.
section_name: Competing interests