Phenotypic variation of Chitala chitala (Hamilton, 1822) from Indian rivers using truss network and geometric morphometrics

Collection of samples

A total of 149 dead specimens of C. chitala (a near threatened species) were collected from commercial riverine catches of seven rivers across India, viz, Son (n = 43), Tons (n = 16), Ken (n = 31), Brahmaputra (n = 16), Ganga (n = 13), Gomti (n = 20) and Gandak (n = 10), at a stretch from May, 2009 to November 2017 with repeated collections from a particular location to achieve a reasonable sample size (Table 1, Fig. S1). These samples form a subset of those studied by Dutta et al. (2020) through mitochondrial markers and had identified a single genetic stock from these rivers. No live animals were used in the experiment performed. All the fish, which were dead, were procured from commercial catches. The protocols followed were approved by the Institutional Animal Ethical Committee (IAEC), ICAR-NBFGR, Lucknow, India vide No. G/CPCSEA/IAEC/2015/2 dated 27 October 2015.

Digitization of samples

The specimens were first cleaned with water, drained, wiped dry and placed on a flat surface having laminated graph sheet on background, which was meant for calibration of the coordinates of digital images (Fig. 1). The methodology of capturing digital images for truss network analysis by Mandal et al. (2020) was followed. The fins were stretched for clear visibility of origin and insertion points. Each specimen, placed on graph paper, was labelled with a unique code and captured using DSC-W300 digital camera (Sony, Japan) with mouth facing towards left, which permits replication of the measurements (Cadrin & Friedland, 1999). Due to varying size range of samples as the images were captured from different heights, reference scale on the graph paper was used for calibration, using the software tpsUtil (Rohlf, 2008a).

Generation of morphometric data

Nine landmarks were identified that extended across the entire fish to represent the full dimensions of the morphology (Fig. 1) and truss data was generated from digitized images through tpsUtil V1.38 (Rohlf, 2008a) program for converting the images secured from JPEG (*.jpeg) to TPS (*.tps). The truss distances were extracted in tps files from the digitised images using tpsDig2 ver. 2.31 in X-Y coordinate data form (Rohlf, 2008b). PAleontological STatistics (PAST) software (Hammer, Harper & Ryan, 2001) was used to construct an inter-linked network of nine landmarks to generate 35 truss measurements.

Statistical analysis

The landmark-geometric morphometric (GM) analysis were performed on truss –morphometric data after eliminating size-effect. The data generated by PAST were log-transformed (Strauss, 1985) for linearization of data and to minimize size-related dissimilarity, all morphometric variables were adjusted by an allometric approach as proposed by Elliott, Haskard & Koslow (1995); M_adj =M*(L_S/L₀)^b, where M = original measurement, M_adj = size adjusted measurement, L₀ = standard length of fish, L_S = overall mean of standard length for all fish from all locations. As standard length (SL, character code 1_4) was used as a basis for transformation, it was excluded from the final analysis (Mamuris et al., 1998) and thus 35 morphometric variables were retained. One-way ANOVA for 35 morphometric characters was performed to evaluate the significant difference among seven rivers. The morphometric characters which showed high significant variations (P < 0.01) were used to obtain the recommended ratio of the number of samples (N) to parameters included (P), i.e., 3.5 to 8, for a stable outcome from multivariate analysis (Johnson, 1981; Kocovsky, Adams & Bronte, 2009). Therefore, 149 fish specimens of C. chitala, which is a subset of specimens collected for genetic stock structure analysis (Dutta et al., 2020); with nine-landmarks and 36 truss variables were found to be optimum for GM analysis.

