Acerentomon italicumNosek, 1969 is the most abundant Protura species in Italy. Two groups of populations characterized by the different position of a pore on tergite VII were analyzed. Specimens from Austria and from the Italian provinces of Brescia and Bolzano (Lombardy and Alto Adige, respectively) possess a postero-sub-medial pore (psm); in specimens from Trento (Trentino), Veneto, Friuli-Venezia Giulia and Liguria the pore is in an antero-sub-lateral (asl) position. An integrative approach was applied to assess other characters distinguishing these populations. The lengths of body, head, foretarsus, claw of anterior leg, foretarsal sensilla a and b were measured; in addition some ratios were calculated: TR (length of foretarsus/length of its claw), BS (distance of sensillum t1 from base of foretarsus/distance of t1 from foretarsus tip), length of a/foretarsus length, length of b/foretarsus length and length of a/length of b. A geometric morphometric analysis was performed on the position of setae on tergite VII, selecting as landmarks the bases of setae A1, A2, A3, P1, P1a, P2, P2a and P3. Nucleotide sequences COI-5P and 28S-D2-D3 of both groups were compared. Significant differences in the foretarsus and anterior claw length, in the arrangement of setae on tergite VII and in the barcode sequences COI-5P between the “psm and asl populations” are shown to exist. Based on the results obtained, the two groups are different, but there is not enough support to differentiate them at the species level.
INTRODUCTION
Acerentomon italicum Nosek, 1969 is the most abundant Protura species in Italy (Galli et al., 2016) and it is widespread in Switzerland, Austria, Corsica and Slovenia. A recent redescription of the species (Galli et al., 2016) highlighted the difference between two groups of populations with regards to the porotaxy of tergite VII. Porotaxy was used as a diagnostic character to distinguish Protura species for about thirty years (e.g. Szeptycki, 1991), and it was already pointed out that the position of a single pore can differ between two closely related species. A good example is the difference in porotaxy between Gracilentulus gracilis Berlese, 1908b and G. europaeus Szeptycki, 1993 highlighted by Szeptycki (1993). Throughout most of its range (including the type area) specimens of A. italicum possess a single antero-sub-lateral (asl) pore close to the base of seta A3. Specimens from Switzerland, Austria, Lombardy, Friuli-Venezia Giulia and Slovenia, however, have a postero-sub-medial (psm) pore close to the base of seta P1 instead. The latter is the typical position of such a pore in the Acerentomidae (Shrubovych, 2014). The remaining porotaxy and chaetotaxy of these two groups is largely identical, with a certain level of variability that characterizes all species. Moreover, DNA barcoding showed a remarkable distance between a group of three Ligurian populations (asl) and one from Austria (psm) (Galli et al., 2016). These results prompted us to look for further characteristics that can distinguish the two “groups of populations” (from here on indicated as asl and psm populations). In this paper the results of an integrative study, based on genetic data, specific morphometric analysis and through a geometric morphometric analysis, are shown. Geometric morphometrics, in particular, is a groundbreaking approach in Protura that deserves more attention in the future.
MATERIAL & METHODS
Morphometry
Abbreviations of depositories: MHNG = Natural History Museum of Geneva, Switzerland; UniGe = Genoa University collection, Italy; ISEAPol = Institute of Systematics and Evolution of Animals, Polish Academy of Sciences, Poland.
Material examined: MHNG; 2 ♀; Austria, Gaisberg. – MHNG; 1 ♂, 1 ♀; Austria, St Pankraz. – MHNG; 7 ♂, 7 ♀; Austria, Plabutsh. – ISEAPol; 1 ♀; Austria, Leopoldsberg, Vienna. – UniGe; 5 ♂, 13 ♀; Italy, Lombardy, Brescia Province, Polaveno. – MHNG; 1 ♂; Italy, Trentino-Alto Adige, Bolzano Province, Brunico. – MHNG; 2 ♀; Italy, Trentino-Alto Adige, Trento Province, Levico Terme. – MHNG; 5 ♂, 4 ♀; Italy, Veneto, Belluno Province, Ponte nelle Alpi. – MHNG; 3 ♂, 2 ♀; Italy, Friuli-Venezia Giulia, Udine Province, Paluzza. – MHNG; 1 ♂; Friuli-Venezia Giulia, Udine Province, surroundings of Tolmezzo. – UniGe; 1 ♂, 5 ♀; Italy, Liguria, Savona Province, Bergeggi. – UniGe; 5 ♂, 5 ♀; Italy, Liguria, Genoa Province, Genoa Quinto. For a map of the localities listed above, see Fig. 1.
Measurements: The presence of a pore in the asl or psm position was detected (Fig. 2). When possible, the following measurements (and ratios of two measurements) of each specimen were taken: body, head, foretarsus and claw of anterior leg lengths (in the following simply referred to as tarsus and claw, as commonly used in the literature on Protura); TR (length of tarsus/length of claw) and BS (distance of sensillum t1 from base of tarsus/distance of t1 from tarsus tip) ratios (see Nosek, 1973); foretarsal sensilla a and b lengths; ratio of length of a/tarsus length, ratio of length of b/tarsus length and ratio of length of a/length of b. Measurements were taken with the aid of an interference contrast microscope (Leica DM LB2), a Leica DFC 295 camera and a Leica Application Suite Vers. 3.8. Data were archived separately based on pore position (asl or psm populations) and sex of specimens.
Fig. 1.
Geographical distribution of examined populations of Acerentomon italicum (grey = with psm pore; black = with asl pore).
Statistical analyses: For each character measured outliers were removed before the application of the following analyses. On the measurements taken One-way ANOVA and Mann-Whitney pairwise test (adopting a Bonferroni correction of p-values) were performed to see differences between sexes within the same “population” and between values recorded in psm and asl specimens within each sex. For statistically significant variables a linear correlation with the geographic coordinates of sampling localities was performed to verify any latitudinal or longitudinal trend. Analyses were made using the software “PAST” (Paleontological Statistics) version 4.02 (Hammer et al., 2001).
Geometric morphometrics
Landmarks placement: The Leica system described above was also used to take photos of tergite VII of 53 specimens belonging to both populations (asl: 15 males and 18 females - psm: 9 males and 11 females). Specimens were selected in which on one side of tergite VII the landmarks mentioned below were coplanar, resulting in all being in focus at the same time. Based on these photos the coordinates of defined sets of landmarks were acquired for each specimen. The sets of coordinates (called landmark configurations) contain the information related to shape and size, orientation and location of setae on tergite VII. In particular, eight type I bidimensional landmarks (Bookstein, 1997) were defined: the bases of setae A1, A2, A3, P1, P1a, P2, P2a and P3 (Fig. 2). The TPS software was used for image storing and landmarks placing (TPSutil, version 1.74, see Rohlf, 2017a and TpsDig, version 2.30, see Rohlf, 2017b).
Further analyses were performed using the “geomorph” package (Adams et al., 2020), written in the “R” scientific computing language (R Core Team, 2020). Landmarks previously digitized with TPS Dig were imported into the R environment using the “readland.tps()” function.
Generalized procrustes analysis: The first step of any geometric morphometric analysis is to perform a Procrustes superimposition of the raw coordinate data (for a general overview of geometric morphometrics see Adams et al., 2013). A Generalized Procrustes Analysis (GPA, see Rohlf & Slice, 1990) translates all specimens to the origin, scales them to unit-centroid size, and optimally rotates them (using a least-squares criterion) until the coordinates of corresponding landmarks align as closely as possible. The resulting aligned Procrustes coordinates represent the shape of each specimen, and are found in a curved space related to Kendall's shape space (Kendall, 1984). To perform this analysis the “gpagen()” function was used. GPA allows to extract the shape information from all landmark configurations by removing the parameters of non-shape variation (i.e. differences in size, location and orientation). Information carried by location is removed by shifting the landmarks configurations until they share the same origin in a common coordinate system (Theska et al., 2020) where each landmark is represented by shape variables called “Procrustes coordinates”. A data-cleaning process using the “plotOutliers()” function led to remove from the resulting analyses four specimens that proved to be outliers: one male and two females belonging to the psm group, and one female belonging to the asl group. Aligned coordinates of the remaining 49 specimens were represented using the “plotAllSpecimens()” function that creates a plot of the landmark coordinates for all specimens. This is useful for examining patterns of shape variation after GPA. Selecting “mean=TRUE”, the mean shape was calculated and added to the plot (Fig. 3). Together with procrustened coordinates, the consensus and the centroid size were retrieved. Centroid size (CS) is the most common and explicit measure of size in geometric morphometrics. It is computed as the square root of the sum of the squared distances of all landmarks from their centroid (Rohlf & Slice, 1990; Goodall, 1991; Dryden & Mardia, 1998).
Fig. 2.
Acerentomon italicum, tergite VII: landmarks used for geometric morphometric analysis of asl (left) and psm (right) populations. Scale bar = 50 µm.
Principal Component Analysis: The Principal Component Analysis (PCA) is probably the most common method used to visualize the general patterns of morphological variation in multidimensional data (in this case the covariance matrix of the procrustened coordinates of landmarks). By definition, the first principal component (PC1) retains the maximum possible variance (Jolliffe, 2002; Abdi & Williams, 2010). The cumulative variance explained by each PC was examined and the relative screeplot was used to select the number of PCs for further analyses.
The Wilcoxon-Mann-Whitney test was performed to highlight differences between groups. Extreme shape variation along the PCs was explored using the “plotRefToTarget()” function. To assess the reliability of the PCA results, a cross validated (using leave-one-out jacknife procedure) linear discriminant analysis (LDA) (MASS package, Ripley et al., 2013) was run on extracted PCs. Discrimination was performed separately between populations and sexes.
Thus scores were considered as the projection of the centered data in the linear space defined by the eigenvectors (Jolliffe, 2002). Therefore they were analyzed to understand shape variations between groups into a dimensionally reduced space. Nevertheless the loadings were used to assess the weight of each variable in the description of the shape variation in every PC using the “fviz_contrib()” function from the “factoextra” package (Kassambara & Mundt, 2017).
Covariation method: To infer differences in shape according to factors or continuous variables, the “procD.lm()” function was used. It performs Procrustes ANOVA with permutation procedures to verify statistical hypotheses describing patterns of shape variation and covariation for a set of Procrustes-aligned coordinates. The mean shape variation among factors (populations and sexes) was tested by means of a Goodall's F test (Goodall, 1991). A second model (multivariate regression) was used to describe shape variation according to size. Since size in geometric morphometrics is represented by centroid size, a simple multivariate regression, where shape represents the response variable and the centroid size represents the independent variable, was fitted to the data. Lastly a multivariate analysis of variance (single factor MANCOVA) was used to test changes in shape between groups, while accounting for shape varying with size.
Genetic distances
Nucleotide sequences COI-5P and 28s-D2-D3 available on the Barcode of Life Data systems (BOLD - https://www.boldsystems.org/) for Acerentomon species belonging to doderoi-group (sensu Nosek, 1973) were analyzed: A. italicum, A. christiani Shrubovych & Resh in Shrubovych et al., 2016, A. maius Berlese, 1908b (Tables 1-2). Alignment of sequences was performed using “MUSCLE”, as implemented in Mega v. 7.0.20 (Kumar et al., 2016). The genetic distances based on Kimura 2-parameter (Kimura, 1980) were calculated within groups and between groups using the appropriate menu on Mega v. 7.0.20. In order to evaluate their taxonomic significance, the specimens belonging to the A. italicum psm population were considered as part of a separate Molecular Operational Taxonomic Unit (MOTU).
RESULTS
Morphometry
Values of the variables measured in asl and psm populations are shown in Table 3. The One-way ANOVA and Mann-Whitney pairwise test (adopting a Bonferroni correction of p-values) showed only two variables with statistically significant differences between psm and asl populations: tarsus length (F = 23.22; p << 0.001) and claw length (F = 12.18; p = 0.0025). Both characters have significantly larger values in psm than asl populations (Table 4, Fig. 4). For none of the characters it was possible to outline a latitudinal or longitudinal trend. Larger sizes of females than males within the same group were not statistically significant.
Table 2.
Accession numbers at Barcode of Life Data systems (BOLD) for specimens belonging to the Acerentomon doderoi-group whose nucleotide sequences 28S-D2-D3 were used in this work. Acerentomon italicum specimens without notes belong to Italian populations with an asl pore.
(continued)
Table 3.
Acerentomon italicum: morphometry of asl and psm populations. Mean values ± standard deviations; lengths in µm.
Geometric morphometrics
The first four components of PCA accounted for more than 73% of the total shape variation (Table 5). Also the screeplot inspection (Fig. 5) confirmed the same number of principal components as a reliable measure. The scoreplot based on the first two PCs shows that the two populations were quite different in shape (Fig. 6). Only for the first two PCs Wilcoxon-Mann-Whitney tests on differences between groups retrieved significant differences (PC1: U = 474, p = 5.83E-06; PC2: U = 360, p = 6.57E-02). When sexes were compared, no clear distinctions arose from between-group inference. Since only PC1 and PC2 accounted for shape variation, extreme shape deformations along these axes were plotted as a thin-plate spline deformation grid (Fig. 7). More than 40% of variation in PC1 was due to the horizontal coordinate of seta A2 (A2x), while PC2 shows a more even situation where P3 and A3 accounted for most of the variation (Fig. 8).
Table 4.
Mann-Whitney pairwise test (adopting Bonferroni correction of p-values). Output for tarsus length and claw length.
Fig. 4.
Boxplot of tarsus and claw length measurements in both sexes of psm and asl populations of Acerentomon italicum.
A linear discriminant analysis performed between populations confirmed the reliability of PCA, showing that the first four PCs have an overall classification accuracy of 83.7% (41 specimens correctly classified from a total of 49, Table 6), while discrimination between males and females did not provide meaningful results (51% overall accuracy).
The Goodall test retrieved significant differences in setae arrangement only between populations (F = 7.2242, p = 1e-04 with 1 and 48 df) but not between sexes (F = 0.7429, p = 6.17E-01 with 1 and 48 df). Despite centroid size showing a strong group structure both between populations (W = 35, p-value = 2.36E-08) and sexes (W = 412, p-value = 2.32E-02), multivariate regression highlighted a barely significant amount of covariation between shape and size (F = 1.91; p = 0.05 with 1 and 48 df).
Table 5.
Results of principal components analysis (PCA). For each PC the amount of explained variance is reported together with the cumulative proportions (both in %).
Genetic distances
Values of Kimura 2-parameter within-group and between-group genetic distances (Kimura, 1980) calculated on the basis of the nucleotide sequences COI-5P and 28s-D2-D3 available at BOLD for Acerentomon species belonging to the doderoi-group are shown in Tables 7-8. With regard to the barcode sequence COI-5P, the distance of 0.292 between the only available specimen of A. italicum belonging to the psm population and those with an asl pore is much higher than any within-group distance (Table 7) and it is comparable with between-group distances (Table 8). On the other hand, in the case of expansion segments 28s-D2-D3, the distance between the group of three A. italicum specimens with a psm pore (from Austria) and those with an asl pore is 0.077 (Table 8). This value is close to the within-group distance calculated for asl specimens (0.075 in Table 7) and lower than the other between-group distances (Table 8).
Fig. 6.
Scoreplots of the principal component analysis performed on 49 A. italicum specimens and 8 landmarks along the first two principal axes explaining 27.24% (PC1) and 24.17% (PC2) of the shape variation.
Fig. 7.
Extreme shape deformation for PC1 and PC2 maximum and minimum values. Arrows show deformations from mean shape (consensus).
Table 6.
Confusion matrix obtained performing LDA on first four PCs. The numbers of specimens correctly identified as belonging to asl (top left) and psm (down right) populations are given.
Table 7.
Genetic distances based on Kimura 2-parameter within Acerentomon species belonging to the doderoi-group (comparing COI-5P and 28S-D2-D3 sequences).
Table 8.
Genetic distances based on Kimura 2-parameter between Acerentomon species belonging to the doderoi-group (comparing COI-5P and 28S-D2-D3 sequences).
DISCUSSION
The traditional morphometric analyses for this paper were limited to few characters. We focused on those characters that, based on a recent redescription of A. italicum (Galli et al., 2016) and on our experience, seemed to be the most promising in order to highlight differences between psm and asl populations. In particular two foretarsal sensilla (a and b) are among the most important characters for Acerentomon species identification (see Shrubovych et al., 2016), especially in those species assumed to be close to A. italicum (see Galli et al., 2016). Therefore the absolute and relative length of these sensilla (to the tarsus and between them) was added to the structures that are traditionally measured. However, their between-group differences are not statistically significant. Just two size variables (tarsus and claw length) differ significantly between the two groups examined. Even though we did not outline any geographical trend for the examined variables, we cannot exclude that this may be due to the small sample examined. This could also be the reason for not having found differences in morphological characters which are of diagnostic importance in Acerentomidae (see Shrubovych, 2014).
Fig. 8.
Barplot showing the contribution of variables (expressed as percentage) along the first and second principal axes. The thick reference lines correspond to the expected value if the contributions were uniform. On the abscissa are the × and y coordinates of selected landmarks (setae A1 to A3 and P1 to P3, see Fig. 2).
Generally geometric morphometric analyses were applied to medium- and large-sized animals, and in hexapods they often dealt with wing venation. In the extensive literature dealing with this method, only few papers are about microscopic animals or about chaetotaxy (Mutanen et al., 2009; Ramírez-Sánchez et al., 2016; Theska et al., 2020). This study therefore is the first one in which a geometric morphometric approach is applied to Protura and one of the first ever to deal with setae arrangement and microscopic animals. Geometric morphometrics allows to analyze shapes without having almost any information about size. This approach avoids potential bias generated by intra- and even inter-population variation that can be caused by sexual dimorphism, which is quite common among arthropods. Both explorative (PCA) and inferential (ANOVA) techniques highlight differences between the arrangement of setae on tergite VII in the two groups examined. Results from PCA showed that the largest variation is found in the arrangement of setae A2, A3 and P3. Shape-related variation apparently does not depend on the variation of size-related elements as shown by regressions. When groups were compared, significant differences between populations were clearly highlighted by the Goodall test.
Both morphometry and geometric morphometrics do not show significant differences between sexes within each group. This confirms what was already highlighted by Galli et al. (2016) who examined all morphometric characters of A. italicum traditionally used in Protura taxonomy. However, Table 3 shows that females are on average slightly bigger than males. This size difference was evidenced also in other species of Acerentomidae like Acerentulus shrubovychae Galli & Capurro, 2013 and Acerentulus tortii Galli et al., 2017.
The existence of species complexes in Protura was already suggested by Fratello & Sabatini (1989) on the basis of differences between kariotypes of different Eosentomon transitorium Berlese, 1908a populations. A species-level diversification can also be assumed for allopatric populations of some species that were likely separated by ancient vicariance events, as it is the case in the amphi-Atlantic Delamarentulus tristani (Silvestri, 1938) (see Galli & Rellini, 2020). Theoretically DNA barcoding could be an excellent tool to distinguish cryptic species, and in some cases it was used for that purpose since the early days (e.g. Herbert et al., 2004). On the other hand, it was soon shown that especially for insufficiently sampled taxa (such as Protura) barcoding performs poorly when used for species identification (e.g. Meyer & Paulay, 2005). Our results must therefore be considered with extreme caution, especially in relation to the conflicting data from the analysis of COI-5P and 28S-D2-D3 sequences. For the latter some studies show that the more variable D2 region is more suitable not only for inferring relationships but also for taxonomic and diagnostic purposes at species level (e.g. Schmidt et al., 2006).
CONCLUSIONS
Our integrative analyses show that psm and asl populations of A. italicum differ for some characters other than the porotaxy of tergite VII. In particular psm specimens are characterized by significantly longer foretarsi and claws. The geometric morphometric analysis of the arrangement of setae on tergite VII allows to discriminate between these populations both in descriptive and inferential terms, leaving no doubt about the shape difference between asl and psm groups. Their DNA barcode sequences COI-5P differ at an interspecific level. These data show a clear difference between the studied groups. However, it is worth pointing out that only a small sample from the available localities has been studied, with few specimens selected for each. It would be better to include more localities and more specimens. Moreover, only a single DNA sequence was available for the psm group. At present there is therefore not enough evidence to differentiate the two groups at the species level with reasonable certainty.
Nevertheless this study is the first one in which a geometric morphometric approach was applied to Protura, and one of the first ever to deal with the arrangement of setae and microscopic animals. Beyond the result in itself, it is important to show that this methodology can be applied to the study of evolutionary patterns and trends of chaeto- and porotaxy in Protura in addition to the traditional morphometric approach and to the genetic one. This is particularly useful when genetic data indicate significant differences between populations which supposedly belong to the same species.
ACKNOWLEDGEMENTS
We wish to thank Dr Peter Schwendinger of the Natural History Museum of Geneva for kindly providing us the specimens we analyzed for this paper. Special thanks go to Dr Yun Bu for valuable suggestions that helped greatly improving the manuscript.
