Morphological variation of African snakehead (Parachanna obscura) populations along climate and habitat gradients in Côte d’Ivoire, West Africa

Intraspecific morphological variation may reflect phenotypic plasticity or adaptive divergence. While adaptive shape divergence may occur more likely among isolated populations with reduced gene flow, phenotypic plasticity may reflect morphological responses to heterogeneous environments, even in spatially connected populations. We evaluated both processes while examining morphological variations among seven wild populations of snakehead fish (Parachanna obscura) along climate and habitat gradients in Côte d’Ivoire, West Africa. Morphological variations were studied by multivariate canonical variate analysis (CVA) as based on geometric morphometrics of 15 fish body landmarks. Correlations between shape variations among populations and climate and habitat characteristics and between morphological and geographic distances were calculated. We found significant morphological variations among the seven populations. The variations in fish shape were concentrated on landmarks related to swimming and feeding, suggesting a contribution of environmental variation to morphological differentiation. However, we did not detect significant effects of climate and habitat variables on fish shape. The trend between geographical and morphological distances was likewise not significant. Therefore, a mechanistic understanding of the factors causing shape variation among P. obscura populations in West Africa could not yet be achieved.


Introduction
In Africa, many researchers have reported that freshwater organisms, especially fish populations, are strongly affected by numerous threats including climate change and anthropogenic perturbations such as dam construction, water extraction, habitat modification, water pollution, and overexploitation (Darwall et al. 2011). For example, in Côte d'Ivoire, the climate has significantly changed over the twentieth century, and the trend is expected to continue in the future. Experts on The Intergovernmental Panel on Climate Change predicted an increase in temperature by 3 °C by 2100, and a decrease of precipitation by 8% during the next hundred years, under the RCP4.5 scenario in Côte d'Ivoire (Bernard 2014). It is predicted that by the 2050s, hydrologic conditions that are significantly different from the existing conditions will be experienced by more than 80% of Africa's freshwater fish species (Darwall et al. 2011), which is likely to have a significant impact on fish diversity.
One of the African freshwater fish species with great economic value is Parachanna obscura (family Channidae), commonly known as African snakehead fish. It is widespread in the intertropical convergence zone in Africa, where the water temperature ranges from 26 to 28 °C. Some studies have indicated a decline in the populations of African snakehead fish, mainly due to overfishing and other anthropogenic threats coupled with the effect of the ongoing climate change (Lalèyè 2020;Amoutchi et al. 2021). However, Channidae species, including P. obscura, can withstand stressful conditions, for example, low oxygen concentrations, notably by the possession of a cavity above the gill chamber, which functions as an accessory respiratory organ (Lee and Ng 1991;Herborg et al. 2007;Kpogue et al. 2012). In addition to these specific organs acquired during phylogeny, variation in body morphology along climate and habitat gradients (Munday et al. 2008;Lema et al. 2019;Del Rio et al. 2019) may increase the resilience of species to the effects of global change.
Intraspecific morphological variation is often produced and maintained by divergent selective regimes, for example, in response to abiotic and biotic factors that differ among the habitats of populations. Under these conditions, morphological variation can reflect phenotypic plasticity of a common gene pool. However, morphological diversity may also reflect adaptive processes resulting in strong genetic divergence among populations (Lazzarotto et al. 2017), with the potential to enable populations to reach new adaptive peaks (Mouseau et al. 2000). While adaptive morphological variation based on genetic divergence develops and is maintained more likely among spatially isolated populations with reduced gene flow (Rundle and Nosil 2005), phenotypic plasticity is predicted to be most advantageous in a heterogeneous environment when future conditions can reliably be predicted (Sultan and Spencer 2002).
Cote d'Ivoire is a West African country characterized by a diversified climate and water network system. The climate divides the country along the latitude into three principal climatic zones such as Guinean in the south, Sudano-Guinean in the middle, and Sudanian in the northern region. In addition, there is a fourth zone being the particular climate of the mountain zone in the western region (Amoutchi et al. 2021). The Guinean zone is characterized by a sub-equatorial climate, the Sudano-Guinean zone by an equatorial transition climate, and the Sudanian zone is characterized by a tropical climate. The water systems of the country are characterized by a vast and complex system including many lakes, wetlands, rivers, and streams. Four major river basins with several tributaries, namely, Comoe, Sassandra, Cavally, and Bandama, and many further coastal rivers are found in the country ( Fig. 1) (Girard et al. 1970). Thus, if morphological variation among P. obscura populations is caused by responses to differing local environmental conditions and magnified by spatial isolation among populations, we would expect a strong morphological variation of this species within Cote d'Ivoire.
In our study, we used geometric morphometrics to examine morphological variation among seven wild populations of P. obscura along the climate and habitat gradients in Côte d'Ivoire. We tested whether shape variation correlates to differences of climate and habitat characteristics and examined whether geographical distance among populations correlates with morphological distance. We expected the variations to be localized primarily in fish body parts associated with ecological functions, for example, in fins and caudal peduncle involved in swimming and maneuvring, thus reflecting phenotypic plasticity. A strong correlation between morphological and spatial distance might be indicative of a stronger effect of spatial isolation on morphological divergence.

