The Relative Effects of Biotic and Abiotic Factors on the Recruitment of Freshwater Mussels ( Margaritifera laevis )

: Freshwater mussels, Unionoida, are endangered across the globe due to recruitment failure. In the present study, with general linear mixed models, we investigated the relative effects of biotic (host ﬁsh density) and abiotic (water depth, ﬁne sediment, water temperature, and water quality) factors on the recruitment of Margaritifera laevis in 10 streams of Hokkaido, northern Japan. We additionally examined the factors regulating the density of the host ﬁsh Oncorhynchus masou masou with general linear models. The proportion of juvenile mussels had a unimodal relationship with the host density, which was the most inﬂuential factor among the others examined. The positive relationship between mussel recruitment and host density can be attributed to an increased host ﬁsh infection rate. The negative correlation between mussel recruitment and host density at high ﬁsh densities may be due to reduced larval growth on host ﬁsh that are in poor physical condition. We also found that host ﬁsh density was negatively affected by nutrient enrichment. Our results suggest that mitigating water quality degradation to recover host ﬁsh density should be prioritized to improve mussel recruitment. Although stock enhancement is effective for increasing the salmon population density, excess stocking can further disturb mussel recruitment.


Introduction
Biodiversity has rapidly declined during the sixth mass extinction [1]. Clarifying the causes of species decline is an essential first step for evidence-based conservation. Previous research has focused primarily on population size as an indicator of current population status e.g., [2]. However, population size does not necessarily reflect actual extinction risks. For example, recruitment can be impaired even in large populations of long-lived species [3]. Such apparently healthy populations (i.e., large populations) of long-lived species will go extinct if recruitment does not recover in the future. Therefore, the causes of recruitment failure must be carefully examined to conserve populations of long-lived species [4,5].
In freshwater ecosystems, unionid mussels (order Unionoida) are among the most long-lived species (often reaching 100 years of age), and they play significant roles in purifying the water and providing habitats for other freshwater organisms [6,7]. As with other long-lived species, many species of freshwater mussels have undergone substantial declines across the globe [8][9][10]. For example, approximately 70% of Japanese unionid mussels are listed on the national Red List as either endangered or threatened species [11,12].
While the acute impacts of physical habitat destruction are well-recognized as a major cause of mussel decline, less understood are the chronic stressors that impair the recruitment processes of unionid mussels. Therefore, an in-depth understanding of the factors limiting mussel recruitment is of critical importance for conserving this highly imperilled group of animals.
Biotic and abiotic factors may be involved in chronic recruitment failure in unionid mussels. Freshwater mussels have a complex life cycle in which larvae called glochidia are obligate parasites of the gills or fins of host fish. After the larval parasitic life stage (which typically lasts for several weeks to months), juvenile mussels shed from the host fish and burrow into the substrate [13][14][15]. Therefore, the recruitment of freshwater mussels is influenced by biotic (i.e., host fish) and abiotic factors during their parasitic and benthic life stages. The host fish density in the surrounding environment may regulate the host fish infection rate (see, e.g., [16][17][18][19]). Moreover, excessive sedimentation leads to recruitment failure in freshwater mussels through the impairment of benthic environments [20][21][22]. There is increasing awareness that a combination of biotic and abiotic factors may determine the recruitment of freshwater mussels [16,21,23]. An approach that incorporates both biotic and abiotic factors may illuminate the relative effects of these factors, providing useful insights into the conservation of unionid mussels. However, the relative effect is rarely tested because biotic and abiotic factors have often been considered separately.
Examining the factors that regulate host fish density is also critically important for the conservation of unionid mussels. This is because habitat degradation indirectly influences the parasitic stage, which is crucial to the initial growth of juvenile mussels, by reducing the host fish density [24,25]. A combined investigation of the determinants of mussel recruitment and host fish density can improve our understanding of the conservation of unionid mussels.
Margaritifera laevis is a Japanese freshwater mussel whose glochidia are obligate parasites of the gills of juvenile masu salmon, Oncorhynchus masou masou. This long-lived species has an estimated lifespan of 80 years and is distributed in cold streams in Sakhalin, Russia, and the Honshu and Hokkaido Islands of Japan [26]. The population of M. laevis has experienced nationwide decline over the past few decades, and M. laevis is now listed as an endangered species on the Red List of Japan [12]. Although northern populations of M. laevis in Hokkaido, Japan seem to remain intact, the disproportionate lack of juvenile mussels implies recent recruitment failures [9]. In the present study, we examined the relative effects of biotic (host fish density) and abiotic (sedimentation and water quality) factors on the recruitment of freshwater mussels. In doing so, we statistically evaluated the relative effects of these potential factors on the proportion of juvenile mussels in the population (i.e., a proxy for juvenile recruitment). Furthermore, we explored the abiotic factors limiting the density of O. masou masou hosting M. laevis larvae. We expected that biotic factors (i.e., host fish density) had the most influence on the recruitment of M. laevis because released glochidia from female mussels cannot begin the juvenile stage if they cannot encounter the proper host fish individuals.

