Ten populations of Monacha pantanellii (Table 1, Fig. 1) were considered in our analysis of their molecular and morphological (shell and genitalia structure) variability, and compared with the M. cantiana s. l. lineages (Pieńkowska et al. 2018a, 2019b). The sequences deposited in GenBank were also considered for the molecular analysis (Table 2). Two other Monacha species were used for morphological and molecular comparison: M. cartusiana (Müller, 1774) and M. parumcincta (Rossmässler, 1834). Another 23 populations of M. pantanellii were studied on a qualitative morphological basis (they were not included in the statistical analysis) (Table 3).

Figure 1. 

Localities of Monacha pantanellii and M. cartusiana populations listed in Tables 1, 3 (M. pantanellii black circles Table 1, grey circles Table 3). Details of localities of other Monacha species and their molecular lineages were provided in previous papers (Pieńkowska et al. 2015, 2018a, 2018b, 2019b).

New material examined is listed as follows, when possible: geographic coordinates (Lat & Long and UTM references) of locality, locality (country, region, site, municipality and province), collector(s), date, number of specimens (sh/s shell/shells; spcm/spcms specimen/specimens), and collection where material is kept in parenthesis (Tables 1, 3). The material is kept in the F. Giusti collection (FGC; Dipartimento di Scienze Fisiche, della Terra e dell’Ambiente, Università di Siena, Italy). The material used for comparison has already been described (see Pieńkowska et al. 2018a: table 1, 2019b: table 1).

Fifty-eight specimens representing ten population of Monacha pantanellii were used for molecular analysis (Table 1). DNA extraction, amplification and sequencing methods are described in detail in our previous paper (Pieńkowska et al. 2018a).

Two mitochondrial and two nuclear gene fragments were analysed, namely cytochrome c oxidase subunit 1 (COI), 16S ribosomal DNA (16S rDNA), an internal transcribed spacer of rDNA (ITS2) and histone 3 (H3), respectively. All new sequences were deposited in GenBank (Table 1). The COI, 16S rDNA, ITS2 and H3 sequences obtained from GenBank for comparison are listed in Table 2.

The sequences were edited by eye using the programme BioEdit, version 7.2.6 (Hall 1999, BioEdit 2017). Alignments were performed using CLUSTALW (Thompson et al. 1994) implemented in MEGA7 (Kumar et al. 2016). The COI and H3 sequences were aligned according to the translated amino acid sequences. The ends of all sequences were trimmed. The lengths of the sequences after trimming were 592 bp for COI, 286 positions for 16S rDNA, 501 positions for ITS2 and 279 bp for H3. The sequences were collapsed to haplotypes (COI and 16S rDNA) and to common sequences (ITS2 and H3) using the programme ALTER (Alignment Transformation EnviRonment) (Glez-Peña et al. 2010). Gaps and ambiguous positions were removed from alignments prior to phylogenetic analysis. Mitochondrial (COI and 16S rDNA) and nuclear (ITS2 and H3) sequences were concatenated (Table 4) before phylogenetic analysis. Finally, the sequences of COI, 16S rDNA, ITS2 and H3 were concatenated (Table 4) for Maximum Likelihood (ML) and Bayesian Inference (BI).

Estimates of evolutionary divergence between the sequences of COI obtained in this study and other sequences from GenBank were conducted with MEGA7 using the Kimura two-parameter model (K2P) (Kimura 1980). The analysis involved 83 nucleotide sequences. All positions containing gaps and missing data were eliminated. There was a total of 615 positions in the final dataset.

