Study area and mussel beds
We conducted this study in the middle course of the Caeté River, an alluvial lowland river ~ 150 km long, located in northeastern Pará state, in the eastern Brazilian Amazon (Fig. 1). There is marked seasonality in river hydrology, which is a feature of Amazonian rivers (Junk, 1997), with an average discharge (± SD) of 48.3 ± 11.5 m3/s in the rainy season and 8.4 ± 2.9 m3/s in the dry season (Simeone et al., 2018). Formed by sandy alluvial sediments and a low channel slope, the Caeté has a predominantly meandering morphology, with irregular and tortuous meanders in the middle course and the formation of central sandbanks in the thalweg. The landscape consists of secondary forest floodplain, with small-scattered human settlements, subsisting mainly by fishing and family-based farming. Since large human settlements are scarce along the Caeté River, fluvial habitats are in relatively natural conditions, and there have been no artificial modifications to the channel along the course of the river.
In the study area, ~ 20 km long, we selected five sites known to have mussel beds (Fig. 1c) (Simeone et al., 2018, 2021a). Castalia ambigua is distributed in outer meander margins, associated with medium sand (ranging from 0.3 to 0.5 mm), low substrate heterogeneity, and stable riverbed in periods of high and low flows: Reynolds number between 500 and 1,800 (Simeone et al., 2021b). In addition, a hydraulic gradient was observed in the study area, with hydrodynamics near the riverbed increasing from upstream to downstream.
Classification of the river hydrological gradient
Reynolds number, which describes flow conditions near the riverbed, was associated with higher densities of mussels in the Caeté River (Simeone et al., 2021b). Thus, we used this variable as a baseline for an objective description of the hydrological gradient and to test our hypotheses. We used the Fisher-Jenks algorithm with the classInt package (Bivand, 2020) in GNU R 4.0.1 (R Core Team, 2020) to classify the hydrological gradient into zones. The Fisher-Jenks algorithm pooled the five sites into three distinct zones (Fig. 1c; Table 1) among which C. ambigua shell shape of both morphotypes were tested for differences. The zones were classified as follows: Zone 1 (sites 1 and 2 with low hydraulic energy); Zone 2 (corresponding to site 3 with medium hydraulic energy); and Zone 3 (sites 4 and 5 with high hydraulic energy).
Table 1 Mean ± SD values of the Reynolds number and C. ambigua shell length (mm) for each of the three hydrological zones classified using the Reynolds number on the basis of the Fisher-Jenks algorithm in the Caeté River, Bragança, Pará, Brazil
Hydrological zones
|
Reynolds number ± SD
|
Mean shell length (mm) ± SD
|
|
|
Morphotype I - Males
|
Morphotype II - Females
|
Zone 1
|
855 ± 94.6
|
36.3 ± 2.3
|
35.7 ± 1.9
|
Zone 2
|
1063 ± 68.3
|
39.2 ± 2.2
|
38.1 ± 1.5
|
Zone 3
|
1604 ± 88.9
|
43.9 ± 1.7
|
42.4 ± 2.1
|
Castalia ambigua shape acquisition and outline extraction
A total of 120 individuals of C. ambigua (60 individuals of each morphotype) were used in the shell shape analysis (Table 1). We obtained 75 individuals between November and December 2019 along a 20 m transect placed at each site using semi-quantitative searches, which provide better spatial coverage on all available hydraulic microhabitats. Mussels already collected from the same sites in November 2015 and conserved in alcohol in the laboratory (n = 45) were also included in the shape analysis. Juveniles were not included since age may mask habitat influence on mussel shell shape (Guarneri et al., 2014). Specimen size ranged from 34.7 to 44.5 mm along the maximum anterior-posterior dimension. As mussel shells may suffer damage due to predator activity and erosion caused by burrowing (Mata et al., 2019), only shells that were free from damage were used in the analysis. High resolution photographs of each individual were taken (Fig. 1d) using a digital camera attached to a table at a standardized height of 20 cm and a reference scale of 1 cm, with the shells in two positions: external view of the right valve (lateral view) and umbo upwards (umbonal view). Standardized imaging equipment is important in morphometric studies to reduce alignment error and analytical misinterpretations (Evin et al., 2020). Since mussel populations are declining worldwide (Böhm et al., 2020) and in the Caeté River (Simeone et al., 2021c), we took photographs of mussels between November and December 2019 in the field, returning the individuals to their original habitat in order to avoid unnecessary transport, stress and potential mortality. In the laboratory, we used the same standardized imaging equipment to photograph the mussels conserved in alcohol.
Photographs of both the lateral and umbonal views were processed with the GNU Image Manipulation Program (GIMP: https://www.gimp.org/) software and converted into black masks on a white background (8-bit gray-scale images). Images were then aligned on the mussel shell and stored in the same folder. Outlines were directly imported from each image and converted into a list of xy-coordinates using the function import_jpeg() in the Momocs package (Bonhomme et al., 2014) with GNU R 4.0.1 (R Core Team, 2020).
Sex determination of C. ambigua
For sex determination, we randomly selected 30 individuals of C. ambigua among those conserved in alcohol in the laboratory and that were used in shell shape analyses (15 from both morphotypes). We sexed each individual and examined the gonadal fluid under a light microscope. Mussels presenting oocytes or larvae (glochidia) in any developmental stage were assigned as females. Specimens with the occurrence of spermatozoa or sperm-morulae were indicated as males (Heard, 1975). The sex ratio was statistically compared with a chi-squared test in GNU R 4.0.1 (R Core Team, 2020).
Castalia ambigua shell shape analysis
We performed all shape analyses with GNU R 4.0.1 (R Core Team, 2020). We carried out an elliptic Fourier transformation of outlines (Haines & Crampton, 2000) using the efourier() function in the Momocs package (Bonhomme et al., 2014), to examine shell shape variation between both C. ambigua morphotypes among the hydrological zones. Fourier transformation decomposes an outline into a number of basic waves, called harmonics, which are powerful enough to extract the geometric information of these outlines, reducing and eliminating redundant information between adjacent xy-coordinates (Bonhomme et al., 2014). In addition, complex shapes can be fitted, outlines smoothed and rotational differences removed (Claude, 2008). After preliminary calibration, we reconstructed the outlines by increasing the number of harmonics (Crampton, 1995), and chose eight harmonics that corresponded to almost of 95% of the total harmonic power. Four coefficients per harmonic, which describe the size, shape and orientation of each harmonic (Haines & Crampton, 2000), were extracted for each shell outline and used as geometric variables. Because Fourier coefficients contain no size information, a standardization of the outlines to the same size is not necessary.
We performed a principal components analysis (PCA) on the matrix of coefficients in order to explore the morphological variability between both C. ambigua morphotypes among the hydrological zones and to define variables capturing most of the morphological variation (Bookstein, 1997). PCA rotates the original variables into new linear variables called principal components. The calculated principal components, which best explain the observed variation in shape, were considered as new shape variables.
We carried out a canonical variates analysis (CVA) based on the new shape variables in the Morpho package (Schlager, 2017), to identify the linear combination of shape features that was best able to account for differences between- and within-groups that were defined a priori (zones and morphotypes). We formally tested for differences among the hydrological zones and between C. ambigua morphotypes with a permutational multivariate analysis of variance (PERMANOVA) with a Euclidian distance matrix and 9,999 permutations. We used morphotypes and hydrological zones as fixed effects.
To understand the contribution of individual variables on the shell shape of C. ambigua, we reconstructed extreme outlines along each principal component (Crampton, 1995). To visualize differences at the extremes of the morphospace between both C. ambigua morphotypes, we generated deformation grids (Claude, 2008) and iso-deformation lines through mathematical formalization in the Momocs package (Bonhomme et al., 2014) using thin plate splines (TPS) analysis, which compares average shapes and visualizes the outline deformations in the morphospace (Bookstein, 1997).
Finally, we used a generalized additive model (GAM) with a cubic regression spline to explain C. ambigua shape variance in response to hydraulic energy: in our case the Reynolds number. A GAM is a generalized linear model with a linear predictor that allows for a flexible response on the covariates by defining regression splines (Wood, 2006). Initial data exploration revealed heterogeneous shape variance in the principal components. Thus, we standardized the shape data prior to analysis since we were not interested in between group (zone or morphotype) heterogeneity (Zuur et al., 2009). The response variable (Reynolds number) did not require any transformation. We ran the model using the mgcv (Wood, 2006) and nlme (Pinheiro et al., 2017) packages.