Quantitative Sampling and Species Detection
The present study was conducted in 10 streams in Hokkaido, Japan (eight streams in northern Hokkaido, one stream in central Hokkaido, and another stream in eastern Hokkaido) (Figure 1). Hokkaido is a large island (83,450 km 2 ) and one of the most famous agricultural prefectures in Japan. The proportion of farmland in the catchment area of each study stream was 0.25 ± 0.24 (measure ± SD) (Table S1). We sampled mussels in the summers of 2016 and 2017. We established a 2-10 km segment within each stream based on accessibility and mussel distribution. Within each segment, we established two or three sampling sites (width of the sites; 2.5-13.0 m) with 30-50 m reach lengths that included riffles and rapids. There were a total of 26 sampling sites. At each site, 15 quadrats were sampled with a Surber net (0.09 m 2 , 0.5 mm mesh size). Each quadrat was thoroughly investigated by excavating to a depth of 10 cm in the substrate and sieving the substrate immediately through a 2 mm mesh sieve. All collected mussels were photographed with a digital camera in the field, and the shell length and height were measured with image analysis software (ImageJ 1.50i) [27] in the laboratory. Margaritifera togakushiensis, which are taxonomically close to M. laevis, are found in Japan [28]. In our study, we distinguished M. laevis from M. togakushiensis by the method based on the ratio of the shell length to the shell height; this method has been shown to provide 80% accuracy for species identification in northern and central Hokkaido and 92% accuracy for species identification in eastern Hokkaido [29]. With this method, we eliminated the possibility of misidentifying M. laevis as M. togakushiensis to the extent possible.
Water 2021, 13, x FOR PEER REVIEW 3 of 11 riffles and rapids. There were a total of 26 sampling sites. At each site, 15 quadrats were sampled with a Surber net (0.09 m 2 , 0.5 mm mesh size). Each quadrat was thoroughly investigated by excavating to a depth of 10 cm in the substrate and sieving the substrate immediately through a 2 mm mesh sieve. All collected mussels were photographed with a digital camera in the field, and the shell length and height were measured with image analysis software (ImageJ 1.50i) [27] in the laboratory. Margaritifera togakushiensis, which are taxonomically close to M. laevis, are found in Japan [28]. In our study, we distinguished M. laevis from M. togakushiensis by the method based on the ratio of the shell length to the shell height; this method has been shown to provide 80% accuracy for species identification in northern and central Hokkaido and 92% accuracy for species identification in eastern Hokkaido [29]. With this method, we eliminated the possibility of misidentifying M. laevis as M. togakushiensis to the extent possible. Locations of the study streams where quantitative sampling was performed (gray circles) and the sampling areas where mussels were collected for age estimation (circles) in Hokkaido. The name and specific locations of the study streams are withheld to protect the habitat of M. laevis.