Maximum Likelihood (ML) analyses were then performed with MEGA7. Monacha cartusiana and Monacha parumcincta were added as outgroup species in each analysis. For ML analysis of concatenated sequences, the following best nucleotide substitution models were specified according to the Bayesian Information Criterion (BIC): HKY+G+I (Hasegawa et al. 1985, Kumar et al. 2016) for COI and 16S rDNA concatenated sequences of 878 positions (592 COI + 286 16S rDNA), T92+G+I (Tamura 1992, Kumar et al. 2016) for ITS2+H3 concatenated sequences of 780 positions (501 ITS2 + 279 H3), and T92+G+I for COI+16S rDNA+ITS2+H3 concatenated sequences with a total length of 1658 positions (592 COI + 286 16S rDNA + 501 ITS2 + 279 H3). Bayesian analysis was conducted with MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003) using the evolution model already used for ML calculation. Four Monte Carlo Markov chains were run for one million generations, sampling every 100 generations (the first 250,000 trees were discarded as ‘burn-in’). This gave us a 50% majority rule consensus tree. In parallel, Maximum Likelihood (ML) analysis was performed with MEGA7 (Kumar et al. 2016) and calculated bootstrap values were mapped on the 50% majority rule consensus Bayesian tree.

One hundred and thirty-four specimens representing M. pantanellii, M. cantiana s. l., M. parumcincta and M. cartusiana were considered to investigate shell variability between these four species (including six molecular lineages of M. cantiana s. l.) (see Table 1 and Pieńkowska et al. 2018a, 2019b); the 43 specimens of nine populations of M. pantanellii (Fio1, Val, Sab, Alt, Lor, Tur2, Tur1, Car and Ani, see Table 1) were also considered to investigate shell variability between specimens of these populations. Shell variability was analysed randomly choosing five adult specimens from each population, when possible. Twelve shell variables were measured to the nearest 0.1 mm using Adobe Photoshop 7.0.1 on digital images of apertural and umbilical standard views taken with a Canon EF 100 mm 1:2.8 L IS USM macro lens mounted on a Canon F6 camera: AH aperture height, AW aperture width, LWfW last whorl final width, LWmW last whorl medial width, LWaH height of adapical sector of last whorl, LWmH height of medial sector of last whorl, PWH penultimate whorl height, PWfW penultimate whorl final width, PWmW penultimate whorl medial width, SD shell diameter, SH shell height, and UD umbilicus diameter (Pieńkowska et al. 2018a: fig. 1).

One hundred and thirty-five specimens of M. pantanellii, M. cantiana s. l. (with its six molecular lineages), M. parumcincta and M. cartusiana were analysed to examine anatomical variability between species; the 50 specimens of ten populations of M. pantanellii were also considered to investigate genital variability between populations of this species (see Table 1 and Pieńkowska et al. 2018a, 2019b). Snail bodies were dissected under the light microscope (Wild M5A or Zeiss SteREO Lumar V12). Anatomical details were drawn using a Wild camera lucida. Acronyms: BC bursa copulatrix, BW body wall, DBC duct of bursa copulatrix, DG digitiform glands (also known as mucous glands), E epiphallus (from base of flagellum to beginning of penial sheath), F flagellum, FO free oviduct, GA genital atrium, GAR genital atrium retractor, OSD ovispermiduct, P penis, PP penial papilla (also known as glans), V vagina, VA vaginal appendix (also known as appendicula), VAS vaginal appendix basal sac, VS vaginal sac (only present in M. cartusiana; see Pieńkowska et al. 2015: figs 11, 12), VD vas deferens. Seven anatomical variables (DBC, E, F, P, V, VS, VA) were measured under a light microscope (0.01 mm) using callipers (see: Pieńkowska et al. 2018a: fig. 2).

Detailed methods of multivariate ordination by Principal Component Analysis (PCA) and Redundancy Analysis (RDA), performed on the original shell and genitalia matrices as well as on the shape-related Z-matrices, are described in a previous paper (Pieńkowska et al. 2018a).

Differences between species for each shell and genital character were assessed through box-plots and descriptive statistics. Overall significance of differences was obtained using the Kruskal-Wallis test; when the test proved significant, multiple comparisons between pairs of species were performed using Dunn’s test. In order to control the false discovery rate (FDR), the Benjamini-Hochberg correction was used to adjust P-values for multiple comparisons. The dunn.test function with the altp = TRUE option and α = 0.01 in the dunn.test R package were used for analysis (Dinno 2017).

