Exploration of phylogeography of Monachacantiana s.l. continues: the populations of the Apuan Alps (NW Tuscany, Italy) (Eupulmonata, Stylommatophora, Hygromiidae)

Abstract Two new lineages CAN-5 and CAN-6 were recognised in four populations of Monachacantiana (Montagu, 1803) s.l. from the Italian Apuan Alps by joint molecular and morphological analysis. They are different from other M.cantiana lineages known from English, Italian, Austrian and French populations, i.e. CAN-1, CAN-2, CAN-3 and CAN-4, as well as from the other Italian Monacha species used for comparisons (M.parumcincta and M.cartusiana). Although a definite taxonomic and nomenclatural setting seems to be premature, we suggest that the name or names for these new lineages as one or two species should be found among 19th century names (Helixsobara Mabille, 1881, H.ardesa Mabille, 1881, H.apuanica Mabille, 1881, H.carfaniensis De Stefani, 1883 and H.spallanzanii De Stefani, 1884).


Introduction
Study of the phylogeography of the Monacha cantiana (Montagu, 1803) s.l. by joint molecular and morphological analysis revealed a number of cryptic lineages, some of which might deserve distinct taxonomic status.
Examination of a first group of English, Italian, Austrian and French populations showed that it consisted of at least four distinct lineages (CAN-1, CAN-2, CAN-3, CAN-4) . One of these lineages (CAN-1) included most of the UK (5 sites) and Italian (5 sites) populations examined. Three other lineages were represented by populations from two sites in northern Italy (CAN-2), three sites in northern Italy and Austria (CAN-3) and two sites in south-eastern France . A taxonomic and nomenclatural setting is only currently available for CAN-1 and CAN-4. The lineage CAN-1 corresponds to the true M. cantiana (Montagu, 1803) because it is the only one that includes topotypical English populations. The lineage CAN-4 is attributed to M. cemenelea (Risso, 1826), for which a neotype has been designated and deposited. A definitive frame for the other two has been postponed because it requires much more research.
We have now studied some populations from the Apuan Alps at the north-western extremity of Tuscany, a well-known hotspot of diversity and endemism (Lanza 1997;Biondi et al. 2013;Garbari and Bedini 2014;Carta et al. 2017;Orsenigo et al. 2018). Molecular study revealed two more lineages , molecularly distinct from each other and from all the others, but morphologically indistinguishable from each other and only slightly distinguishable from all the other lineages of M. cantiana.

Taxonomic sample
Four new populations of Monacha cantiana s.l. were considered in our analysis of their molecular and morphological (shell and genitalia structure) variability (Table 1) and compared with the other M. cantiana lineages . The sequences deposited in GenBank were also considered for the molecular analysis. Two other Monacha species were used for molecular comparison (Monacha cartusiana (Müller, 1774)) and for morphological and molecular comparison (M. parumcincta (Rossmässler, 1834)).

