Specimens from ten French and two Dutch populations of Monacha cantiana were considered for analysis of the variability of their molecular and morphological (shell and genitalia) features (Table 1, Fig. 1). Specimens from four new British and one new Italian population were used for comparative molecular analysis with other populations of M. cantiana s.l. (Table 1, Fig. 1). Sequences deposited in GenBank for M. cantiana s.l. from other populations (Manganelli et al. 2005; Duda et al. 2011; Kruckenhauser et al. 2014; Cadahia et al. 2014; Pieńkowska et al. 2015, 2018b, 2019a, 2019b; Razkin et al. 2015; Neiber and Hausdorf 2017; Čejka et al. 2020, 2022) and three other Monacha species (M. cartusiana: Pieńkowska et al. 2015, 2022; Neiber and Hausdorf 2017; Caro et al. 2019; Čejka et al. 2020; M. pantanellii (De Stefani, 1879): Pieńkowska et al. 2020; M. parumcincta (Rossmässler, 1834): Pieńkowska et al. 2018b) were also selected for molecular analysis (Suppl. materials 1–4) and supplemented with several new sequences of mitochondrial (16SrDNA) and nuclear (ITS2 flanked with 5.8SrDNA and 28SrDNA) genes (Table 1). Sequences of Trochulus hispidus (Linnaeus, 1758) deposited in GenBank by Neiber et al. (2017), Neiber and Hausdorf (2017), Caro et al. (2019) and Proćków et al. (2021) were used as an outgroup to construct phylogenetic trees (Suppl. materials 1–4). The localities for reference populations of M. cantiana s.l. CAN-1 – CAN-6, M. pantanellii, M. cartusiana, and M. parumcincta were shown on maps published in our previous papers (Pieńkowska et al. 2018b: fig. 63, 2020: fig. 1).

The material examined originated from the populations listed in Table 1 with the following data: geographic coordinates, country and region, short description of collection site, name of collector, date, number of specimens studied and the collection where the material is stored (in brackets). The origin of the material used for comparison has been described in previous publications (Pieńkowska et al. 2015: appendix 1; Pieńkowska et al. 2018b, 2019a, 2019b, 2020, 2022: table 1).

Sixty-six specimens of the six lineages of M. cantiana s.l. (CAN-1, CAN-2, CAN-3, CAN-4, CAN-5, and CAN-6) (Pieńkowska et al. 2018b, 2019a) and five specimens suitable for morphological analysis of the French populations were considered for shell variability (Table 1). Twelve shell variables were measured to the nearest 0.1 mm using ADOBE PHOTOSHOP 7.0.1 on digital images of standard apertural and umbilical views taken with a Canon EF 100 mm 1:2.8 L IS USM macro lens mounted on a Canon F6 camera (see also Pieńkowska et al. 2018b: fig. 1):

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.

Sixty-four specimens of the six lineages of M. cantiana s.l. (CAN-1, CAN-2, CAN-3, CAN-4, CAN-5 and CAN-6) (Pieńkowska et al. 2018b, 2019a) and seven adult specimens of the French populations were analysed for anatomical variability (Table 1). Snail bodies were dissected under a light microscope (Wild M5A or Zeiss SteREO Lumar V12). Anatomical details were drawn using a Wild camera lucida. Abbreviations/acronyms are as follows (see also Pieńkowska et al. 2018b: fig. 2):

Figure 2. 

Shells of Monacha cantiana from France. Specimen Esc1 from Oise, Escales-Saint-Pierre (A), specimen Ble1 from Pas-de-Calais, Blecquenecques n. Marquise (B), specimen Pie1 from Seine-Maritime, Pierrepont (C) and specimen Bet1 from Seine-Maritime, Béthencourt n. Grandcourt (D).

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,

GAR genital atrium retractor,

P penis,

PP penial papilla,

SOD spermoviduct,

V vagina,

VA vaginal appendix (also known as appendicula),