DNA sequencing resulted in 53 and 57 sequences of mitochondrial COI and 16S rDNA as well as 43 and 57 sequences of nuclear ITS2 and H3 gene fragments, respectively. They were all deposited in GenBank as MT380011–MT380063 (COI), MT376031–MT376087 (16S rDNA), MT376088–MT376130 (ITS2) and MT385776–MT385832 (H3) (Table 1). Thirty-six COI (COI 1–COI 36) and 25 16S rDNA (16S 1–16S 25) haplotypes, as well as 15 ITS2 (ITS2 1–ITS2 15) and 15 H3 (H3 1–H3 15) common sequences were recognised among them (Table 1). They were used for phylogenetic analysis with appropriate sequences representing M. parumcincta (PAR) and M. cartusiana (CAR), as well as six molecular lineages of M. cantiana s. l. (CAN-1–CAN-6) obtained from GenBank (Table 2). ML trees for concatenated sequences of mitochondrial COI and 16S rDNA (Fig. 2, Table 4) and of nuclear ITS2 and H3 (Fig. 3, Table 4) gene fragments, as well as the Bayesian Inference (BI) phylogenetic tree of concatenated sequences of COI+16S rDNA+ITS2+H3 gene fragments (Fig. 4, Table 4) clustered the concatenated sequences in one clade (PAN) separated from all other clades hitherto recognised for M. cantiana s. l. (CAN-1, CAN-2, CAN-3, CAN-4, CAN-5, CAN-6), M. parumcincta (PAR) and M. cartusiana (CAR) populations (Pieńkowska et al. 2018a, 2018b, 2019b).

Figure 2. 

Maximum Likelihood (ML) tree of concatenated COI and 16S rDNA haplotypes of Monacha pantanellii (see Table 4). New COI and 16S rDNA sequences of M. pantanellii (Table 1) were compared with COI and 16S rDNA sequences of M. cantiana s. l. and M. parumcincta obtained from GenBank (Tables 2, 4). Numbers next to branches indicate bootstrap support above 50% calculated on 1000 replicates (Felsenstein 1985). The tree was rooted with M. cartusiana concatenated sequences obtained from GenBank (Table 2).

K2P genetic distances between COI haplotypes are summarised in Table 5. Differences in COI haplotypes of M. pantanellii are rather small (up to 4.5%). Three varied somewhat more (COI 8 from populations from Vallonina [Val] and Carsoli [Car], COI 30 from Valle del Turano [Tur1] and COI 32 from Carsoli [Car]), bringing the mean for all populations to 0.2–6.7%. It was not possible to differentiate one population from the others. It is noteworthy that haplotypes of M. pantanellii are very different (15.5–22.0%) from the others representing M. cantiana s. l. (i.e., M. cantiana CAN-1–CAN-3, M. cemenelea (Risso, 1826) CAN-4, and M. sp. CAN-5–CAN-6) as well as from M. parumcincta (18.1–21.4%) and M. cartusiana (16.6–18.3%) (Pieńkowska et al. 2018a, 2019b).

Figure 3. 

Maximum Likelihood (ML) tree of concatenated ITS2 and H3 common sequences of Monacha pantanellii (see Table 4). New ITS2 and H3 sequences of M. pantanellii (Table 1) were compared with ITS2 and H3 sequences of M. cantiana s. l. and M. parumcincta obtained from GenBank (Tables 2, 4). Numbers next to branches indicate bootstrap support above 50% calculated on 1000 replicates (Felsenstein 1985). The tree was rooted with M. cartusiana concatenated sequences obtained from GenBank (Table 2).

Figure 4. 