Field sampling and data collection
The study was conducted on several Ivorian inland watersheds (Fig. 1, Table 1). Bia River (KRIN), Sassandra River (SBR), and San-Pedro Lake (SANP) were sampled in the southern region (Guinean climate zone). Kan Lake (KAN), located centrally, and Buyo Lake (SAS) and Nzo River (NZO), located in the center-western part of Cote d'Ivoire, were sampled from the Sudano-Guinean zone. Finally, the Bagoue River (BAG) from the northern region was sampled in the Sudanean climate zone.
A total of 167 fish from seven P. obscura populations were collected (Table 1). P. obscura fishermen from each sampling site were employed for fish  (Table 1) from which P. obscura specimens were collected, and position of Côte d'Ivoire in West Africa sample collection using traps and fish nets. Within 1 h from capture, the left side of each specimen was photographed using Nikon (D5600) camera (calibrated against a graduated ruler to standardize measurements) for geometric morphometric analysis. The short phase between capturing and photographing fish avoided unwanted morphological distortions, for example, from rigor mortis. Specimens' length and weight were also measured. Then, specimens were fixed in 10% formaldehyde solution and transported to the laboratory for further analysis. Physico-chemical parameters of each site including water temperature (°C), pH, total dissolved solids (ppm), electrical conductivity (µS/cm), and redox potential (mV) were measured in situ with a handheld multi-parameter monitoring device (WQM-243/PHTK-243, TEKCOPLUS, Hong Kong, China). Water depth (cm) and water transparency (cm) were also measured (Appendix Table 4). The canopy cover of vegetation (%) along the shoreline was estimated visually. Measurements were done over the dry and rainy seasons at each site in September 2020 and January, June, and December 2021. The values were later averaged to characterize the situation of the locations during both seasons (Appendix Table 4 ). The latitude and longitude of each site were measured using a handheld global positioning system. Additional climate data (mean air temperature, maximum and minimum air temperature, relative humidity, and mean precipitation) were obtained for each sampling site from NASA POWER Data Access (https:// power. larc. nasa. gov/ data-access-viewer/) by entering the latitude and longitude of sampling sites (Appendix Table 4 ).

Morphological variation
To perform geometric morphometric analysis, landmarks were converted to shape coordinates by the generalized least square Procrustes superimposition (GPA; Rohlf and Slice 1990). This method transforms landmark configuration by preserving all information about shape differences among specimens but removes non-shape variability in raw coordinates such as location, orientation, scale, and rotation, and then standardizes each specimen to unit centroid size (CS). Centroid size represents a dimensionless size measure computed as the square root of the summed squared Euclidean distances from each landmark to the specimen centroid (Bravi et al. 2013). The shape data obtained were further size-corrected using a regression of shape (Procrustes coordinates) on size (log centroid size) for each population separately. The obtained regression residuals were used as standardized shape data in further analyses on shape variations (see for a similar approach Scharnweber et al. 2013). By using a set of 21 polymorphic microsatellite markers, we have recently demonstrated that the populations of snakehead fish in Côte d'Ivoire including those studied here were genetically differentiated (Amoutchi et al. 2023). Therefore, we applied canonical variate analysis (CVA), a multivariate statistical method used to distinguish variation among multiple pre-defined groups (populations) of specimens, to quantify inter-population shape variations. CVA results are statistically reported as Mahalanobis distance (Md), which represents a measure of distance among populations relative to the within-population variation. The statistical significance among Mahalanobis distances were determined by bootstrapping the Md matrix (permutation test with 10,000 runs). All the subsequent analyses were performed in MorphoJ version 1.07a software. Probabilities < 0.05 were considered significant. Group assignments were cross-validated by a jackknife resampling test, performed in R package Morpho (Schlager 2017).

Influence of geographical isolation or distance on morphological variations
To test for the effect of geographic distance on the morphological variation among P. obscura populations, the Mantel correlation between the pairwise Mahalanobis distances (Appendix Table 5) and the natural logarithm (ln) of the pairwise Euclidean geographic distances (km) was calculated. The geographical distance matrix was calculated by QGIS 3.22.6 using the global positioning system (GPS) coordinates of the sampling locations. Mantel correlation was run in vegan package (Oksanen et al. 2022) in R, with 1000 permutations.

Climate and habitat variability and its influence on P. obscura morphology
To assess the environmental heterogeneity among collection sites of P. obscura populations, principal component analysis (PCA) on all physico-chemical, climatic, and habitat variables (Appendix Table 5 ) was applied. The PCA was calculated by R.4.1.2 software using the factoextra' package (Kassambara and Mundt 2020). For testing the link between fish body shape and climate and habitat variables, Pearson correlations were calculated between the first two principal components explaining the highest percentage of variability in climate and habitat variables among sampling sites and the population centroids of body shape from CVA. Pearson correlation was computed in the stats package in R4.1.2 (R Core Team 2021).

Morphological variation among P. obscura populations
Average length and weight of specimens varied among populations (Table 1). Fishes of KAN and SBR were the largest, while NZO, KRIN, and SAS populations were composed of smaller specimens. However, the length range of specimens suggested that only adult individuals were included in the study.
Canonical variate analysis indicated significant (p < 0.05) Mahalanobis distances (Md) among all population pairs (Appendix Table 5 ), with the highest Md values observed between KRIN, KAN, and the other populations. The classification rates obtained from the CVA based on the jackknife cross-validation test resulted in an overall classification rate of 60% of the specimens into their original groups. The highest classification rates were recorded for NZO and SANP, while the lowest classification rate was recorded for BAG (Table 2), indicating high intra-group variability within this population. The first two CVA factors contributed 38.4% (CV1) and 21.9% (CV2) of the total variation. CV1 reflected variation in the caudal body part (landmarks 6, 7, 8, 9, and 10), dorsal fin length (landmark 5), and upper body depth (Fig. 3). CV2 described variation in the anterior region of the body, notably in the snout orientation (landmarks 1 and 2), lower body depth (landmark 11), the shape of the upper region of the head (landmark 4), eyes position (landmark 3), and the positions of the pectoral fin and posterior edge of opercula (landmarks 13, 14). Overall, however, the morphological variation between individuals and populations was relatively weak, and in no region of the fish body was particularly strongly accentuated.
Despite the statistically significant differentiation among populations, individuals of KAN and BAG were strongly scattered in the CVA plot, and hence, their population ellipses overlapped substantially with those of other populations (Fig. 4). Apart from these two populations, the other five populations were grouped distinctly from each other on the scatterplot. NZO and SANP were the most distinct groups, with only few individuals overlapping with other groups. Individuals of KAN and KRIN populations with negative CV1 and CV2 values had a deeper and slightly shorter caudal peduncle and dorsal fin, less deep body, a posterior edge of their operculum and pectoral fin in a downwardly sloping position, the snout slightly sloping upwards, and eyes positioned in the bottom part of the head towards the snout, in contrast to individuals from NZO. Individuals of SANP with negative CV2 and positive CV1 values were characterized by a slightly longer caudal peduncle narrowing at the insertion bases of the dorsal and anal fins, a slightly longer dorsal fin, a posterior edge of their operculum and pectoral fin in a downwardly sloping position, the snout slightly sloping upwards, and eyes positioned in the bottom part of the head towards the snout. These individuals were morphologically opposed to individuals of the SAS population. Finally, SBR individuals, with negative CV2 values and central CV1 values, had a posterior edge of their operculum and pectoral fin in a downwardly sloping position and the snout slightly sloping upwards with eyes positioned in the bottom part of the head towards the snout. Individuals of the BAG population were morphologically strongly diverse, and hence, their positions were scattered throughout the diagram.

Influence of geographical distance on morphological variations
There was no significant correlation between pairwise morphological distance and geographical distance among the seven populations (Mantel test, P = 0.23, Fig. 5).
Climatic and habitat characteristics of collection sites of P. obscura populations PCA of 13 climate and habitat variables extracted two principal components explaining 76.2% of the total variation (Fig. 6). PC1 (47.8% of total variation) mainly reflected the variation in climatic variables such as relative humidity (RH), maximum air temperature (T_max), mean air temperature (T_Mean), minimum air temperature (T_min) and mean precipitation (precipit), as well as the habitat variables like water temperature, redox potential (ORP), electrical conductivity (cond), and total dissolved solids (TDS). T_max, T_Mean, ORP, and water temperature were negatively correlated with PC1, while all other variables were positively correlated with PC1 (Table 3). PC2, with a loading value of 28.5%, mainly reflected the variation of habitat variables, including water depth (Depth), transparency (Transp), and canopy. Except canopy, which was positively correlated with PC2, all other variables were negatively correlated with PC2. The habitat of the KRIN population was characterized by a high canopy, low water transparency, acidic pH, and shallow water depth. The habitats of NZO, SAS, and SBR populations were characterized by high water transparency, neutral or basic pH, low water and air temperature, a short canopy, moderately high precipitation and humidity, and deeper water. The collection sites for BAG and KAN populations had higher mean and maximum air temperatures, higher water temperatures, high redox potential, and low precipitation and relative humidity. The SANP population collection site was   Table 1), and habitat and climate variables (arrows) on the first two principal components explaining 76.2% of the total variation. T_Mean, mean air temperature (°C); T_max, maximum air temperature (°C); T_min, mini-mum air temperature (°C); RH, relative humidity (%); precipit, mean precipitation (mm/day); T_water, water temperature (°C); Depth, water depth (cm); Cond, electrical conductivity (µS/cm); TDS, total dissolved solids (ppm); ORP, redox potential (mV); Transp, transparency (cm); Canopy, canopy (%) the most polluted one with the highest values for chemical parameters such as electrical conductivity and total dissolved solids. Precipitation and humidity were also higher at this site.

