Geometric morphometric analysis of gonopods in Bicoxidens flavicollis populations (Diplopoda, Spirostreptida, Spirostreptidae)

Male gonopods are useful in taxonomic diagnoses and descriptions of millipedes, although they may vary intraspecifically in shape and size. To assess this intraspecific variation, we used geometric morphometric analysis to compare gonopod morphology among eight isolated populations of the colour-polymorphic southern African millipede Bicoxidens flavicollis . Our results showed that gonopod shapes vary significantly across the examined populations, and elucidated subtle variations. CVA cross-validation test indicates an average classification rate of 75% for the five populations for which we had more than one specimen. Although we had a small number of replicates for three populations, our results still illustrate the importance of applying quantitative approaches to what would otherwise be qualitative and subjective gonopod shape categories in millipedes. As such, the taxonomic assignment of the populations of B. flavicollis may require further investigation, and further revisions would be required with an integrative approach, including molecular data, in order to re-evaluate the taxonomic diversity and distribution data of this species. Finally, we highlight the conservation potential of divergent populations as evolutionary insurance against a dynamic and unpredictable climate, whether or not they undergo full speciation. Our results support the position that population divergence and variation in male genitalia of B. flavicollis could be coupled. Bond et al. (2003) concluded based on gonopod morphology that the classification of millipedes may overlook taxonomic diversity due to OTU lumping. Although this might be true in some genera, it does not seem to be the case in B. flavicollis . In order to better assess the taxonomic importance of male gonopod shape variation in tropical millipedes, geometric morphometrics could be included routinely in these studies. Our results have conservation implications for B. flavicollis in particular, and millipedes in general. Populations of the same species in different environments represent opportunities for maintaining a large adaptive genetic base that may, or may not, proceed to full speciation. Conserving such populations and the processes that maintain them is akin to procuring an evolutionary insurance against the vagaries of unpredictable changes in the environment.


Introduction
Morphological differences are important in taxonomy to delimit taxa (Jacob et al. 2004;Schlick-Steiner et al. 2007;Foley et al. 2017). In millipedes, gonopod morphology is central in species description because the male gonopods are divergent and speciesspecific. Recent studies have demonstrated that speciation may occur without changes in gonopod morphology in several invertebrates (see Bond and Sierwald 2002;Novo et al. 2010). Bond et al. (2003) reported that speciation occurred unaccompanied by gonopod divergence in a spirobolidan millipede, Anadenobolus excisus Karsch, 1881 species complex, and genetic change is uncorrelated to gonopod change. Although intraspecific variation in gonopods is common and differences in gonopod structure among taxa may be subtle (Pimvichai et al. 2011), male gonopods remain the key sources of traits that are used in millipede systematics (cf. Enghoff 2017;Frederiksen et al. 2012). Features of gonopods with consistent application in millipede taxonomy include the sternum, coxites and telopodite (Pimvichai et al 2011;Enghoff 2017), although Krabbe and Enghoff (1985) reported a case of considerable variability in gonopod coxite within a localised population of Orthoporus antillanus, suggesting that the variation in the character may not be exclusively geographic in origin.
The genus Bicoxidens Attems, 1928 comprises nine species which are endemic to southern Africa, occurring in diverse habitats such as Miombo woodland, riverine, and montane vegetation (Mwabvu et al. 2007). The species Bicoxidens flavicollis Attems, 1928 sensu lato is the most widely distributed species, with records from Zimbabwe and one from western Mozambique ( Fig. 1) (Mwabvu et al. 2007). Although the male body size and gonopod shape in B. flavicollis populations appear identical, these observations have been based only on qualitative comparisons. According to Mwabvu et al. (2007), the pattern of body colour is related specifically to the locality of a determined population, with black, brown, or orange-yellow specimens having been recorded only in Zimbabwe, for instance. Nonetheless, quantitative comparisons of the gonopods are necessary in order to clarify the circumscription of this putative species complex. Furthermore, investigating intraspecific variations in gonopod morphology in populations of B. flavicollis could identify isolated populations and may lead to the recognition of undescribed diversity and a revision of distributional data.
Given the colour polymorphism and high genetic divergence (18% for CO1 and 6% for 16S rRNA) recorded in two populations of B. flavicollis (Mwabvu et al. 2013), we assessed the variation of gonopod morphology in differentiating colour polymorphic populations. We compared the shapes of gonopods from eight populations of B. flavicollis, predicting that gonopod shape would be highly conserved, and despite the variation in body colour, populations of B. flavicollis would not be differentiated using traditional comparative morphology of gonopods.
From this perspective, we employed in this study for the first time a geometric morphometric analysis to assess the taxonomic value of gonopod morphology in separating populations, and to test the hypothesis that gonopod morphology may underestimate the taxonomic diversity of a spirostreptidan species.