Age Estimation
To assess the recruitment status of the study populations, we estimated the proportion of juvenile mussels (the number of juvenile mussels divided by the total number of mussels) at each study site. In this study, mussels under 10 years old, which are considered to be sexually immature, were treated as juvenile mussels [30]. Juvenile mussels were classified based on the following growth models. The age of mussels can be estimated by counting the number of growth rings engraved onto the shell ligament. However, this method requires the sacrifice of live mussels and is therefore inappropriate for use with Figure 1. Locations of the study streams where quantitative sampling was performed (gray circles) and the sampling areas where mussels were collected for age estimation (circles) in Hokkaido. The name and specific locations of the study streams are withheld to protect the habitat of M. laevis.

Age Estimation
To assess the recruitment status of the study populations, we estimated the proportion of juvenile mussels (the number of juvenile mussels divided by the total number of mussels) at each study site. In this study, mussels under 10 years old, which are considered to be sexually immature, were treated as juvenile mussels [30]. Juvenile mussels were classified based on the following growth models. The age of mussels can be estimated by counting the number of growth rings engraved onto the shell ligament. However, this method requires the sacrifice of live mussels and is therefore inappropriate for use with endangered species [31]. Therefore, we modeled a shell length-age relationship to minimize the number of sacrificed individuals. We developed separate regional shell length-age relationships for northern Hokkaido, central Hokkaido, and eastern Hokkaido given the geographical distance between the study streams ( Figure 1). We collected 150 mussels from the study streams (50 mussels in northern Hokkaido, 28 mussels in central Hokkaido, and 72 mussels in eastern Hokkaido) and measured their shell length with image analysis software (ImageJ 1.50i) [27]. We estimated the age of the mussels by counting the number of growth rings under a dissecting microscope. Since the top of the ligament was often damaged, it was often not possible to read the exact age. Therefore, the growth rings of the damaged portion were estimated using the length of the ligaments and the number of growth rings of undamaged individuals [32]. We modeled the shell length-age relationships in the three regions with the following four nonlinear growth models (see, e.g., [31,33,34]): von Bertalanffy growth curve : Gompertz function : Hyperbolic saturation function : Logistic function : where L t is the shell length (mm) of a mussel at age t, L ∞ is the theoretical maximum length (mm), k is a growth coefficient, t 0 is the theoretical age at zero length, and a is a constant. It must be noted that the classical von Bertalanffy growth model used in previous studies with freshwater mussels is modified here, following Urban [34], and is not the "generalized" von Bertalanffy growth model. The performance of each model was assessed on the basis of their sum of squared residuals (RSS) values, and the model with the minimum RSS value for each region was applied.

Abiotic Factors
We measured five abiotic factors: water depth, proportion of fine sediment, summer maximum water temperature, total nitrogen concentration (TN), and total phosphorous concentration (TP). The water depth and proportion of fine sediment were measured at each study site, whereas the summer maximum water temperature, TN, and TP were measured at a single location within each study segment.
We measured the water depth (cm) at the corners and center of each quadrat. The mean water depth was calculated as the average of all the measurements at each study site. We took three sediment samples (right bank side, left bank side, centerline of the site; 300 mL) from the substrate with a trowel and sieved them with 2 mm and 0.042 mm mesh sieves to identify sediments fine enough for a mussel to burrow in. The sieved sediments that contained <2 mm and >0.042 mm particles were brought back to the laboratory in a cooler box with ice packs and stored in a freezer. In the laboratory, we sieved these sediments with a 0.25 mm mesh sieve (i.e., fine enough to impede interstitial water flow) to classify the sediments into fine sediments (<0.25 mm) that went through the sieve and coarse sediments (>0.25 mm) that remained on the sieve. After drying the sediments at 60 • C for 24 h, we weighed and calculated the proportion of fine sediment (fine sediment/fine + coarse sediments). The mean proportion of fine sediment was calculated as the average for the three samples at each study site. The maximum summer water temperature (in • C) was measured using data loggers (CO-U20L-04, Onset Computer Corporation, U.S.) located in each study segment. We defined the maximum summer water temperature as the maximum water temperature during July and August of 2018 based on hourly records of water temperature. To measure TN and TP (10 −2 g L −1 ), we collected surface water samples from each study segment in autumn 2017 and filtered them immediately through precombusted glass-microfiber filters (Whatman GF/F, GE Healthcare UK Ltd., Buckinghamshire, UK) using syringes. The filtered water samples were quickly brought back to the laboratory in a cooler box with ice packs and stored in a freezer at −15 • C until further analyses. The TN and TP of each sample were analyzed by the colorimetric method using a continuous flow injection analyzer (AACS-4, BL-Tech Co. Ltd., Osaka, Japan) after NaOH-K 2 S 2 O 8 digestion. The index of nutrient enrichment was calculated by principal component analysis of TN and TP because these values showed a strong correlation with each other (r > 0.7).