Materials examined
New materials examined are listed as follows, when possible: geographic coordinates of locality, locality (country, region, site, municipality and province), collector(s), date, number of specimens with name of collection where materials are kept in parenthesis ( Table 1). The materials are kept in the F. Giusti collection (FGC; Dipartimento di Scienze Fisiche, della Terra e dell'Ambiente, Università di Siena, Italy). The materials used for comparison have already been described (see Pieńkowska et al. 2018: table 1) and is now supplemented with some new nucleotide sequences (Table 2).

DNA extraction, amplification and sequencing
DNA extraction, amplification and sequencing methods are described in detail in our previous paper .

Phylogenetic inference
Two mitochondrial and two nuclear gene fragments were analysed, namely cytochrome c oxidase subunit 1 (COI), 16S ribosomal DNA (16SrDNA), histone 3 (H3) and an internal transcribed spacer of rDNA (ITS2), respectively. All new sequences were deposited in GenBank (Tables 1, 2). The COI, 16SrDNA, H3 and ITS2 sequences obtained from GenBank for comparisons are listed in Table 3.
The sequences were edited by eye using the programme BIOEDIT, version 7.2.6 (Hall 1999). Alignments were performed using CLUSTAL W (Thompson et al. 1994) implemented in MEGA 7 (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 591 bp for COI, 355 positions for 16SrDNA, 315 bp for H3 and 496 positions for ITS2. The sequences were collapsed to haplotypes (COI and 16SrDNA) and to common sequences (H3 and ITS2) 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 16SrDNA) and nuclear (H3 and ITS2) sequences were combined (Table 4) before phylogenetic analysis. Finally, the sequences of COI, 16SrDNA, H3 and ITS2 were combined (Table 4) for Maximum Likelihood (ML) and Bayesian inference (BI). Before doing so, uncertain regions were removed from 16SrDNA alignment with the GBLOCKS 0.91b (Castresana 2000;Talavera and Castresana 2007) using parameters for relaxed selection of blocks. This procedure shortened 16SrDNA sequences from 355 to 275 positions.
The sequences of COI obtained in this study together with other sequences from Gen-Bank were analysed by the genetic distance Neighbour-Joining method (Saitou and Nei 1987) implemented in MEGA7 using the Kimura two-parameter model (K2P) for pairwise distance calculations (Kimura 1980). Maximum Likelihood (ML) analyses were then performed with MEGA 7. Monacha cartusiana and Monacha parumcincta were added as outgroup species in each analysis. For ML analysis of combined sequences, the following best nucleotide substitution models were specified according to the Bayesian Information Criterion (BIC): HKY+G (Hasegawa et al. 1985;Kumar et al. 2016) for COI and 16SrDNA combined sequences of 879 positions (591 COI + 288 16SrDNA), TN92+G   (Ronquist and Huelsenbeck 2003) using the same evolution model as for ML calculation. The GTR substitution model (Nei and Kumar 2000;Kumar et al. 2016), assuming a gamma distributed rate variation (+G) allowing for some sites to be evolutionarily invariable (+I), was identified as the best-fit substitution model using JMODELTEST2 (Darriba et al. 2012). 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. The haplotype network was inferred with NETWORK 5.0.0.1 to reflect all relationships between COI and 16SrDNA haplotypes. During the analysis, a median-joining calculation implemented in NETWORK 5.0.0.1 was used (Bandelt et al. 1999).

Morphological study
Seventy-eight specimens of seven clades (six lineages of M. cantiana s.l.: CAN-1, CAN-2, CAN-3, CAN-4, CAN-5 and CAN-6; one lineage of M. parumcincta) were considered for shell variability (see Table 1 and . Shell variabil- Table 3. GenBank sequences used for comparison in molecular analysis. ity 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, UD umbilicus diameter (see : fig. 1). Seventy-five specimens of seven clades (all lineages of M. cantiana s.l. plus one lineage of M. parumcincta) were analysed for anatomical variability (see Table 1 and . 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, E epiphallus (from base of flagellum to beginning of penial sheath), F flagellum, FO free oviduct, GA genital atrium, OSD ovispermiduct, P penis, V vagina, VA vaginal appendix (also known as appendicula), VAS vaginal appendix basal sac, VD vas deferens. Six anatomical variables (DBC, E, F, P, V, VA) were measured using a calliper under a light microscope (0.01 mm) (see : fig. 2).

Species
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 Z-matrices (shape-related matrices), are described in a previous paper .
Differences between species for each shell and genital character were assessed through box-plots and descriptive statistics. Significance of differences (set at p ≤ 0.01) was obtained using analysis of variance (ANOVA); when the test proved significant, an adjusted posteriori pair-wise comparison between pairs of species was performed using Tukey's honestly significant difference (HSD) test. All variables were log transformed before analysis.

Figure 3.
Bayesian 50% majority-rule consensus tree of the combined data set of COI and 16SrDNA haplotypes, and H3 and ITS2 common sequences (see Table 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 and M. parumcincta combined sequences obtained from GenBank (Table 4). in pairs (Table 5). The clades CAN-5 and CAN-6 differed considerably (12.4-14.3%). The clade CAN-5 differed to a similar degree from CAN-3 and CAN-4 clades (13.3-15.4%). Differences between these two clades (CAN-3 and CAN-4) and the clade CAN-6 were even larger (14.3-16.8%). Both CAN-5 and CAN-6 were also separated by very large genetic distances from all other clades (16.5-21.3%).

Morphological study: shell
The RDA with lineage constraint on the shape and size matrix (Fig. 16) showed that RDA 1 (44%, p < 0.001) separated the groups CAN-1, CAN-2, CAN-3, CAN-4, CAN-5 and CAN-6 from PAR. The preliminary classic PCA revealed size as the first major source of morphological variation, since PC1 (74%) was a positive combination of all variables. On the contrary, RDA 2 (7%, p < 0.01) separated CAN-1, CAN-2 and CAN-3 from CAN-4, CAN-5 and CAN-6 with PAR in intermediate position. In this regard, PC2 (11%) accounted for a contrast between LWmH vs LWaH and PWH variables.  RDA on the shape (Z) matrix (Fig. 17) showed no separation of lineages, confirming that size is a major source of morphological variation. Shape-related PCA indicated that LWfW and LWmW vs SH, LWaH and PWH were the two principal shape determinants on PC1 and AD vs UD on PC2.
Box plots (Fig. 18) proved the poor discriminating value of shell characters in distinguishing lineage pairs. The best discriminant character was UD that distinguished 13 clade pairs according to Tukey's honestly significant difference test, followed by LWmH and LWmW that distinguished seven clade pairs each. The most recognizable pairs were CAN-1 vs PAR, CAN-3 vs PAR, CAN-6 vs PAR, CAN-2 vs PAR and CAN-5 vs PAR (12, 11, 10, 8 and 7 significant characters, respectively). Five significant shell characters distinguished CAN-3 vs CAN-4, four CAN-4 vs CAN-6, two CAN-1 vs CAN-4, CAN-1 vs CAN-5 or CAN-3 vs CAN-5 and only one CAN-1 vs   (Table 6).

Morphological study: anatomy
The bodies (generally pinkish or yellowish white) and mantle (with sparse brown or blackish spots near the mantle border or on the lung surface, a larger one close to the pneumostomal opening) of CAN-5 and CAN-6 are very similar to those of the other lineages of M. cantiana s.l. and M. parumcincta studied so far ). The    internally jagged and with a sort of solid pilaster on one side in M. cantiana s.l.; central canal not connected to external wall, internally smooth or slightly jagged and almost completely filled by large invagination in M. parumcincta).
RDA with lineage constraint on the shape and size matrix (Fig. 42) showed that RDA 1 (36%, p < 0.001) separated the M. cantiana s.l. (CAN-1, CAN-2, CAN-3, CAN-4, CAN-5 and CAN-6) from PAR. The preliminary classic PCA revealed size as the first major source of morphological variation, since PC1 (54%) was a positive combination of all variables. On the contrary, RDA 2 (12%, p < 0.001) separated the group CAN-1, CAN-2, CAN-3, CAN-4 and PAR from the group CAN-5 and CAN-6. In that regard, PC2 (17%) accounted for a contrast between P and DBC vs F.  RDA with species constraint on the shape (Z) matrix (Fig. 43) showed that RDA 1 (28%, p < 0.001) separated the group CAN-1, CAN-2 and CAN-3 from the group CAN-5 and CAN-6 with PAR and CAN-4 in intermediate position and that RDA 2 (11%, p < 0.001) separated PAR from CAN-4 with the large group CAN-1, CAN-2, CAN-3, CAN-5 and CAN-6 in intermediate position. Shape-related PCA indicated that P, E and DBC vs VA and F were the two principal shape determinants on PC1 and DBC and VA vs E, V and F on PC2. In the latter case, removing the size effect altered the overall relationship patterns.
Box plots (Fig. 44) for anatomical characters showed that F and VA have the best discriminating value (they distinguished 11 and 8 clade pairs, respectively, according to Tukey's honestly significant difference test), followed by E and V (five and four  (Table 6).  found that M. cantiana, as usually conceived, actually consists of four distinct lineages (CAN-1, CAN-2, CAN-3 and CAN-4). Examination of a group of four additional populations from the Apuan Alps revealed two more lineages . From a molecular point of view, they are quite distinct from each other and from all the others but from a morphological point of view they are indistinguishable from each other and only slightly distinguishable from the others.

Discussion
Our present results confirm that lineages CAN-1, CAN-2 and CAN-3 can be distinguished by analysis of mitochondrial gene (COI and 16SrDNA) sequences (Figs 1, 4, 5) but not by nuclear gene (H3 and ITS2) sequences (Fig. 2). On the other hand, analysis of both nucleotide sequences (of mitochondrial and nuclear genes) showed that the CAN-4, CAN-5 and CAN-6 lineages are distinct from all the others (Figs 1-3). Moreover, these gene sequences clearly separated M. cantiana lineages from M. parumcincta.
Based on their studies of lepidopteran relationships, Hebert et al. (2003a, b) suggested that nucleotide sequences of the mitochondrial COI gene could be a universal tool for species distinction. This so called "barcode method" has since been widely used (Tautz et al. 2003;Hebert et al. 2004Hebert et al. , 2013Hajibabaei et al. 2007;Packer et al. 2009;Goldstein and Desalle 2011;Čandek and Kuntner 2015;Dabert et al. 2018;Yang et al. 2018, but see e.g.: Moritz and Cicero 2004;Taylor and Harris 2012). It has also been used to solve taxonomic problems in different gastropod families (Hershler et al. 2003;Remigio and Hebert 2003;Rundell et al. 2004;Elejalde et al. 2008;Duda et al. 2011;Delicado et al. 2012;Breugelmans et al. 2013;Proćków et al. 2013Proćków et al. , 2014. However, a 3% threshold was established arbitrarily by Hebert et al. (2003a, b) as a marker of species distinction, and in several stylommatophoran families it proves to be much higher (Davison et al. 2009;Hausdorf 2010, 2012;Scheel and Hausdorf 2012). Moreover, we have always stressed (Pieńkowska et al. 2015 that molecular features alone are insufficient to define species but need to be supported by anatomical features. In light of the above, we underline that the interspecific genetic distances in COI sequences between both, CAN-5 and CAN-6, and all other lineages of M. cantiana s.l. (CAN-5 vs CAN-1/CAN-2/CAN-3/CAN-4 -13.3-18.2%, CAN-6 vs CAN-1/ CAN-2/CAN-3/CAN-4 -14.3-19.2%; Table 5) are an order of magnitude greater than Hebert's 3% threshold (Hebert et al. 2003a, b). It is also an order of magnitude greater than intraspecific divergence ("barcode gap", see Hebert et al. 2004;Čandek and Kuntner 2015) within CAN-5 and CAN-6 lineages, 1.3% and 1.6%, respectively ( Table 5). The analysis of mitochondrial COI and 16SrDNA sequences (Figs 1, 4, 5) are supported by the results of nuclear ITS2 and H3 sequences (Fig. 2). This suggests that CAN-5 and CAN-6 lineages taken together create a taxon separate from the other lineages of M. cantiana s.l. Despite CAN-5 differs from CAN-6 at a similarly high level  there are no morphological differences between specimens of both lineages. The speciation of CAN-5 and CAN-6 lineages therefore seems to emerge more promptly in molecular (mitochondrial gene sequences) than in morphological (shell, genitalia) features, probably because of a rapidly evolving mitochon-drial genome (Thomaz et al. 1996;Remigio and Hebert 2003). As mentioned above, molecular data alone cannot be used to distinguish species. It must be supported by morphological features of shells and/or genital anatomy before any decision is made about taxonomy or nomenclature.
Statistical analysis of 12 shell and six anatomical characters showed that CAN-5 and CAN-6 cannot be distinguished from each other by morphology (no character shows statistically significant differences according to Tukey's honestly significant difference test). They are only marginally distinct from CAN-1, CAN-2, CAN-3 and CAN-4, but clearly distinct from M. parumcincta, used for comparison: two or three characters distinguish the group CAN-5 plus CAN-6 from CAN-1, CAN-2 and CAN-3; one character distinguishes CAN-5 from CAN-4; five characters distinguish CAN-6 from CAN-4; 11-14 characters distinguish the group CAN-5 plus CAN-6 from PAR. It is possible that the small sample available for lineages CAN-4 and CAN-6 (one population for each) biased comparison of these two lineages. The best discriminant characters separating the group CAN-5 plus CAN-6 from all the other lineages are umbilicus diameter (UD) and flagellum length (F). In both cases the lineages CAN-5 and CAN-6 have the highest values (Table 7).
As in the case of other lineages, the greatest bias of morphological analysis was the small sample available for lineages CAN-2, CAN-3, CAN-4 and CAN-6, which prevented a realistic account of their variability. As far as we know, this newly recognised group only occurs in the Apuan Alps and consists of two differentiated lineages (CAN-5 and CAN-6). Although examination of additional populations is desirable, intra-Apuan differentiation is also known for other organisms such as plants (Bedini et al. 2011) and animals (Zinetti et al. 2013).
Six available names have been introduced for Monacha cantiana s.l. from northwestern Tuscany (see Appendix 1). The oldest, Helix anconae, was established by Issel (1872) for specimens reported from a wide area extending northward to Arenzano in Liguria and southward to island of Elba and the Maremma of Tuscany. However, all the localities quoted are in coastal and lowland Liguria and Tuscany, while the populations including the group CAN-5 plus CAN-6 are from mountain sites. This would exclude a relationship of this nominal taxon with these lineages. F mean ± S.D. 8.5 ± 1.5 7.6 ± 1.0 8.0 ± 1.2 10.2 ± 1.1 14.9 ± 2.5 14.9 ± 1.7 6.9 ± 1.0 range 5.3-12.2 6.2-8.9 6.2-10.0 8.9-11.4 10.8-18.7 13.6-17.8 5.4-8.2 number of specimens 23 7 9 5 15 5 11 All dimensions in mm.
All the other names were established by Mabille (1881) and De Stefani (1883-1888 for specimens collected in the Apuan Alps. Syntypes of the three nominal taxa introduced by Mabille (1881) are in Bourguignat's collection at the Muséum d'histoire naturelle, Genève (Switzerland) (Figs 45-47). Syntypes of the two nominal taxa established by De Stefani (1883-1888 are not known and probably lost. Umbilicus diameter of the shells of the syntypes of Mabille's species and the specimens illustrated by De Stefani is consistent for at least five of these nominal taxa with that of Monacha of the group CAN-5 plus CAN-6 (Helix sobara Mabille, 1881, Helix ardesa Mabille, 1881, Helix apuanica Mabille, 1881, Helix carfaniensis De Stefani, 1883and Helix spallanzanii De Stefani, 1884). Mabille's three nominal taxa have precedence over those of De Stefani, and because the former were published simultaneously in the same paper, their relative precedence can only be determined by the first revisor (ICZN 1999: Art. 24). Although all three match Monacha of the group CAN-5 plus CAN-6, the best correspondence is with Helix sobara. Nevertheless, the availability of these names for the lineages CAN-5 and CAN-6 is somewhat difficult and not immediate. These nominal taxa were only established on shell characters, but no shell character shows statistically significant differences between CAN-5 and CAN-6. Their relationships could only therefore be established by molecular study of topotypes, but unfortunately Mabille (1881) did not quote any precise collection site. In some cases, the identity and relationships of extinct taxa have been addressed and clarified through study of ancient DNA from dried tissue (e.g. Villanea et al. 2016;Vogler et al. 2016). Unfortunately, this approach is not applicable for Mabille's syntypes because they consist only of shells devoid of any dried tissue. Thus, the case can only be solved by appeal to article 75.5 of the Code (ICZN 1999).

Helix spallanzanii
Type locality: Apuan Alps, Vagli. De Stefani (1884: 231) stated that the type is from Apuan Alps and depicted a shell from Vagli.