Materials and methods
Copulating and non-copulating males and females of B. flavicollis were collected by hand from eight geographically separate populations in Zimbabwe as summarised in Table 1. All the specimens were preserved in 96% ethanol. The voucher specimens from each population were deposited in the Natural History Museum (NHMZ), Bulawayo, Zimbabwe, and the KwaZulu-Natal Museum (KZNM), Pietermaritzburg, South Africa. The total sample size on which our analysis was based was 40 (Table 1). The specimens were identified by T. Mwabvu after examining the male gonopods under a Carl Zeiss (Stemi DV4) dissecting microscope and using the species identification keys in Mwabvu et al. (2007). Male gonopods were photographed using Auto Montage software (Leica Microscope MZ12s with 3 CCD Toshiba Camera).
All the samples were photographed in the same oral view and magnifications. In total, eighteen landmarks for 40 specimens were obtained from the examined gonopods (Table 2; Fig. 2) using a suite of TPS programmes (TPS Util and TPSDig, Rohlf 2015). Considering the bilateral symmetry and to avoid any difference regarding gonopodal position only the right side was digitized. Procrustes standardisation was performed on the landmark coordinates to remove the effect of size, rotation, and position. To minimise the possibility of committing a type II error, we assessed measurement error on digitised landmarks using the intraclass correlation coefficient approach, also termed repeatability (Fisher 1958). To measure the repeatability of placing landmarks, an individual specimen from three randomly selected populations  Apex of median lobe of apical proplica I 2 Apex of axe-shaped process (or apical axe-shaped process) I 3 Tip of lobe of axe-shaped process I 4 Base of lobe of axe-shaped process I 5 Lateral mid-proplica II (semi-landmark; midpoint between landmarks 4 and 6) 6 Apex of sternite (or apical sternite) I 7 Base of sternite (or basal sternite) I 8 Base of proplica (or basal proplica) II (semi-landmark; the most distal part of the basal prolica) 9 Basal lateral edge of paracoxite II (semi-landmark; furthest point on the lateral edge of paracoxite 10 Distal lateral paracoxite I 11 Paracoxite apex I 12 Telopodite (at midlength between knee and femoral lobe) II (semi-landmark; midpoint between knee and femoral lobe) 13 Base of rounded lobe on proplica I 14 Base/pit of telopodite knee I 15 Apex of telopodite knee I 16 Proximal edge of lateral process I 17 Apex of lateral process I 18 Apical groove between lateral process and median lobe I (14 of the landmarks were type I; 4 were type II sensu Bookstein 1997a, b).
(Chihota, Mazowe and Sahumani) was used to provide six replicates of gonopod images in a total of 18 images for each one. We then landmarked each replicate and performed a one-way ANOVA on the Procrustes standardised landmark coordinates using population as the categorical variable. Repeatability was computed according to the formula: where R is the repeatability value, S A 2 is the among-individuals component of variance, and S W 2 is the within-individuals variance component (Sokal and Rohlf 1995;Arnqvist and Martensson 1998;Fruciano 2016). For shape analysis, the variance-covariance matrix derived from the Procrustes coordinates were subjected to exploratory principal component analysis to reduce the dimensionality of the data to the most significant shape variables. As the first 15 principal components (PCs) accounted for about 95% of observed variance, the corresponding PC scores were selected as variables for a canonical variates analysis (CVA) to assess shape discrimination across populations. This was followed with cross-validation tests to assess the rate of correct classification of specimens into their population groups. We tested the null hypothesis that gonopod shapes are not different among the eight populations of B. flavicollis using the PC scores and population as dependent and predictor variables, respectively. All analyses were implemented in PAST (Hammer et al. 2001) and MorphoJ (Klingenberg 2011).

Results
The assessment of measurement errors arising from landmark digitisation demonstrates that variation between repeated digitisations of the same specimen is significantly lower than variation among different specimens, with an intraclass correlation coefficient of 0.99. The thin plate spline deformation grids (Fig. 3) depict some patterns of gonopod shape variations among the 18 landmarks, as well as among the various populations of B. flavicollis. Landmarks 6, 7 and 11 displayed the most visually perceptible variability. The overall shape changes along the first two principal component axes also identified the same sets of landmarks, in addition to semi-landmarks 8 and 9, as showing pronounced variability across all populations. Although a scatter plot of specimens along PC1 and PC2 (Fig. 3) does not show sufficient discrimination of the populations, the canonical variate analysis indicates a significant effect of population on gonopod shape (Wilk's λ = 0.0029, F = 4.17, df = 60 and 73, P = 7.36 × 10 -9 ). PC1 and PC2 represent 34.4% of the total shape variation and correspond largely to the sternite/metaplica and the paracoxites, respectively. The gonopod shape space occupied by each specimen, as shown on the CV1 and CV2 axes, suggests clustering along population trajectories (Fig. 4).

Discussion
For any morphological structure in an organism, there are almost limitless possibilities of theoretical shapes that could evolve. However, only a few are evolutionarily viable due to functional, and thus adaptive constraints placed on such structures by natural selection (Mcghee 2007;Adebowale et al. 2012). These adaptive constraints may explain in part why comparatively only a limited set of shape configurations are encountered across the biological world, even among lineages with diverse genetic heritage.
Given this background and the significance of male gonopods in the reproductive success of arthropods, some levels of morphological stasis in gonopods would be expected in a species as widely distributed as B. flavicollis. Our results demonstrate that gonopod shape in B. flavicollis varies significantly among the eight populations, with landmarks 6, 7 and 11 showing pronounced shape changes along the first two PCs. However, it is conceivable that extensive sampling of Chihota and Muterere specimens would have blurred the distinction between them and their respective nearest neighbour population groups. The combination of landmarks 6 and 7, corresponding to the apex and base of the sternite, defines a spatial relation between two points that could be discriminatory at interspecific level of comparison. These three landmarks (6, 7, and 11) are regarded as type I (see Bookstein 1997a), and as such could be expected to hold some taxonomic or evolutionary significance for B. flavicollis. This is consistent with the finding of Enghoff (2017), who reported that the shapes of the sternite on gonopods of Spirostreptidae species are sufficiently variable to be of taxonomic value. While it is tempting to simplify the diagnostic attributes of the gonopod to the sternite, the proplica and the processes on the metaplica (however useful these are), it is equally noteworthy that the next eight PCs (PC3-PC10) capture 51.5% of shape variation. This indicates a more nuanced aspect to gonopod shape variation, and it is spread across most of the landmarks.
Based on the similarities of gonopods in a widely distributed African millipede genus Doratogonus Attems 1914, Hamer and Slotow (2000) reported that the diagnoses of the species could be too inclusive. Studies in other invertebrate groups, such as scorpions (Jacob et al. 2004), seem to support the assertion that morphology underestimates taxonomic diversity. However, these conclusions were based on qualitative assessments of gonopods, which is susceptible to subjective interpretations. A major strength of geometric morphometric methods is the quantification and analysis of shape information in a reproducible manner. This approach allows a more objective evaluation of important but cryptic shape variations that may be imperceptible to the human eye. Such cryptic variations, if heritable, could be the difference between operational taxonomic units (OTU). Furthermore, recent methodological advances in geometric morphometrics have seen the extraction of shape data and their incorporation into a phylogenetic matrix to infer evolutionary relationships (Smith and Hendricks 2013;Díaz -Cruz et al. 2021).
Although observations of conservatism in morphology have been reported in a number of animal taxa including spiders (Huber et al. 2005), copepods (Lajus et al. 2015) and rotifers (Fontaneto et al. 2007), our results support that gonopod shape differs among populations in B. flavicollis. A high degree of genetic divergence documented between two populations of B. flavicollis (see Mwabvu et al. 2013) suggests that gonopod morphology and genetic divergence could be coupled. Our results further demonstrate that gonopod morphology, when rightly quantified, could serve as an additional tool for differentiating populations or taxa. Besides the clear discontinuity in colour patterns, the eight populations in the present study occupy distinct geographic areas, which suggests a limited gene flow between populations. Thus, the variation in environmental conditions at these localities could have influenced the subtle but significant divergence in genital morphology. In addition to molecular data and body colour patterns, gonopods have reliable diagnostic traits that identify unique populations in the B. flavicollis. The evidence suggests that B. flavicollis could be an inclusive complex requiring further studies to clarify the taxonomic assignment of the populations. On the other hand, the species could have a genetic machinery that allows for some degree of morphological latitude in gonopod divergence in differing environments.
Future work on B. flavicollis should include the use of nuclear markers to compare the levels of genetic divergence among the populations. Based on the similar pattern of body colour between populations, the high levels of genetic divergence reported by Mwabvu et al. (2013), and the differences in gonopod shapes observed in our study, B. flavicollis could be a species complex. Importantly, Pimvichai et al. (2011) reported identical gonopods in some species of Thyropygus Pocock, 1894 (Harpagophoridae) that had high DNA sequence divergence, and Frederiksen et al. (2012) observed identical gonopods in millipede populations of Julidae (Julida) that differ in body colouration and DNA sequences. These observations indicate that genetic diversity may not necessarily be matched by morphological changes, and strengthens the argument for thorough molecular analyses of B. flavicollis populations.

Conclusion
Even though three out of eight populations in our study had one specimen, our study highlights the importance of using quantitative methods in taxonomy. Our results support the position that population divergence and variation in male genitalia of B. flavicollis could be coupled. Bond et al. (2003) concluded based on gonopod morphology that the classification of millipedes may overlook taxonomic diversity due to OTU lumping. Although this might be true in some genera, it does not seem to be the case in B. flavicollis. In order to better assess the taxonomic importance of male gonopod shape variation in tropical millipedes, geometric morphometrics could be included routinely in these studies. Our results have conservation implications for B. flavicollis in particular, and millipedes in general. Populations of the same species in different environments represent opportunities for maintaining a large adaptive genetic base that may, or may not, proceed to full speciation. Conserving such populations and the processes that maintain them is akin to procuring an evolutionary insurance against the vagaries of unpredictable changes in the environment.