Biotic Factors
We assembled abundance data for the host fish, O. masou masou, across the 10 study streams. Only 1-year-old and younger fish were observed in the streams. In one or two survey reaches (survey reach length; 40.6 ± 23.2 m (measure ± SD)) within 2 km of each mussel survey study segment, one-pass host fish sampling was conducted with a combination of cast-net sampling and electrofishing from 2005 to 2016. We used the cast-net method as an efficient method to collect swimming fishes in riffles, which are difficult to collect with electrofishing. In the analysis, we used data from five or six non-consecutive years during the sampling period that were one or two years apart (Table S2). We calculated the mean host fish density during the sampling period and used that value in the following statistical analyses. The data were provided by the Hokkaido Regional Development Bureau and the Salmon and Freshwater Fisheries Research Institute.

Statistical Analysis
First, we developed general linear mixed models (GLMMs; n = 26) with a binomial error distribution via the R package "lme4" [35] to examine the relative effects of biotic and abiotic factors on the recruitment success of M. laevis. We used the proportion of juvenile mussels as a response variable and the biotic (logarithmically transformed host fish density) and abiotic factors (water depth, proportion of fine sediment, maximum summer water temperature) as explanatory variables. Since the host fish density and the index of nutrient enrichment showed a strong correlation (r < −0.7), we did not add these variables to the model simultaneously to account for multicollinearity. Thus, we used the host fish density as the explanatory variable after a preliminary analysis that compared the Akaike information criterion (AIC) values of the four GLMMs developed with each variable (Table S3). All explanatory variables were normalized (i.e., scaled to a mean of 0 and a standard deviation of 1). Stream identity was included as a random effect to account for the effect of unmodeled factors attributable to differences among the streams. To identify plausible models, we selected the best model with the lowest AIC via the R package "MuMIn" [36]. We performed model averaging when multiple best models (∆AIC < 2) were found. We interpreted the significance of the parameters in the best model or the averaged model on the basis of a 95% confidence interval that did not include 0. We evaluated the relative effect of each factor based on estimates of the standardized partial regression coefficients of the explanatory variables. Because the standard partial regression coefficient does not depend on the unit of each variable, it was possible to compare the degree to which the explanatory variable affected the objective variable.
Second, we used general linear models (GLMs; n = 10) with a Gaussian error distribution to examine the abiotic factors affecting the host fish density. We used the host fish density as a response variable and the four abiotic factors (water depth, proportion of fine sediment, summer maximum water temperature and index of nutrient enrichment) as explanatory variables. The same procedures used in the analyses of the recruitment success of M. laevis were adopted for the model selection. The model selection was based on the AIC for a small sample situation (AICc). All analyses were performed in R software v. 3.6.3 [37].

