How fish traits and functional diversity respond to environmental changes and species invasion in the largest river in Southeastern China

Study area and field sampling

The Min River is located in Southeastern China, ranging from 116.38°E to 119.72°E and 25.38°N to 28.32°N. The river runs through Fujian Province with a length of 584 km and a basin area of about 61,000 km². Both sampling in 1979 and 2015 were conducted by traditional gill nets fishing (20 × 2.5 m, 3 layers with mesh sizes of 1, 2, and 3 cm) at the same 24 sites in May from the upstream to the lower reach of the river (Fig. 1). Three gill nets were placed in the left, middle, and right parts of the river at each site. The lead line settled the net on the river bottom, and the float line kept it vertical in the water column, allowing it to capture fish as they swim. Gill nets were set before sunset and lifted after sunrise to cover maximal peaks of fish activity. According to the experts who participated in the sampling in 1979, electrofishing was also used in 1979, but detailed information was not documented. As suggested by these experts, additional electrofishing (backpack electrofishing unit, Model: CWB-2000 P, Yufengda, China; 12-V import, 250-V export) which lasted 30 min and extended 50 m along the river course was applied in 2015 in wading areas at each site to assure that sampling effort in 2015 was at least as exhaustive as that in 1979. Species in 1979 and 2015 were identified and validated according to FishBase to avoid invalid species, synonyms, and homonyms (Froese & Pauly, 2019).

Permission to conduct the study was granted by the Animal Care and Ethics Committee of Ocean University of China (No. 20150305003). Field experiments in 2015 were approved by Research Council of Ocean University of China (No. 20150305004).

Environmental variables