Bayesian 50% majority-rule consensus tree of the concatenated data set of COI and 16S rDNA haplotypes, and ITS2 and H3 common sequences (see Table 4). Sequences of M. pantanellii were compared with appropriate sequences of M. cantiana s. l. and M. parumcincta obtained from GenBank (Tables 2, 4). Posterior probabilities (left) and bootstrap support above 50% from ML analysis (right) are indicated next to the branches. Bootstrap analysis was run with 1000 replicates (Felsenstein 1985). The tree was rooted with M. cartusiana concatenated sequences obtained from GenBank (Table 2).

Monacha pantanellii has a globose to sub-globose shell, variable in size, colour, and presence of paler subsutural and peripheral bands, with roundish to oval slightly descending aperture, a brownish peristome and a very small to small umbilicus (Figs 5–31).

Figures 5–14. 

Shell variability in Monacha pantanellii from Monte Fionchi, summit [Fio1] (FGC 8140) (5, 6), Monte Fionchi, Torrecola [Fio2] (FGC 38944) (7), road to Montenero Sabino [Sab] (FGC 41552) (8–10) and Turania [Tur1] (FGC 42971) (11–14).

Figures 15–22. 

Shell variability in Monacha pantanellii from Lago del Turano [Tur2] (FGC 41654) (15–17), Via Salaria, Ornaro Alto [Alt] (FGC 41553) (18–20) and Vallonina [Val] (FGC 25345) (21, 22).

Figures 23–31. 