Results
The best-fitting models for estimating the shell length-age relationship were the Gompertz function in northern and central Hokkaido and the hyperbolic saturation function in eastern Hokkaido ( Figure S1, Table S4). The estimated shell length at 10 years of age was 49.76 mm (in northern Hokkaido), 17.00 mm (in central Hokkaido), and 28.56 mm (in eastern Hokkaido).
The number of individuals and the proportion of juvenile mussels varied among streams ( Figure 2). Among the 2279 mussels collected (shell length: 8.5-148.7 mm), 388 mussels were estimated as juveniles based on the shell length-age relationships. The proportion of juvenile mussels was 0.15 ± 0.17 (measure ± SD), and no juveniles were found at five sites in the three streams. The mean water depth and the proportion of fine sediment of the studied sites were 29.8 ± 10.2 cm and 0.18 ± 0.11, respectively, and the maximum summer water temperature, TN, TP, and host fish density of the studied segment were 21.5 ± 2.79 • C, 0.41 ± 0.34 10 −2 g L −1 , 0.02 ± 0.01 10 −2 g L −1 , and 0.24 ± 0.19 ind. m −2 , respectively (Table S1). The first principal component explained variations in the observed TN and TP values, which had a strong positive correlation with both variables (r > 0.9).
The number of individuals and the proportion of juvenile mussels varied among streams (Figure 2). Among the 2279 mussels collected (shell length: 8.5-148.7 mm), 388 mussels were estimated as juveniles based on the shell length-age relationships. The proportion of juvenile mussels was 0.15 ± 0.17 (measure ± SD), and no juveniles were found at five sites in the three streams. The mean water depth and the proportion of fine sediment of the studied sites were 29.8 ± 10.2 cm and 0.18 ± 0.11, respectively, and the maximum summer water temperature, TN, TP, and host fish density of the studied segment were 21.5 ± 2.79 °C, 0.41 ± 0.34 10 −2 g l −1 , 0.02 ± 0.01 10 −2 g l −1 , and 0.24 ± 0.19 ind. m −2 , respectively (Table S1). The first principal component explained variations in the observed TN and TP values, which had a strong positive correlation with both variables (r > 0.9).
Host fish density was the most influential factor of all the potential influencing ones. The model selection process revealed that the proportion of juvenile mussels was affected only by the host fish density (Table 1, Marginal R 2 : 0.246, Conditional R 2 : 0.250). The host fish density had a unimodal relationship with the proportion of juvenile mussels ( Figure  3a). The other tested factors had no relationship with the proportion of juvenile mussels (Figure 3b-d). The host fish density was affected only by the index of nutrient enrichment (Table S5, R 2 : 0.540), and this relationship was negative (Figure 3e).  Host fish density was the most influential factor of all the potential influencing ones. The model selection process revealed that the proportion of juvenile mussels was affected only by the host fish density (Table 1, Marginal R 2 : 0.246, Conditional R 2 : 0.250). The host fish density had a unimodal relationship with the proportion of juvenile mussels (Figure 3a). The other tested factors had no relationship with the proportion of juvenile mussels (Figure 3b-d). The host fish density was affected only by the index of nutrient enrichment (Table S5, R 2 : 0.540), and this relationship was negative (Figure 3e).  Figure 3. Relationships between the proportions of juvenile mussels and investigated factors: host fish density, water depth, proportion of fine sediment and summer maximum water temperature (a-d), and the relationship between host fish density and index of nutrient enrichment (d). The points represent the observed values. The solid line represents the predicted values (Table 1, Table S5).