We considered eight environmental variables, covering climatic variables (temperature and precipitation), human factor (human population), geographic variables (altitude and slope), river structure variables (Strahler order of stream segment and downstream link), and presence or absence of species invasion. Temperature (°C) and precipitation (mm) in 1979 and 2015 were derived from “China surface climatological data set (V3.0)” downloaded from China meteorological data service network (https://data.cma.cn/). Geostatistical spatial analyses were used to generate the inverse distance weight interpolation (IDW) map and describe the spatial variability in the mean annual temperature and precipitation at each site (Lu & Wong, 2008). Human population (thousand) at the county scale in 1979 and 2015 was derived from the China Population Statistics Yearbook (National Bureau of Statistics, 1980) and (National Bureau of Statistics, 2016). Altitude (m) and slope at 30 arc-second resolution were extracted from EarthEnv (https://www.earthenv.org/topography) (Amatulli et al., 2018). Two categorical variables were applied as the surrogates for spatial river structure, including Strahler order of stream segment (Order) (Strahler, 1957) and downstream link (Dlink) (Osborne & MJ, 1992). For Dlink variable, the magnitude of a link is firstly defined, which is the number of first-order segments upstream of a given point on a channel, then the Dlink at any point is the magnitude of the link below the next downstream confluence at any point. Species invasion was represented by a binary factor, where “1” refers to presence of non-native species, and “0” refers to absence of such species.

Measurement of functional traits

Ten morphological and ecological traits were selected for FD calculation based on two reasons: (1) these traits are associated with different functional groups, including feeding habit, trophic, swimming capability, habitat preference, and life cycle of the fish; and (2) these traits capture the response of fish species to abiotic and biotic environments, i.e., they are response traits (Violle et al., 2007). Morphological traits were relative head length, relative eye size, relative snout length, relative head depth, relative pectoral fin length, swimming factor, and relative body depth, while ecological traits were water flow preference, position in the water column, and maximum life span (Table 1) (Gatz Jr, 1979; Pool, Grenouillet & Villeger, 2014; Toussaint et al., 2016; Villeger, Grenouillet & Brosse, 2013; Watson & Balon, 1984).

Morphological traits were derived from basic morphological measurements, including standard length, head length, pectoral fin length, snout length, body depth, head depth, eye diameter, caudal peduncle depth, and caudal fin depth (Fig. S1). To mitigate the influence of intraspecific trait variations, at least 30 individuals of each species were required in both 1979 and 2015, and for those species with less than 30, all samples were measured. Samples subjected to measurements were individuals of similar size for each species, which also exceeded the minimum length of sexual maturity extracted from FishBase (Froese & Pauly, 2019) and “The Fish of Fujian Province” (Zhu, 1984). To mitigate the impact of deformations caused by dehydration in preserved samples on the accuracy of data in 1979, only intact individuals were measured and values of morphological traits for calculating FD were ratios of these basic morphological measurements. Ecological traits were extracted from FishBase (Froese & Pauly, 2019) and other ichthyological literatures (Chen, 1998; Chu, Zhen & Dai, 1999; Wu & Zhong, 2008; Yue, 2000; Zhu, 1984).

Functional diversity metrics

FRic, FEve, and FDiv were used because they measure size, regularity, and divergence of species in the functional space, respectively (Laliberté & Legendre, 2010; Villeger, Mason & Mouillot, 2008). Values of FEve and FDiv vary from 0 to 1, while values of FRic had no limit. Sørensen dissimilarity (FDsor) was used to calculate functional beta diversity (Villeger, Novack-Gottshall & Mouillot, 2011). It is decomposed into two components, turnover (FDsim) and nestedness (FDsne) (Villeger, Grenouillet & Brosse, 2013). By assuming that fish morphological traits in the Min River did not evolve at the temporal scale of decades, trait matrices in 1979 and 2015 were combined into a single matrix with mean morphological trait values of species. Then, this single matrix was subjected to Principal Coordinate Analysis (PCoA) with Gower’s distance using R package “ade4” (Anderson & Willis, 2003; Gower, 1966). Because FRic and FDiv both require more species than traits (Laliberté & Legendre, 2010) and the smallest species richness in this work was six, the first five PCoA axes, which explained 82.64% of cumulative variances, were preserved for FD calculations. Finally, R package “FD” was used to calculate FRic, FEve, and FDiv (Laliberté & Legendre, 2010), and “betapart” to calculate FDsor, FDsim, and FDsne (Baselga et al., 2020).

Paired t-test was used to test the difference of species richness, FRic, FEve, FDiv, FDsor, FDsim, and FDsne between 1979 and 2015. When paired t-test detected no significant difference, Levene’s test was used to test the difference of variance (Fox & Weisberg, 2019).

Random forest regression

Random forest (RF) regression was applied to explore the impacts of environmental factors on FRic, FEve, FDiv, FDsor, FDsim, and FDsne given that RF does not require a priori determination of probability density distribution (Breiman, 2001). Variance inflation factor (VIF) was applied to test multicollinearity among environmental factors with a stepwise procedure (Dormann et al., 2013). All selected environmental factors had a VIF below 5. Response variables in RF were resampled 5000 times to generate a large number of decision trees. The importance of an environmental factor to the model was calculated as how much the mean standard error (MSE) increased when this variable was permuted. The more the MSE increased, the more important this variable was. Ten repeats of five-fold cross-validation were applied to check the performance of the model, then the most important variables which yielded the least cross-validation error were identified. Partial dependence plots were used to visualize the response of FD to these important variables, which were classified into five categories: positive, negative, unimodal, bimodal, and irregular (Goldstein et al., 2015). R package “randomForest” was used for modeling (Liaw & Wiener, 2002), “usdm” for stepwise VIF testing (Naimi et al., 2014), and “pdp” for partial dependence plotting (Brandon, 2017).

Fourth-corner combined with RLQ

The association between functional traits and environmental variables was tested using the fourth-corner combined with RLQ approaches (Dray et al., 2014; Dray & Legendre, 2008). Species occurrence matrix (L) was ordinated by correspondence analysis (CA). Environmental variables matrix (R) was ordinated by Hill-Smith analysis for a mix of qualitative and quantitative variables (Hill & Smith, 1976). Traits matrix (Q) was ordinated by PCoA analysis with Gower’s distance as was used in FD calculation. Fourth-corner analysis was performed to correlate fish functional traits and environmental variables with the first two axes of RLQ. Following the proposal of Dray et al. (2014), this analysis combined permutation models 2 and 4 and used 9,999 repetitions to test the null hypothesis (H₀) that the species traits (Q) are not related to the environmental variables (R). In permutation model 2, sites were randomized to test the relationship between species occurrence (L) and the environmental variables (R), and in permutation model 4, species were randomized to test the relationship between species occurrence (L) and their traits (Q) (Dray & Legendre, 2008). H₀ is rejected when significant relationships are found in both permutation models. P values were adjusted using the false discovery rate (FDR) procedure to mitigate type I error (Benjamini & Hochberg, 1995). The R package “ade4” was used for the fourth-corner combined with RLQ analysis (Dray & Dufour, 2007).

Temporal changes of functional diversity

There were 125 species in 1979 and 82 species collected in 2015, 56 of which were shared by both periods. Species richness in 1979 was significantly higher than in 2015 (Table S1). FRic in 1979 was significantly higher than that of 2015 (Fig. 2A). Mean FRic in 1979 was 25.32 (± 3.35 SE), with the highest value at site 15, while mean FRic in 2015 was 7.00 (± 5.29 SE), with the highest value at site 21. FEve in 2015 was not significantly different from that of 1979, but Levene’s Test for equality of variances confirmed that variance of FEve in 2015 was significantly higher than that in 1979 (Fig. 2B). Paired t-test confirmed that FDiv in 1979 was significantly higher than that in 2015, except that FDiv at sites 4 and 13 were lower in 1979 than that in 2015 (Fig. 2C).

Decomposition of functional beta diversity revealed the main component driving FD dissimilarities in the river. FDsor (Fig. 3A) and FDsne (Fig. 3B) in 2015 were significantly higher than those in 1979, respectively, but FDsim (Fig. 3C) in 2015 was not statistically different from that in 1979. These results indicate that difference in FDsor between 2015 and 1979 was dominated by changes in FDsne.

Environmental factors driving functional diversity

RF models for FRic, FDsor, FDsne, and FDsim explained 83.83%, 88.78%, 83.42%, and 35.75% of the total variances, respectively (Table 2). Models for FEve and FDiv were not explanatory which explained less than 10% of the total variance. Ranking by their relative importance, precipitation, temperature, species invasion, and human population were the most important factors for modeling FRic, FDsor, FDsim, and FDsne, except that human population was not explanatory for FDsim. Partial dependence plot showed that responses of FRic to precipitation, temperature, species invasion, and human population were negative, while responses of FDsor, FDsim, and FDsne to precipitation, temperature, species invasion, and human population were positive, except that response of FDsim to human population was irregular (Fig. 4).

Correlations between functional traits and environmental factors

Global significance test rejected H₀ and confirmed that fish functional traits were significantly correlated with environmental factors (p <0.05 for both Model 2 and Model 4, Table 3). Axis 1 and Axis 2 of the RLQ analysis explained 90.40% and 6.54% of the total co-inertia, respectively. The first RLQ axis was negatively correlated with precipitation, temperature, and presence of invasive species, and positively correlated with absence of invasive species (Fig. 5A). Significantly associated traits with higher precipitation, higher temperature, and presence of invasive species were higher swimming factor (minimum caudal peduncle depth divided by maximum caudal fin depth), higher relative eye size (eye diameter divided by head depth), and lower relative pectoral fin length (pectoral fin length divided by standard length) (Fig. 5B). Traits associated with absence of invasive species, lower temperature, and lower precipitation were eurytopic water flow preference, higher relative pectoral fin length, lower swimming factor, and lower relative eye size. No significant associations were detected for the second axis, neither for environmental factors nor for traits.

Species associated with higher swimming factor and relative eye size were endemic species Sinibrama macrops, Pseudobagrus tenuis, and Pseudolaubuca sinensis, etc. (Fig. 5C). These species occurred in sites with higher precipitation, higher temperature, and presence of invasive species. Species associated with eurytopic water flow preference and higher relative pectoral fin length were nilem carp Osteochilus vittatus, Gobiobotia longibarba, and catfish Clarias fuscus, etc. These species were correlated with the opposite environmental conditions.

Sites in 2015 and 1979 were clearly separated on the ordination diagram (Fig. 5D). Sites in 2015 were characterized by environmental conditions with presence of invasive species, higher precipitation, and higher temperature, and associated with fish species with higher swimming factor, higher relative eye diameter, and lower relative pectoral fin length. Sites in 1979 were characterized by opposite environmental conditions and associated with species with eurytopic water flow preference and higher pectoral fin length. Meantime, sites in 2015 were more disaggregated than in 1979 on the ordination diagram, suggesting that trait compositions in 2015 were more varied than in 1979.

Temporal changes of functional diversity and community assembly

The Min River ever possessed diverse and distinctive fish species accounted by its heterogeneous microhabitats in history (Zhu, 1984). For example, Chinese sturgeon Acipenser sinensis and Chinese sucker Myxocyprinus asiaticus dwelled in the river according to historical records (Kang et al., 2014). Chinese sturgeon A. sinensis is listed as a critically endangered species on IUCN red list (IUCN, 2020). Although there have been some conservative attempts such as artificial propagation and release since 1980 to protect this species in the Min River (Kang, 2002), no individual was sampled in 2015. Chinese sucker M. asiaticus is an endemic species and quite popular in local aquariums in the Min River (Editorial Committee of Fujian Province Annals, 1999). This fish was collected at sites 11, 12, and 13 in 1979, but was absent from all these sites in 2015. Those 69 species that were collected in 1979 but absent in 2015 are all native to the Min River. For example, yellowcheek Elopichthys bambusa, long spiky-head carp Luciobrama microcephalus, and Ochetobius elongates are endemic species with high economic values (Editorial Committee of Fujian Province Annals, 1999). They have become rare nowadays in major rivers and lakes in China due to overfishing and habitat changes (Feng et al., 2014; Zhang et al., 2010).

FRic is known to be positively correlated with species richness (Villeger, Mason & Mouillot, 2008). Loss of native species with distinctive functional traits diminishes the functional diversity and leads to functional homogenization (Chua, Tan & DCJ, 2019). The temporal decrease of FRic in this work was directly linked to the reduction of species richness in 2015 comparing with 1979. Temporally homogeneous and stable habitat supports relatively higher FEve (Massicotte et al., 2014). On the other hand, intense land-use/cover changes under human interferences in some streams have been associated with a decrease in FEve (Barbosa, Pires & Schulz, 2020). Variance of FEve in 2015 was significantly higher than that in 1979, which indicates that fish functional traits in 2015 were distributed more unevenly than in 1979, and this unevenness possibly suggests temporally unstable habitat conditions.

FDsne refers to functional dissimilarity resulted from trait loss (Villeger, Grenouillet & Brosse, 2013). Higher FDsne in 2015 than 1979 reflects that smaller communities became functional subsets of larger ones as a consequence of trait loss. At local scales, ecosystems suffering from intense interferences manifest higher FDsor. For example, high-impact environments show higher FDsor and FDsne but lower FRic in cladoceran communities in Northeastern and Southern Brazil (Simoes et al., 2020). The dominance of FDsne in this work further indicates that loss of distinctive functional traits as a result of species loss might be accountable for the temporal pattern of functional beta diversity.

Species loss is associated with habitat deterioration in many cases, which strengthened environmental filtering. For example, local-scale factors related to habitat degradation were associated with fish trait convergence and decreases of fish species richness in a Mediterranean river (Radinger, Diego Alcaraz-Hernandez & Garcia-Berthou, 2019). FDiv in 2015 was significantly lower than that in 1979, indicating that more trait convergence, i.e., strengthened environmental filtering, occurred in 2015 than in 1979. Patterns in community structure and functional diversity of freshwater fish communities at global scales are correlated with abiotic or biotic filters (Comte et al., 2016). At local scales in the Min River, species loss and increased trait convergence indicate that temporally strengthened environmental filtering was the main temporal community assembly mechanism, which reconciles the fact that the Min River has been exposed to multiple interferences, including climate change, human disturbances, and species invasion.

Importance of environmental factors driving functional diversity

Global climate change impacts freshwater fish FD worldwide (Buisson et al., 2013). High fish FD loss in Southeast Asia, Indo-Burma, and South-Central China is accompanied by global climate change (Macdonald et al., 2020). Intensity of anthropogenic activities firmly associates with fish FD. For instance, increasing human pressure in Neotropical reservoirs is negatively correlated with fish FD (Dias, De Oliveira & Baumgartner, 2021). Species invasion also undermines regional fish FD by outcompeting local species (Milardi et al., 2019). There are over 400 species of non-native freshwater fish species in China, most of which were introduced exotically by aquarium and aquaculture industries (Xiong et al., 2015). Similar factors have been recognized in the Min River. For example, temperature and precipitation were increased by around 1.5 °C and 390 mm on average in 2015 compared with 1979 in the basin, respectively; human interferences, such as land-use/cover changes under urbanization, dam construction, and agriculture, have been intensified during the past decades (Ma et al., 2020; Ren et al., 2016; Zhao et al., 2016); several invasive species have been identified in the basin (Deng et al., 2020). The synergic effect of these factors might contribute to the loss of habitat suitability for fishes in the Min River.

There has been a growing amount of evidence that climatic change is a strong determinant for fish FD (Kuczynski et al., 2018). Fish habitat suitability is expected to shift under scenarios of global climate change for physiological, geological, or hydrochemical reasons (Marras et al., 2015; Xing et al., 2019). In the Min River, precipitation and temperature were the most influential factors driving the patterns of FRic, FDsor, FDsim, and FDsne. Decadal changes in temperature and precipitation possibly alter fish habitat conditions, including water temperature, duration of flood or drought season, water velocity, and water depth, which influence species distribution (Chang & Jung, 2010). For example, variations of temperature and precipitation in spring and summer changed fish communities in the Lower Colorado River Basin (Pool et al., 2010). This supports the result as previously discussed that environmental filtering is the dominant mechanism driving fish community assembly in the Min River temporally.

Several invasive species were collected in 2015, which occurred at 11 out of 24 sites in the river. While in 1979, none of these species were collected at these sites. These species are grass carp Ctenopharyngodon idella, Wuchang bream Megalobrama amblycephala, Mozambique tilapia Oreochromis mossambicus, Chinese catfish Hemibagrus macropterus, warmouth Lepomis gulosus, and Chinese loach Sinibotia superciliaris (Deng et al., 2020). For example, invasive species like Mozambique cichlid O. mossambicus occurred at site 1, site 3, and site 10, and meantime grass carp C. idella at site 10. Mozambique cichlid O. mossambicus is a successful and vagile invader. Grass carp C. idella, which originated from the Yangtze River in China, was a pest because of its damage to submerged vegetation (Deng et al., 2020). Introduction of these species possibly excludes endemic species with similar niches. For example, endemic species mud carp Cirrhinus molitorella, a near-threatened fish on IUCN red list (IUCN, 2020), was collected at sites 1 and 10 in 1979 but was absent in 2015. This species grazes on algae, phytoplankton, and detritus, which is a similar dietary with Mozambique cichlid O. mossambicus (Froese & Pauly, 2019). The importance of species invasion revealed by RF models indicates that these negative biotic interactions also contribute to community structure in the Min River.

Human population was the least important factor which was explanatory for FRic, FDsor, FDsne, and FDsim. Many anthropogenic activities associate with human population. Urbanization stimulates natural resource consumption and freshwater supplies for cities. Increased demand for water supply indicates more water reservoirs to build, and water pollution is aggravated in areas with high human population density (De Mello et al., 2020). People’s demand for food is also influential. For example, citizens prefer certain fishes such as yellow catfish Tachysurus fulvidraco in local cuisine (He et al., 2012). Human population along the Min River had increased by 36% from 1979 to 2015, which indicates strengthened temporal human influences. Nevertheless, ranking the relative importance of different factors confirms that environmental filtering is the dominant mechanism driving fish community assembly in the Min River temporally.

Responses of functional traits to environmental variables

Ordination of environmental factors and fish functional traits based on fourth-corner combined with RLQ supports previous findings that precipitation, temperature, and species invasion were the most important factors driving fish functional diversity. Meantime, sites in 2015 were more disaggregated than in 1979 on the ordination diagram, suggesting that trait compositions in 2015 were more varied than in 1979. This supports the previous finding that functional beta diversity in 2015 was significantly higher than that in 1979.

Higher swimming factor and higher relative eye size were significantly correlated with higher precipitation, higher temperature, and presence of invasive species. Swimming factor is directly associated with fish swimming capability. Higher values of swimming factor indicate relatively higher water propulsion per swing of the caudal fin (Villeger, Grenouillet & Brosse, 2013). This possibly requires that fish species living in habitats with higher precipitation and temperature adapt to fast-flowing water. For example, endemic species Pseudobagrus tenuis and Pseudolaubuca sinensis with higher swimming factor and relative eye diameter both occurred at sites 3 and 10 with a relatively higher average temperature of 19 °C and annual precipitation of about 2,000 mm. These two sites were also inhabited by invasive omnivorous species Mozambique cichlid O. mossambicus, which possibly incurred intensified interspecific competition for food with P. tenuis and P. sinensis. Meanwhile, relative eye diameter is correlated with the ability for food detection (Beston et al., 2019). Higher swimming capability and relative eye diameter of these two species assure that they can secure adequate food resources successfully in an ecosystem with fast-flowing water and intense interspecific competition (Starrs et al., 2017).

Eurytopic water flow preference and higher relative pectoral fin length were significantly correlated with absence of invasive species, lower precipitation, and lower temperature. Pectoral fin length is assumed to be correlated with the capability of low-speed maneuvering of fishes (Jones et al., 2020; Watson & Balon, 1984). For example, endemic species Gobiobotia longibarba has a higher relative pectoral fin length. It was collected at sites 12, 13, and 14 in 1979 but absent in 2015. Annual precipitation at these three sites was about 1,680 mm in 1979 but was 2040 mm in 2015. Lower precipitation indicates lower water velocity during the flood season. Endemic species like G. longibarba thus do not need to resist fast-flowing water nor to frequently swim fast to compete for food resources. On the other hand, higher relative pectoral fin length gives these fishes higher slow-speed maneuvering performance in relatively slow-flowing water, which is possibly important to detect food in local habitats (Jones et al., 2020).