Shell variability in Monacha pantanellii from Via Salaria, Poggio San Lorenzo [Lor] (FGC 41551) (23–25), Valle dell’Aniene, Roccagiovine [Ani] (FGC 42974) (26–29) and Carsoli [Car] (FGC 41651(30, 31).

RDA with species or molecular lineage constraint on the shape and size matrix (Fig. 32) showed that RDA 1 (33%, P < 0.001) separated all the species or molecular lineages from PAR. The preliminary classic PCA revealed size as the first major source of morphological variation, since PC1 (70%) was a positive combination of all variables. On the contrary, RDA 2 (7.3%, P < 0.001) slightly separated CAN-1, CAN-2 and CAN-3 from CAN-4, CAN-5, CAN-6 and PAN with PAR in intermediate position. In this regard, PC2 (13%) mostly accounted for contrast between LWmH vs. LWaH and PWH.

RDA on the shape (Z) matrix (Fig. 33) showed a hazier separation of species or molecular lineages, confirming that size is a major source of morphological variation, although both RDA axes proved to be significant. In particular, RDA 1 separated CAR, CAN-5, CAN-6 from PAR, CAN-1 and CAN-3, with the other groups in a more or less intermediate position. Conversely, RDA 2 separated PAR and CAR from all the other species or molecular lineages. Shape-related PCA indicated that SH, LWaH and PWH vs. LWfW were the principal shape determinants on PC1 and PWmW, AH and AD vs. UD on PC2.

Figures 32, 33. 

Principal component analysis (PCA) and Redundancy analysis (RDA) with species or molecular lineage constraint applied to the original shell matrix (32) and shape-related Z-matrix (33).

Box plots (Fig. 34) proved that the shell characters only have discriminating value in distinguishing Monacha pantanellii from other species or molecular lineages in a few cases. In fact, according to Dunn’s test with Benjamini-Hochberg adjustment (α = 0.01), no character significantly distinguished PAN from CAN-1, CAN-2 and CAN-4, only one distinguished it from CAN-5 (UD), only two from CAR (LWah, PWH), four from CAN-6 (SD, LWmW, LWfW, UD), six from CAN-3 (SH, AH, SD, AD, LWmW, PWmW) and eight from PAR (SH, AH, SD, AD, LWmW, PWfW, LWfW, UD) (Table 6).

Figure 34. 

Box plots for shell characters of the nine Monacha species or molecular lineages investigated. The lower and upper limits of the rectangular boxes indicate the 25^th to 75^th percentile range, and the horizontal line within the boxes is the median (50^th percentile).

RDA with population constraint on the shape and size matrix (Fig. 35) showed that RDA 1 (53.6%, P < 0.001) separated them into two groups, the first, including populations from Via Salaria, Ornaro Alto [Alt], Valle dell’Aniene, Roccagiovine [Ani], Monte Fionchi, summit [Fio1], Via Salaria, Poggio San Lorenzo [Lor], Montero Sabino [Sab] and Lago del Turano [Tur2] was separate from the second consisting of populations from Carsoli [Car], Turania [Tur1] and Vallonina [Val]. On the contrary, RDA 2 (4.0%, P > 0.05) showed no significant separation of populations. The preliminary classic PCA revealed size as the first major source of morphological variation, since PC1 (67.0%) was a negative combination of all variables.

RDA on the shape (Z) matrix (Fig. 36) showed no significant separation between populations, again confirming that size is a major source of morphological variation. Shape-related PCA indicated that LWaH and PWH vs. LWfW were the two principal shape determinants on PC1 and AH vs. LWmH on PC2.

Figures 35, 36. 

Principal component analysis (PCA) and Redundancy analysis (RDA) with population constraint applied to the original shell matrix (35) and shape-related Z-matrix (36) of specimens of Monacha pantanellii.

Monacha pantanellii has distal genitalia very similar to those of the Monacha cantiana group. The most remarkable features are the usually short vaginal appendix with mid or proximal vaginal insertion, the long flagellum and the penial papilla with thick external wall bordering a central duct without strips joining it to the external wall and with a lumen filled by many variably sized pleats (Figs 37–63).

Figures 37–40. 

Genitalia (proximal parts excluded) (37), internal structure of distal genitalia (38), transverse sections of medial epiphallus (39) and apical penial papilla (40) of Monacha pantanellii from Monte Fionchi summit [Fio1] (FGC 8140).

Figures 41–44. 

Genitalia (proximal parts excluded) (41), internal structure of distal genitalia (42), transverse sections of medial epiphallus (43) and apical penial papilla (44) of Monacha pantanellii from Monte Fionchi, Torrecola [Fio2] (FGC 38944).

Figures 45–48. 

Genitalia (proximal parts excluded) (45), internal structure of distal genitalia (46), transverse sections of medial epiphallus (47) and apical penial papilla (48) of Monacha pantanellii from Vallonina [Val] (FGC 25345).

Figures 49–52. 

Genitalia (proximal parts excluded) (49), internal structure of distal genitalia (50), transverse sections of medial epiphallus (51) and apical penial papilla (52) of Monacha pantanellii from Valle dell’Aniene, Roccagiovine [Ani] (FGC 42974).

Figures 53–56. 

Genitalia (proximal parts excluded) (53), internal structure of distal genitalia (54), transverse sections of medial epiphallus (55) and apical penial papilla (56) of Monacha pantanellii from Carsoli [Car] (FGC 41651).

Figures 57–59. 

Genitalia (proximal parts excluded) (57), internal structure of distal genitalia (58) and transverse section of apical penial papilla (59) of Monacha pantanellii from Lago del Turano [Tur2] (FGC 41654).

Figures 60–63. 

Internal structure of distal genitalia of Monacha pantanellii from Valle dell’Aniene, Roccagiovine [Ani] (FGC 42974) (60), Via Salaria, Ornaro Alto [Alt] (FGC 41553) (61), Turania [Tur1] (FGC 42971) (62) and road to Montenero Sabino [Sab] (FGC 41552) (63).

RDA with species or molecular lineage constraint on the shape and size matrix (Fig. 64) showed that RDA 1 (27%, P < 0.001) separated M. cantiana s. l. (CAN-1, CAN-2, CAN-3, CAN-4, CAN-5 and CAN-6) from PAN, PAR and CAR. The preliminary classic PCA revealed size as the first major source of morphological variation, since PC1 (43%) accounted for VS vs. all the other variables. On the contrary, RDA 2 (22%, P < 0.001) separated CAN-5, CAN-6 and PAN from CAR and PAR. The group CAN-1, CAN-2, CAN-3 and CAN-4 was in intermediate position. In that regard, PC2 (20%) accounted for a contrast between E and VA vs. F, P, V and VS.

RDA with species or molecular lineage constraint on the shape (Z) matrix (Fig. 65) showed that RDA 1 (43%, P < 0.001) separated PAN from the group CAN-1, CAN-2, CAN-3, CAN-4 and CAN-6, with CAN-5, PAR and CAR in intermediate position, and that RDA 2 (20%, P < 0.001) separated CAR from all the others. Shape-related PCA indicated that VA and E vs. V and P were the principal shape determinants on PC1 and VS and V vs. DBC and F on PC2. In the latter case, removing the size effect altered the overall relationship patterns.

Figures 64, 65. 

Principal component analysis (PCA) and Redundancy analysis (RDA) with species or molecular lineage constraint applied to the original genital matrix (64) and shape-related Z-matrix (65).

Box plots (Fig. 66) for anatomical characters showed that VA, F and P have the best discriminating value in distinguishing PAN: they distinguished 6 (VA) and 5 (F and P) species or molecular lineage pairs, respectively, according to Dunn’s test with Benjamini-Hochberg adjustment (α = 0.01), followed by E and V with four and three species or molecular lineage pairs, respectively (Table 6). The most recognisable pairs were PAN vs. CAR and PAN vs. CAN-1 (four significant characters), PAN vs. CAN-2, PAN vs. CAN-3, PAN vs. CAN-5 and PAN vs. PAR (3 significant characters). Only two characters significantly distinguished PAN vs. CAN-4 and only one PAN vs. CAN-6 (Table 6). Anatomical characters have high discriminating value as testified by very low p values after Dunn’s test: in most cases (19 of 22) P < 0.001 (Table 6).

Figure 66. 

Box plots for genital characters of the ten Monacha species or molecular lineages investigated. The lower and upper limits of the rectangular boxes indicate the 25^th to 75^th percentile range, and the horizontal line within the boxes is the median (50^th percentile).

RDA with population constraint on the shape and size matrix (Fig. 67) showed that RDA 1 (64%, P < 0.001) separated populations Carsoli [Car], Monte Fionchi, Torrecola [Fio2], Turania [Tur1] and Vallonina [Val] from populations Via Salaria, Ornaro Alto [Alt], Valle dell’Aniene [Ani], Via Salaria, Poggio San Lorenzo [Lor] and Lago del Turano [Tur2], with Monte Fionchi, summit [Fio1] and Montero Sabino [Sab] in intermediate position. The preliminary classic PCA revealed size as the first major source of morphological variation, since PC1 (65%) was a positive combination of all variables. On the contrary, RDA 2 (13%, P < 0.001) separated population Val from populations Sab and Fio1, with all the other populations (Car, Fio2, Tur1, Alt, Ani, Lor and Tur2) in intermediate position. In that regard, PC2 (16.5%) accounted for a contrast between V and VA vs. F and DBC.

RDA on the shape (Z) matrix (Fig. 68) showed a less clear separation between populations. RDA 1 (43%, P < 0.001) separated the population Sab from the group of populations Tur2, Val, Alt and Lor, with Ani, Fio1, Fio2 and Car in intermediate position. Shape-related PCA indicated that V vs. F were the two principal shape determinants on PC1 (39.5%). RDA 2 (14%, P < 0.001) separated Tur2 from Alt, Lor and Fio1, with all the other populations in a more or less intermediate position. In that regard, PC2 (24.5%) accounted for a contrast between PD and VA.

Figures 67, 68. 

Principal component analysis (PCA) and Redundancy analysis (RDA) with population constraint applied to the original genital matrix (67) and shape-related Z-matrix (68) of specimens of Monacha pantanellii.