Discussion
In the present study, we examined the relative effects of biotic and abiotic factors on the recruitment success of endangered M. laevis and established that host fish density was the most influential factor. In addition, we clarified the factors regulating the host fish density, which is closely related to the mussels' life cycle. Freshwater mussels are among the most threatened groups of animals in freshwater ecosystems [38,39]. Improving recruitment will be the key for the conservation of this taxon. Our findings provide critical information for the prioritization of conservation measures for target freshwater mussel species.
The results revealed that the host fish density, which showed a unimodal relationship with the proportion of M. laevis juveniles, was the most important factor influencing the recruitment of M. laevis. At low to moderate host fish densities, the proportion of juveniles increased with host fish density. Such a positive relationship was also observed in previous studies (see, e.g., [17]) and may be attributable to the high infection rate at moderate host densities [16]. However, the proportion of juvenile mussels decreased from moderate to the highest host fish densities. One possible reason for this negative correlation is decreased larval growth in host fish that are in poor physical condition. Among stream-dwelling Salmo species, density-dependent competition has been well-documented. Previous studies showed that increasing fish density caused a decline in individual growth rates or body mass (see, e.g., [40,41]). Hasegawa et al. [42] clarified that intraspecific competition reduces the stomach contents and growth rate of our target species, O. masou masou. Decreased larval growth in M. margaritifera was also confirmed when the host fish could not grow sufficiently [43]. Hence, we hypothesize that poor host fish conditions were induced by high host fish density, which finally deteriorated mussel recruitment by reducing larval growth. Elucidating the details of this mechanism will require further research.
Abiotic factors had no direct impact on mussel recruitment. However, we should note that physicochemical environmental factors may indirectly affect recruitment by decreasing the host fish density (see, e.g., [24]). In this study, nutrient enrichment was found to have a negative effect on the host fish O. masou masou. Similar patterns were observed in Ishiyama et al. [44], in which the population of O. masou declined with increasing nutrient enrichment in the Hokkaido stream. In that study, using data on Ephemeroptera, Plecoptera, and Trichoptera (EPT) abundance, the authors suggested that the primary cause of the salmon population decline was the decline in the abundance of preferred food resources.
Understanding the distribution and age composition of populations is a key issue for the conservation of M. laevis. However, there have been very few studies that quantitatively investigated juveniles of M. laevis in Japan (but see [9,45]). Our quantitative study has a meaningful outcome and provides the following management implications for improving the health of M. laevis populations. We found that conserving host fish density is the most important measure for the recruitment of this species. In the study streams, water quality improvements would help to maintain a high host fish density because the index of nutrient enrichment negatively affected the host fish density. In many regions of Hokkaido, catchments have been converted to farmlands, which have decreased freshwater biodiversity (see, e.g., [46,47]). A positive relationship between the index of nutrient enrichment and the proportion of farmland in the basin was shown in our study streams ( Figure S2). Therefore, excessive agricultural land use should be avoided to improve water quality and conserve host fish populations. Restoring riparian forests is an effective measure for improving water quality and promoting fish conservation because forests around streams filter nutrients supplied from catchments (see, e.g., [48,49]). On the other hand, the negative effect of excessive increases in host fish density must also be considered. The host fish, O. masou masou, is one of the key fishery resources in northern Japan, but excessive O. masou stocking programs may increase the risk of M. laevis recruitment impairment.
Few conservation measures targeting M. laevis have been implemented, and more efforts to promote the recruitment success of this species, as well as ex-post-evaluations of these efforts are needed [11]. Our findings regarding high-priority measures will contribute to the future management of the recruitment success of this species. However, none of the potential factors were examined in the present study. For example, the flow regime can greatly affect the recruitment success of M. laevis, because higher flows may cause bedload movement that is detrimental to juvenile mussels [50]. Further studies are needed to gain a deeper understanding of conservation measures. Moreover, we should note that it will take a long time (i.e., several decades) to accurately evaluate conservation outcomes for this species because the mortality rate of benthic juvenile mussels is high and because the shells of young juveniles are too small to observe for the first several years after they start to inhabit streambeds [30]. We believe that elucidating the relative effects of biotic and abiotic factors on mussel recruitment could be effective for understanding the ecology of and conservation priorities for other endangered freshwater mussels.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/w13091289/s1, Figure S1: Relationships between shell length and age in three regions, Figure S2: Relationship between the index of nutrient enrichment and proportion of farmland, Table S1: Abiotic and biotic factors used in the statistical analysis, Table S2: Host fish density in the study streams (A-J) from 2005 to 2016, Table S3: Results of preliminary analysis (GLMMs) for selecting explanatory variables from the host fish density and the index of nutrient enrichment, Table S4: Results of the estimations of four nonlinear growth models in the three regions, Table S5: Results from the GLMs examining the abiotic factors affecting the host fish density in ascending order of AICc.

Data Availability Statement:
The data that support the findings of this study are available from the Supplementary Materials or corresponding author upon reasonable request.