VAS vaginal appendix 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) (Pieńkowska et al. 2018b: 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 Z-matrices (shape-related matrices), are described in our previous papers (Pieńkowska et al. 2018b, 2019a).

We used 95% confidence interval ellipses to evaluate the uncertainty of the estimate of the population mean (centroid) of the data sample. The function ordiellipse with standard errors in the package vegan (Oksanen et al. 2022) was used. Convex hulls (function ordihull in vegan) were used to visually enclose the individuals forming each clade as a measure of data spread. All analyses were performed with RStudio (R version 4.2.1; R Core Team 2021).

Eighty-eight specimens representing 26 populations of the four lineages of M. cantiana s.l. (CAN-1, CAN-2, CAN-3, and CAN-4; Pieńkowska et al. 2018b, 2019a) were used for molecular analysis (Table 1). Molecular methods including DNA extraction, amplification and sequencing are described 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 (16SrDNA), histone 3 (H3) and an internal transcribed spacer 2 of rDNA (ITS2) flanked by the 3’end of 5.8SrDNA and the 5’end of 28SrDNA, respectively. Sequences were edited by eye using BioEdit, v. 7.0.6 (Hall 1999; BioEdit 2017) and aligned using ClustalW, implemented in BioEdit (Thompson et al. 1994). Fragments of COI were amplified using two pairs of primers: F01/R04 (Dabert et al. 2010) or bcsmF1/bcsmR1 (Proćków et al. 2013). Fragments of 16SrDNA were amplified using 16Scs1/16Scs2 primers (Chiba 1999). Sequences containing the 3’end of 5.8SrDNA, complete sequence of ITS2 and 5’end of 28SrDNA were amplified using two sets of primers: LSU1/LSU3 (Wade and Mordan 2000) and NEWS2/ITS2-RIXO (Almeyda-Artigas et al. 2000). Products of the two PCR reactions were aligned and used to assemble single sequences. Fragments of H3 gene were amplified using the primers H3F and H3R (Colgan et al. 1998). The protein coding sequences were aligned according to the translated amino acid sequences. The ends of all sequences were trimmed. After trimming, the lengths of sequences were 615 bp for COI, 804–821 bp for 16SrDNA, 303 bp for H3, and 749–753 bp for ITS2 flanked by the 3’end of 5.8SrDNA and 5’end of 28SrDNA (including 45 bp 5.8SrDNA + 489–493 bp ITS2 + 215 bp 28SrDNA). The borders of ITS2 sequence were searched using ITS2-Database (http://its2.bioapps.biozentrum.uni-wuerzburg.de) (Eddy 1998; Koetschan et al. 2010). The sequences were collapsed to haplotypes using the programme ALTER (Alignment Transformation EnviRonment) (Glez-Peña et al. 2010). The following alignments were made for phylogenetic inference: 591 bp long for COI, 292 or 809 positions long for 16SrDNA, and 775 positions long for ITS2 flanked by the 3’end of 5.8SrDNA and 5’end of 28SrDNA. Finally, the sequences of COI, 16SrDNA, ITS2, and H3 were concatenated. Three sets of concatenated sequences were created: 1) COI16S of 1444 positions in length (615 COI + 829 16SrDNA); 2) H3ITS2 of 1054 positions in length (279 H3 + 775 ITS2 with flanks); 3) CS of 2498 positions in length (615 COI + 829 16SrDNA + 279 H3 + 775 ITS2 with flanks).

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

To infer the phylogenetic relationships the following programmes were used: MEGA7 (Hasegawa et al. 1985; Nei and Kumar 2000; Kumar et al. 2016), IQ-Tree (http://iqtree.cibiv.univie.ac.at/) (Trifinopoulos et al. 2016), RAxML v1.0.0 (Stamatakis 2014) and MrBayes 3.2.6 (Ronquist et al. 2012). For phylogenetic inference Neighbour-Joining, Maximum-Likelihood and Bayesian Inference methods were used.

For each alignment file, best nucleotide substitution models were specified according to the Bayesian Information Criterion (BIC) (see captions to figures). Phylogenetic analyses performed with IQ-Tree, RAxML and MrBayes for three sets of concatenated sequences were done dividing the data set into 2 or 4 partitions: (1) COI, (2) 16SrDNA or (1) COI, (2) 16SrDNA, (3) H3, (4) 5.8SrDNA + ITS2 + 28SrDNA. Best substitution models were inferred according to the Bayesian Information Criterion (BIC) for each of the partitions by MODELFINDER (Kalyaanamoorthy et al. 2017) implemented in IQ-TREE. Bayesian analysis were conducted with four Monte Carlo Markov chains running for 1 million generations, sampling every 100 generations (the first 25% of trees were discarded as ‘burn-in’).

The robustness of the NJ and ML trees generated by MEGA7 were assessed by bootstrap analysis with 1000 replicates (Felsenstein 1985). ML trees built by RAxML were tested by bootstrap analysis with 100 replicates. ML trees obtained with IQ-Tree were constructed under 1000 ultrafast bootstrap replicates (Minh et al. 2013). Finally, BI trees were supported by posterior probability (PP) values. Bootstrap support values from NJ and ML analysis as well as posterior probability (PP) values obtained on 50% majority rule consensus Bayesian tree were mapped onto the ML tree obtained by MEGA7. All the resulting trees were rooted with Trochulus hispidus sequences obtained from GenBank.

Shells of French specimens of M. cantiana (Fig. 2A–D) are globose-subglobose in shape, variable in size and usually whitish or pale yellowish in colour, with slightly descending, roundish to oval aperture, similar to those of the other populations of the lineage CAN-1 (Pieńkowska et al. 2018b: figs 8–11).

RDA with French specimens and “lineage” constraint on the shape and size matrix (Fig. 3B, C) showed that RDA 1 (22.2%, P < 0.01) separated CAN-6 from CAN-4, with CAN-5 and the large group CAN-1, CAN-2, CAN-3, and FRA in intermediate position, as confirmed by 95% confidence interval ellipses (Fig. 3B). The convex hull measure of data spread showed considerable overlap of some clusters. In both cases, FRA specimens fell within CAN-1 variability (Fig. 3B). The preliminary classic PCA showed that size was the first major source of morphological variation, since PC1 (69%) was a positive combination of all variables (Fig. 3A). On the contrary, RDA 2 was not significant (p > 0.05) and accounted for little morphological variation (2.6%). PC2 (15%) mostly reflected a contrast between LWaH and PWH versus LWmH and UD.

Figure 3. 

Analysis of French specimens with “lineage” constraint on the original matrix (A–C) and Z-matrix (shape-related) (D–F) of selected shell sections. Principal component analysis (PCA) (A, D) and redundancy analysis (RDA) with groups shown as ellipses representing 95% confidence intervals with standard errors (B, E) and as convex hull polygons (C, F).

RDA on the shape (Z) matrix (Fig. 3E, F) showed that RDA 1 (34%, P < 0.001) clearly separated CAN-5 and CAN-6 from the group CAN-1, CAN-2, CAN-3, CAN-4, and FRA, as confirmed by the 95% confidence interval ellipses (Fig. 3E) and the convex hulls (Fig. 3F). On the contrary, the RDA2 axis was not significant (P > 0.05), reflecting little morphological variation (5.3%). Shape-related PCA indicated that SH, LWaH and PWH vs LWmW, SD, LWfW, LWmH, and UD were the principal shape determinants on PC1, and AH and AW vs PWH, LWmH, and UD on PC2 (Fig. 3D).

French specimens of M. cantiana have distal genitalia (Figs 4–6) resembling the other populations assigned to CAN-1, which are in turn similar to those of the populations belonging to the CAN-2, CAN-3 and CAN-4 lineages (Pieńkowska et al. 2018b: figs 20–30).

Figure 4. 

Distal genitalia of Monacha cantiana from France. Specimen Bet1 from Seine-Maritime, Béthencourt n. Grandcourt (A–C) and specimen Ble1 from Pas-de-Calais, Blecquenecques n. Marquise (D–F). Distal genitalia (A, D), transverse sections of medial epiphallus (B, E) and apical penial papilla (C, F).

Figure 5. 

Distal genitalia genitalia of Monacha cantiana from France. Specimen Esc1 from Oise, Escales-Saint-Pierre (A–C) and specimen Epa1 from Somme, Épagne-Épagnette, roadside (D). Distal genitalia (A), transverse sections of medial epiphallus (B), apical penial papilla (C) and internal structure of distal genitalia (D).

Figure 6. 

Distal genitalia of of Monacha cantiana from France. Specimen Pie1 from Seine-Maritime, Pierrepont. Distal genitalia (A), transverse sections of medial epiphallus (B), apical penial papilla (C) and internal structure of distal genitalia (D).

RDA with French specimens and “lineage” constraint on the shape and size matrix (Fig. 7B, C) showed that RDA 1 (24.3%, P < 0.001) separated CAN-2 and CAN-6 from FRA and CAN-5, with CAN-1, CAN-3, CAN-4 in intermediate position, as confirmed by 95% confidence interval ellipses (Fig. 7B). The preliminary classic PCA showed that size was the first major source of morphological variation, since PC1 (48.3%) was a positive combination of all variables (Fig. 7A). On the other hand, RDA 2 (21.7%, P < 0.001) clearly separated the group CAN-1, CAN-2, CAN-3, CAN-4 and FRA from CAN-5 and CAN-6. PC2 (17.9%) reflected a contrast between P, VA and DBC vs F and V. Differences between clusters were confirmed visually by 95% confidence interval ellipses (Fig. 7B) and convex hulls (Fig. 7C).

Figure 7. 

Analysis of French specimens with “lineage” constraint on the original matrix (A–C) and Z-matrix (shape-related) (D–F) of selected genital sections. Principal component analysis (PCA) (A, D) and redundancy analysis (RDA) with groups shown as ellipses representing 95% confidence intervals with standard errors (B, E) and as convex hull polygons (C, F).

RDA on the shape (Z) matrix (Figs 7E, F) showed that RDA 1 (33.7%, P < 0.001) separated the 95% confidence interval ellipses of CAN-5, CAN-6 and CAN-4 from the large group CAN-1, CAN-2, CAN-3, and FRA; RDA 2 (8%, P < 0.001) separated CAN-5 and the group CAN-1, CAN-2, CAN-3, FRA from CAN-6 and CAN-4 (Fig. 7E). Convex hulls showed some overlaps, especially in the data spread of CAN-1 (Fig. 7F). Shape-related PCA indicated that P and E vs VA and F were the two principal shape determinants on PC1 and DBC and VA vs V and F on PC2 (Fig. 7D).

Although sequences of all the genes analysed (COI, 16SrDNA, H3, and ITS2 with 5.8SrDNA and 28SrDNA) were not obtained from all 88 specimens (Table 1), as a result of molecular analysis, 272 new sequences were deposited in GenBank. These were 56 new sequences of COI: OR918493–OR918548, 77 of 16SrDNA (long): OR918363–OR918439, 75 of H3: OR939858–OR939932 and 64 of ITS2 (with flanking fragments of 5.8SrDNA and 28SrDNA): OR917347–OR917410 (Table 1). Eleven haplotypes of the COI gene were identified (COI 1 – COI 11), 32 of 16SrDNA (16S 1 – 16S 32), 10 of H3 (H3 1 – H3 10), and 34 of ITS2 with flanking sequences (ITS2 1 – ITS2 34) (Table 1). These haplotypes were used for phylogenetic analysis based on single gene sequences and concatenated mitochondrial and nuclear gene data sets of sequences.

The phylogenetic analysis of COI sequences obtained from the specimens and comparative sequences derived from GenBank is shown in Fig. 8. The results are consistent with previously published findings (Pieńkowska et al. 2018b, 2019a, 2019b, 2020, 2022), distinguishing six lineages (CAN-1 – CAN-6) in M. cantiana s.l. that clustered separately from COI sequences of other species including M. parumcincta, M. pantanellii and M. cartusiana. The new COI sequences (haplotypes 1–10) from France, the Netherlands and England clustered in the CAN-1 lineage. Only the COI 11 haplotype obtained from two specimens of the Italian population from Sasso di Simone (population no. 23 in Table 1) grouped with the CAN-2 lineage.

Figure 8. 

Maximum Likelihood (ML) tree of COI haplotypes of Monacha cantiana. New COI sequences of M. cantiana (Table 1) were compared with COI sequences of M. cantiana s.l., M. parumcincta, M. pantanellii and M. cartusiana obtained from GenBank (Suppl. material 1). Sequences were cut to 591 bp. HKY+G+I was the best nucleotide substitution model according to the Bayesian Information Criterion (BIC). The tree was rooted with Trochulus hispidus sequences obtained from GenBank (Suppl. material 1).

K2P genetic distances (Table 2) showed small genetic differentiation between COI sequences of particular CAN-1 populations (infra-populational distances ranged from 0.2% in Dutch populations to 1.1% in French populations). The K2P distances between these populations were also small (in the range 0.5–1.2%). The K2P distances between French, Dutch, English and Italian populations of CAN-1 and CAN-2 were also small (in the range 3.5–4.1%) while the distance separating the CAN-1 populations from the CAN-3 and CAN-4 populations was much larger (in the range 18.0–18.8%). In turn, the distance separating the CAN-3 and CAN-4 populations was 5.6–6.1%.

Results similar to those of COI analysis were obtained for other single gene analyses (Suppl. materials 8, 9 for 16SrDNA, Suppl. material 10 for the ITS2 gene with flanking 5.8S and 28S gene fragments). Note that the newly obtained 16SrDNA sequences in Suppl. material 8 were trimmed to 292 positions in alignment length because GenBank lacks the reference long 16SrDNA sequences of the 809 positions used to construct the tree in Suppl. material 9. Analysis of newly obtained longer sequences (i.e. ITS2 flanked by 5.8SrDNA and 28SrDNA gene fragments) (ITS2 1 – ITS2 34 haplotypes) and the only comparable sequence of Neiber and Hausdorf (2017) showed that this gene did not differentiate the CAN-1, CAN-2 and CAN-3 lineages. Similar results were obtained previously using ITS2 gene sequences without flanking fragments of 5.8SrDNA and 28SrDNA (Pieńkowska et al. 2018b: fig. 64). Only in the case of sequences assigned to the CAN-4 lineage were they distinct from CAN-1, CAN-2 and CAN-3, as shown in Pieńkowska et al. (2018b: fig. 64).

The phylogenetic tree for concatenated sequences were similar in ML analyses obtained with different software. The tree for mitochondrial gene sequences (COI+16SrDNA) in Fig. 9 shows that the sequences obtained from specimens of the French, Dutch, and English populations (see also Suppl. material 5) grouped with the reference sequences for CAN-1. In a tree of concatenated nuclear genes (Fig. 10: H3+ITS2 with flanks), the sequences from the French populations grouped with CAN-1, CAN-2, and CAN-3 lineages, only sequences of the CAN-4 lineage being distinguished. However, note that the bootstrap and posterior probability values weakly supported the results of the concatenated H3+ITS2 gene sequences. The tree for the concatenated sequences of all the genes analysed in this paper (Fig. 11, see also Suppl. material 7) showed that concatenated sequences CS 1–CS 25 from northern French populations clustered together with CS 26–CS 34 and CS 35–CS 38 sequences obtained from English and Italian specimens, respectively. They all belonged to the CAN-1 lineage. The CAN-1, CAN-2, CAN-3, and CAN-4 lineages grouped separately.

Figure 9. 

Maximum Likelihood (ML) tree of concatenated COI and 16SrDNA haplotypes of Monacha cantiana. New COI and 16SrDNA sequences of M. cantiana (Table 1, Suppl. material 5) were compared with concatenated COI and 16SrDNA sequences of M. cantiana s.l. and M. cartusiana obtained from GenBank (Suppl. materials 1, 2, 5). Length of sequences was 1444 positions (615 of COI + 829 of 16SrDNA). The Bayesian Information Criterion (BIC) specified T92+G+I the best nucleotide substitution model in MEGA7, or HKY+F+G4 for COI and TIM2+F+G4 for 16SrDNA partition in IQ-Tree, RAxML and MrBayes. Numbers next to main branches indicate (left to right): bootstrap supports above 50% calculated by NJ-MEGA7 (Saitou and Nei 1987), ML-MEGA7 (Kumar et al. 2016), IQ-Tree (Trifinopoulos et al. 2016), RAxML (Stamatakis 2014), and posterior probabilities by BI (Ronquist et al. 2012). The tree was rooted with Trochulus hispidus concatenated sequences obtained from GenBank (Suppl. material 5).

Figure 10. 

Maximum Likelihood (ML) tree of concatenated H3 and ITS2 (flanked with 5.8S and 28SrDNA) haplotypes of Monacha cantiana. New H3 and ITS2 sequences of M. cantiana (Table 1) were compared with concatenated H3 and ITS2 sequences of M. cantiana s.l. obtained from GenBank (Suppl. materials 3, 4). Length of sequences was 1054 positions (279 of H3 + 775 of ITS2). The Bayesian Information Criterion (BIC) specified T92+G+I the best nucleotide substitution model in MEGA7, or K2P+I for H3 and K3P+I for ITS2 partition in IQ-Tree, RAxML, and MrBayes. Numbers next to main branches indicate (left to right): bootstrap supports above 50% calculated by NJ-MEGA7 (Saitou and Nei 1987), ML-MEGA7 (Kumar et al. 2016), IQ-Tree (Trifinopoulos et al. 2016), RAxML (Stamatakis 2014) and posterior probabilities by BI (Ronquist et al. 2012). The tree was rooted with Trochulus hispidus concatenated sequences obtained from GenBank (Suppl. material 6).

Figure 11. 

Maximum Likelihood (ML) tree of concatenated COI, 16SrDNA, H3, and ITS2 (flanked with 5.8S and 28SrDNA) haplotypes of Monacha cantiana. COI, 16SrDNA, H3, and ITS2 sequences of M. cantiana were compared with sequences of M. cantiana s.l. and M. cartusiana obtained from GenBank (Suppl. materials 1–4, 7). Length of sequences was 2498 positions (615 of COI, 829 of 16SrDNA, 279 of H3, and 775 of ITS2). Bayesian Information Criterion (BIC) specified GTR+G+I the best nucleotide substitution model in MEGA7, or HKY+F+G4 for COI, TIM2+F+I for 16SrDNA, TIM3e+I+G4 for H3, and K3P+I+G4 for ITS2 partition in IQ-Tree, RAxML, and MrBayes. Numbers next to main branches indicate (left to right): bootstrap support above 50% calculated by NJ-MEGA7 (Saitou and Nei 1987), ML-MEGA7 (Kumar et al. 2016), IQ-Tree (Trifinopoulos et al. 2016), RAxML (Stamatakis 2014), and posterior probabilities by BI (Ronquist et al. 2012). The tree was rooted with Trochulus hispidus concatenated sequences obtained from GenBank (Suppl. material 7).

The authors have declared that no competing interests exist.

No ethical statement was reported.

DB and GM were funded under the National Recovery and Resilience Plan (NRRP), Mission 4 Component 2 Investment 1.4 – Call for tender No. 3138 of 16 December 2021, rectified by Decree n.3175 of 18 December 2021 of Italian Ministry of University and Research funded by the European Union – NextGenerationEU; Award Number: Project code CN_00000033, Concession Decree No. 1034 of 17 June 2022 adopted by the Italian Ministry of University and Research, CUPB63C22000650007 Project title “National Biodiversity Future Center - NBFC”. AL received support from Adam Mickiewicz University, Poznań (Poland), funded project 526000/REZ_PROJEKT “Powrót do badań”.

Conceptualization: AL, FG and GM; Methodology, Formal analysis, Investigation, Data Curation on shell and genitalia: FG, DB and GM; Methodology, Formal analysis, Investigation, Data Curation on molecular data: AL, JRP, KS and MP; Writing - Original draft & Writing - Review and Editing: AL, FG and GM; Supervision: FG, AL and GM; Funding Acquisition: AL and GM.

Joanna R. Pieńkowska https://orcid.org/0000-0003-0372-121X

Giuseppe Manganelli https://orcid.org/0000-0002-8453-280X

Małgorzata Proćków https://orcid.org/0000-0003-2240-7306

Debora Barbato https://orcid.org/0000-0003-1105-1711

Katarzyna Sosnowska http://orcid.org/0000-0002-7506-4231

Folco Giusti https://orcid.org/0000-0001-8722-4653

Andrzej Lesicki https://orcid.org/0000-0002-1924-1934

All of the data that support the findings of this study are available in the main text or Supplementary Information.