Correlation between morphological variation and climate and habitat variables
Overall, there was only a weak correspondence between morphological variation and environmental conditions among populations (Fig. 7). The highest correlation coefficient (r = 0.68, P = 0.09, n = 7, df = 5) was obtained between CV1 and PC1. With declining temperature, increasing precipitation, humidity, and conductivity, and at higher productivity (TDS) (high PC1), the caudal peduncle of snakehead fish was elongated and narrows at the insertion bases of the dorsal and anal fins, the dorsal fin becomes longer, the body becomes deeper in its anterior part, and the posterior edge of the operculum and pectoral fin shift towards a downwardly sloping position (high CV1). All other combinations of PC and CV axes were less strongly correlated.

Discussion
We analyzed intra-specific morphological variation of P. obscura populations along climate and habitat gradients in freshwater habitats of Côte d'Ivoire. Significant morphological variation was observed among all populations, although the within-population variability was also high, and was particularly strongly expressed in the population with the lowest number of studied individuals (KRIN). The morphological variations were principally documented in the head, caudal peduncle, and pectoral and pelvic fins. Other studies on fish species such as Gobiocypris rarus (Shao et al. 2007), Hemigrammus coeruleus (Lazzarotto et al. 2017), Nothobranchius orthonotus ( Bravi et al. 2013;Vrtílek and Reichard 2016), Channa punctatus (Kashyap et al. 2016), or other Cyprinidae family fish (Jacquemin and Pyron 2016) corroborate our findings that intra-specific morphological variation in fishes can be strongly significant. However, we did not find potential causes of intra-specific variability. The isolation by distance test showed no significant correlation between morphological variations and geographical distance. Furthermore, there was only a weak correspondence between morphological variation and differences of environmental conditions among populations. Accordingly, we did not find significant evidence for environmental-driven phenotypic plasticity or for the effects of spatial isolation on the expression of intraspecific morphological variation in this species. CVA results showed that most shape variations were localized in the head (shape, eye, and snout positions), caudal peduncle, and pectoral and pelvic fins. The study by Mouludi et al. (2019) on two populations of Channa gachua, another species belonging to the Channidae family, also showed body shape differences in the position of the snout and caudal peduncle length. Similar morphological variations mostly localized in fins and head region were observed among populations of Lavinia symmetricus (Brown et al. 1992). Variation in eye diameter and head length was observed among the population of Trachurus trachurus through linear morphometrics (Bektas and Belduz 2009). A study on Phycis phycis populations also showed that morphological variation was mainly found in fin regions as well as in the position of the posterior limit of the operculum (Vieira et al. 2016). These regions of fish body are mainly associated with Table 3 Loadings of the 13 climate and habitat variables on the first two principal components explaining 76.2% of total variation and with Eigenvalues greater than 1. Highest loading values are bolded T_Mean mean air temperature (°C), T_max maximum air temperature (°C), T_min minimum air temperature (°C), RH relative humidity (%), precipit mean precipitation (mm/day), T_ water water temperature (°C), Depth water depth (cm), Cond electrical conductivity (µS/cm), TDS total dissolved solids (ppm), ORP redox potential (mV), Transp transparency ( functions such as fish locomotion (pelvic and pectoral fins and caudal peduncle) and feeding (eyes and snout), vital functions that are essential for the survival of the species. Fin and caudal peduncle dimensions are related to swimming ability, which includes speed and maneuverability (Sampaio et al. 2013;Li et al. 2021). Thus, variation in these regions may help the fish coping with environmental variability between different locations. However, in our study, no significant associations were found between fish shape and climate and habitat variables. Shape variation of many species can respond to temperature, dissolved oxygen, radiation, water depth, food availability, current or hydrology, or correlated with their feeding mode, predation risk, and habitat use (Swain et al. 2005;Bunje et al. 2007;Spoljaric and Reimchen 2011;Antonucci et al. 2012;Bravi et al. 2013;Scharnweber et al. 2013;Jacquemin and Pyron 2016). Jacquemin and Pyron (2016) demonstrated that changes in the morphology of Cyprinid fishes can be related to hydrology, where more fusiform individuals with downward turned and extended caudal peduncle dominated during periods of lower discharge. Scharnweber et al. (2013) demonstrated that Rutilus rutilus cyprinids from lakes with high predation pressure were characterized by a more streamlined body and caudally inserted dorsal fin, attributes that facilitate escape from predators. Phenotypic plasticity of body morphology is regarded as a key strategy for populations living in changing environments, although morphological variation may even play a role in adaptive diversification (Svanbäck 2004). We did not find significant evidence of morphological variations in the African snakehead fish in response to variation in environmental conditions. However, there was a relatively high correlation coefficient obtained between CV1 and PC1 (r = 0.68), which supports that morphological variation between snakehead fish populations may in part be affected by environmental variation. Thus, further investigations including more populations may shed light on the factors that induce intra-specific morphological variation in this fish species.
Geographical isolation of populations can lead to morphometric variations between populations, and this morphometric variation can provide a basis for Fig. 7 Pearson correlation plots between the first two CVA axes (CV1 and CV2), representing morphological variation among populations, and the first two PCA axes (PC1 and PC2), which represent variation in habitat and climate varia-bles between sites. Upper left plot-CV1 and PC1, upper right plot-CV1 and PC2, lower left plot-CV2 and PC1, and lower right plot-CV2 and PC2. The correlation coefficient (r) and probability (p) are shown on each plot population differentiation (Bookstein 1991). Spatial isolation by unconnected river systems or by connectivity barriers within river systems may prevent gene flow among populations, and hence, adaptive morphological diversification may have contributed to shaping the morphological variation among populations. In contrast, connectivity of river systems facilitates the movement of organisms, nutrient transfer, and energy flow across the landscape (Opperman et al. 2010), and is used by fish for reproductive migrations, recruitment and recolonization, and access to floodplains (Pettit et al. 2017). Spatial connectivity also allows genetic connectivity among populations and would facilitate weak genetic divergence, such that morphological variation would reflect primarily phenotypic plasticity. In our study, no significant association was observed between spatial and morphological distances, despite the fact that we found significant morphological variations (as indicated by Mahalanobis distances) among the populations obtained from hydraulically unconnected sites. Furthermore, we also found differentiation among populations from the same basin (SBR and SAS, both from Sassandra River), likely caused by a dam constructed on Sassandra River. Indeed, the SBR population was collected from the downside of the Sassandra River downstream of the dam, while the SAS population was collected from the upstream side of the river in the artificial Lake Buyo which is the reservoir impounded by the dam. Finally, fishes of the NZO population that were collected from Nzo River, a tributary of Lake Buyo, were also morphologically different to SBR and SAS populations. A significant effect of spatial isolation on genetic differentiation of snakehead fish populations including those studied here (Amoutchi et al. 2023) supports that spatial arrangement of rivers and potential connectivity among fish populations can be important. Larger sample sizes and a spatially refined sampling scheme be needed to detect the spatial effects on the morphometry of the snakehead fish.

Conclusion
The current study demonstrates significant morphological variation among seven P. obscura populations from sites that strongly differed in climate and habitat variables. However, the correlations between fish shape and climate and habitat variables were weak, and spatial isolation likewise did not contribute substantially to the differences of fish shape among populations. The differences in shape as obtained by geometric morphometrics may reflect phenotypic variation that may aid the fish in surviving in a variety of environments and may also contribute to enhanced fitness under the various threats of freshwater ecosystems in Africa, such as climate change, gold mining, water withdrawal for human needs, overexploitation, industrial waste discharge, pesticide use for agricultural purposes along watersheds, and obnoxious fishing practices (Amoutchi et al. 2021). However, a functional and mechanistic understanding of the adaptive value of shape variations between populations would require common garden experiments and extended analyses of functional genomics, data that are currently not yet available.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.