For assessing distribution pattern of variance associated with variables principal component analysis (PCA) was applied on data matrix of 149 × 35 by reducing redundancy among the morphometric variables and eventually eliminating abundant independent variables (Samaee, Patzner & Mansour, 2009). Kaiser’s (1960) criterion of retaining eigenvalues greater than one (Jolliffe, 2002) was applied to determine PCA components. The canonical discriminant function analysis (CDFA) and discriminant analysis of principal components (DAPC) was also applied on truss data, to identify important truss variables associated with discriminant functions and the distribution pattern of specimens over geographical locations through discriminant functions. CDFA was applied on truss data for identification of important discriminant functions, important truss variables associated with discriminant functions and the distribution pattern of specimens over locations through discriminant functions (Cadrin, 2000). PCA identifies variations within specimens while CDFA identifies variations in specimens over locations. Therefore, a technique (DAPC) was also employed to simultaneously identify the within and between location variations. Further, to assess sensitivity of discriminant function in discrimination of specimen over different locations, the shape variation through receiver operating characteristics (ROC) curve analysis over (1-specificity) vs. sensitivity as axis with values varying in range from (0, 0) to (1, 1) was also carried out. In this analysis, area under curve of ROC performance graph of truss variables, identified by DAPC, indicated the role of discriminant functions in differentiating among locations, rivers, environments. ArcMap 10.8.1 platform (2020) was used to represent the coefficient of variation (CV) of significant top four truss variables, identified by DAPC, on geo-spatial scale. The river network shape file was prepared using the Google Earth platform and river basin shape file was obtained from NRSC India-Water Resource Information System (https://indiawris.gov.in/wris/). To understand geometric shape variations and distribution of landmark configurations, geometric shape variations in mean shape of specimens were analysed through thin-plate-spline (TPS) image analysis (Rohlf, 2008c). This was employed through procrustes superimposition method to extract shape information and infer shape changes, based on the anatomical landmark coordinates. Along with this, generalised procrustes shape variations over principal components of landmark coordinates were also performed to understand overall body shape variation across all landmarks simultaneously. Set of nine landmarks identified were partitioned into two subsets to check for shape variation through procrustes shape modularity.

All statistical analyses for truss-morphometric data were performed using R (ver.4.1.2) packages Adegenet 2.15 & ggplot2, SPSS version 16 (SPSS Inc. SPSS for Windows), SAS version 9.3 (SAS Inc: 2020), PAST (ver. 3.0), MorphoJ (Klingenberg, 2011) and Excel (Microsoft Corporation, 2018).

Principal component analysis (PCA)

PCA on 35 morphometric characters on (n = 149) individuals from seven locations, resulted in identifying 6 significant principal components that contributed up to 93.23% of total variation (Table 2) and eight truss variables with significant and highest loadings on PC-1 and PC-2 (Figs. 2, 3, Table 3). The two PCs cumulatively accounted for 56.46% of total variation with PC-1 accounting for 32.90% that led to the identification of two most important truss variables 5_9 (distance between point perpendicular to dorsal fin origin and point perpendicular to posterior end of maxilla), and 5_7 (distance between point perpendicular to dorsal fin origin and pectoral fin insertion), with significant loadings (Tables 2 and 3).

Relationship between truss and morphometric variables

Morphometric variation between significant truss measurements was identified after normalizing the data by transformation. As Elliot transformation accounts for only size related dissimilarity, an effort was made to compare variation using log transformed truss measurement too as it considers both shape and size dissimilarity. Comparison of the mean for truss variables in M-trans (Elliot transformation) ranged from 0.97–1.38 while in log-morphometric data it varied from 1.18–1.53 (Table 3). For instance, the truss variable, 5_9, in log-morphometric data exhibited higher standard deviation (0.15) and higher range value (0.68) than M-trans data (0.02, 0.15), with same value of mean (1.32). Thus, improved morphometric based identification using the truss variables in log-morphometric data is possible rather than M-trans due to higher standard deviation and range value indicating its greater utilization and application at field level (Table 3, Fig. 3).

Canonical discriminant function analysis (CDFA)

To identify important truss variables associated with discriminant functions, CDFA provided four canonical discriminant functions, which controlled 90.27% of variation with Wilks’ Lambda values ranging from 0.05–0.12, with significant chi-square value (p < 0.05). Two functions (1 & 2) out of the four, cumulatively explained 58.16% of total variation, where function-1 contributed the highest variation (33.06%) followed by function-2 (25.10%) (Fig. 4, Table 4). CDFA identified four most important truss variables with discrimination coefficient for 2_3 and 1_5 as 309.52 and 212.01 respectively for function-1 and for 3_9 and 4_8 as 582.49 and (-) 1034.82 for function-2, respectively. CDFA displayed significant differences among specimens from rivers Son, Ken and Ganga by distinct centroids, but similarity of specimen by merged centroids for river Tons with Brahmaputra and Gomti with Gandak (Fig. 4).

Discriminant analysis of principal components (DAPC)

To further check for distinctness and similarity between the rivers DAPC was employed that provided five significant canonical functions that control 96.04% variation (Tables S1, S2). The discriminant scores for each river (centroid) was observed as different and PCA eigenvalue in DAPC analysis (Table 5, Fig. 5) that represents variations between locations and DAPC analysis of DA & PCA eigenvalues provides seven distinct locations (Fig. 6). Probability score-based membership of 149 specimens among seven locations has higher probability score for correct assignment for specimens from Gandak followed by Gomti and lowest for Son (Fig. 7).

Comparison of performance of identified truss variables

The comparison of performances of identified eight truss variables by principal components (Table 3) and from DAPC (Table 6) pointed out improvement in coefficient of truss variable in DAPC analysis that resulted into seven distinct centroids for seven rivers, under present study. This was further verified from the variation in performance of coefficients in PCA, CDFA and DAPC of truss variables 5_9 and 5_7, (identified through PCA with highest loadings: Figs. SA2–SAC). Comparative study of the three graphs indicated that truss variable (5_9 and 5_7) identified through PCA had higher loadings, but had little role with small coefficients in functions from CDFA, while had higher role as coefficients in functions from DAPC analysis of data. Higher coefficients of these two significant truss variables, 5_9 and 5_7 identified maximum distinctness of each river. Performance analysis of DAPC coefficients identified 5_9, 5_7, 4_7 and 4_5, as top truss variables, which distinguished specimens from different rivers (Table 3). The variation in coefficient of variance of these four identified variables between the rivers are depicted in Fig. 8 (Table S3).

Receiver operating characteristics (ROC) curve analysis

The shape variation through ROC curve over (1-specificity) vs. sensitivity of discriminant function, analyzed for all seven rivers, showed that the value of area under curve (AUC) in each ROC ranged from lowest (0.9389) in Ganga to highest in Gandak (0.9971) (Fig. 9). AUC is a measure of sensitivity and as AUC values ranged between 0.90–1.00 for seven locations, it indicated an “outstanding performance” of identified truss variables in shape-based discrimination.

Functional distribution of truss variables over seven rivers

The functional distribution of the truss variables, which were identified from PCA and DAPC, when plotted to depict the shape variations over seven rivers, revealed that the distribution of PCA displayed seven graphs for seven rivers with similar position of peak (Fig. 10A), while from DAPC analysis, different and distinct distribution of specimens for seven rivers were observed (Fig. 10B). Thus, centroid of each location was dissimilar, though the observations were similar, indicating clear distinction and variation between locations (Fig. 10C).

Thin-plate spline and warp density in shape analysis

Geometrical shape variation between the locations was gaged by evaluating various parameters under thin plate spline shape analysis. Maximum and minimum warp score from PCA on principal components; PC-1 was highest for Ken River (5.42) and least for Ganga (2.32) (Table 7). Similarly, the maximum and minimum warp score on PC-2 was highest for Gomti (1.56) and least for Ganga (1.15). Thin plate spline variation over both mean score and principal score ratio was assessed. The geometric shape (thin plate spline) variation, over mean score, through PCA (relative warps) revealed that the specimens from Tons and Gandak rivers had higher relative warps (deformation) as 1.27 and 1.23, respectively, whereas other river basins had relative warps (deformation) in the range of 1.01–1.06 (Table 7, Fig. 11). The thin plate spline shape analysis through PCA (relative warps PC1, PC2), over principal score ratio (PSR) on PC-1 denoted as TPS-PSR-PC1, indicated that relative warp was observed highest for Gandak and least for Tons (0.87) (Table 7). However, TPS-PSR-PC2 distinguished the seven rivers with scores in the range from 0.89–1.11 (Table 7, Fig. S3). As specimens from Gandak showed high relative warp over principal score ratio, it implies that the geometrical shape of these are clearly distinct from other locations. But there was some minor similarity in geometrical shape between locations viz., Son and Brahmaputra on preliminary basis (Fig. 12). Thus, to capture small-scale deformations and variations between locations, warp density scores on principal components (PC1 and 2) was also assessed. The warp-density-score (WDS) for seven rivers over PC1 had values in range from (-) 50.50 to 21.40 while WDS values scores over PC2 ranged from (-) 33.67 to 20.20 (Table S4). The WDS over principal components (PC 1, 2) indicated the extent of shape variation of specimens collected from seven rivers (Table S4, Fig. 12). WDS indicated clear score for all locations indicating distinction between the specimens over locations shows shape-based variations.

Procrustes geometric morphometric analysis

The procrustes analysis of variance (ANOVA) of all specimens over the nine landmark-coordinates, to determine centroid size & shape variation (Tables S5, S6) indicated that there existed significant (p < 0.05) centroid size-based differences and significant procrustes shape variations for fish specimens from seven rivers. Highest procrustes distance (highest variation w.r.t. mean shape within location) was observed for Gandak (0.07), followed by Ken (0.02), Gomti (0.00), Son. (−0.01), Ganga (−0.02) and the least for Tons and Brahmaputra (−0.03) (Figs. 13 and 14) from generalized procrustes shape. This indicates that maximum distinct shape from generalized procrustes shape was observed in Gandak specimens while it was least for Tons and Brahmaputra.

Generalised procrustes shape variations over principal components analysis and between locations

To visualise the shape variation, landmark-coordinate data matrix was subjected to principal component analysis that identified seven principal components, which controlled 97.32% of total variance. The shape variation of 149 fish specimens has been observed on maximum & minimum scores on procrustes-principal components PC1 (PC2) (Table S7, S8, Fig. 15). The differences in mean shape over procrustes distances were significant for Tons & Brahmaputra (procrustes distance: 0.01903507, Mahalanobis distance: 3.9625, p = 0.0011) and between Gandak & Gomti (procrustes distance: 0.08225527, Mahalanobis distance: 5.2176, p = 0.003). The shape variations for location Tons & Brahmaputra was also observed over shape differences through lollipop graph, discriminant scores graph and cross-validation graph. Similar graphical study was also performed for locations; Gandak & Gomti (Fig. S4). It was noted that two locations (Tons & Brahmaputra) though have shape differentiation through lollipop graph (Fig. S5), does not indicate much shape differentiation over landmarks as the tail length at all landmarks were minimal while significant procrustes distance proves shapes differentiation. Similarly, procrustes distance was recorded as 0.08225527, while, Mahalanobis distance was found to be 5.2176 at p-value 0.0003 between Gandak and Gomti indicating significant shape distinction. Cross validation and shape differentiating lollipop graph indicates clear differentiation in this case between locations as the tail length at almost all landmarks is high (Figs. S6A, S6B).

Procrustes Shape Modularity over subset of landmarks

A set of nine landmarks identified were partitioned into two subsets to check for shape variation through procrustes shape modularity. Landmarks 1, 8, and 9 formed a subset that were related to body feature; hump, and landmarks 4, 5, and 6 formed another subset related to fins. The shape modularity hypothesis of landmark subsets (Figs. S6A, S6B) indicated that both the subsets identified had higher relative variance coefficient (0.85, 0.85) (Table S9) which indicated the significant role of landmarks in body shape variations (Figs. S7A, S7B). The utility of procrustes shape modularity in shape-based clustering of specimens over subset of landmarks (1,8,9 vs. 2,3,4,5,6,7) and subset of landmarks (4,5,6 vs. 1,2,3,7,8,9.) with the clustering criteria as relative variance (RV) coefficients greater or less than 0.85 was carried out (Figs. S8A, S8B).