Unveiling cryptic diversity among Müllerian co-mimics: insights from the Western Palaearctic Syntomis moths (Lepidoptera: Erebidae: Arctiinae)

Accurate species delimitation is of primary importance in biodiversity assessments and in reconstructing patterns and processes in the diversification of life. However, the discovery of cryptic species in virtually all taxonomic groups unveiled significant gaps in our knowledge of biodiversity. Mimicry complexes are good candidates to source for cryptic species. Indeed, members of mimicry complexes undergo selective pressures on their habitus, which results in strong resemblance even between distantly related species. In this study, we used a multi-locus genetic approach to investigate the presence of cryptic diversity within a group of mimetic day-flying moths whose systematics has long been controversial, the Euro-Anatolian Syntomis. Results showed incongruence between species boundaries and the currently accepted taxonomy of this group. Both mitochondrial and nuclear markers indicate the presence of four, well-distinct genetic lineages. The genetic distance and time of divergence between the Balkan and Italian populations of S. marjana are the same as those found between S. phegea and S. ragazzii, the last two being well-distinct, broadly sympatrically occurring species. The divergence between the two lineages of S. marjana dates back to the Early Pleistocene, which coincided with substantial changes in climatic conditions and vegetation cover in Southern Europe that have likely induced geographic and ecological vicariance. Syntomis populations belonging to the taxa kruegeri (s. str.), albionica and quercii are now considered a separate species from marjana s. str. and are thus distinguished as Syntomis quercii Verity, 1914, bona sp., stat. nov. Our results show that the species richness of mimicry complexes inhabiting temperate regions might still be severely underestimated.


Introduction
Accurate species delimitation is of primary importance for understanding tempo and modes of species diversification (Wiens 2007). Historically, new species have mainly been recognized using morphological trait variation (Wiens and Servedio 2000). However, species can also evolve in absence of conspicuous morphological trait divergence, as shown by the discovery of the so-called cryptic species (Bickford et al. 2007). These had originally been spotted based on subtle ecological or behavioural features (Mayr 1963), though the use of genetic markers has nowadays boosted their discovery Knowlton 1993;Beheregaray and Caccone 2007;Carstens et al. 2013). The recognition of cryptic species virtually across all taxonomic groups (Pfenninger and Schwenk 2007), including many well-studied taxa (Roca et al. 2001;Fennessy et al. 2016), unveiled significant gaps in our knowledge of biological diversity (Beheregaray and Caccone 2007;Struck et al. 2018). These gaps, in turn, constrain our understanding of the mechanisms involved in biological diversification and in the establishment of interactions at the community level, limiting the deployment of measures in many areas of biodiversity management, from conservation biology to pest or disease control (Beheregaray and Caccone 2007;Bálint et al. 2011;Robuchon et al. 2019).
Although cryptic species are widespread throughout the tree of life, they are over-reported in some groups, such as freshwater fishes, deep-sea clams, polychaetes, frogs, mites, parasitic insects and nematodes (reviewed in Pérez- Ponce de León and Poulin 2016). Cryptic species seem to be commoner in animals using non-visual mating signals (e.g. many frogs) or under selection promoting morphological stasis or convergent evolution (e.g. parasites) (Bickford et al. 2007;Struck et al. 2018). They are also common in cases of recent and/or rapid speciation when short divergence time did not allow the accumulation of detectable phenotypic differences (Reidenbach et al. 2012;Gustafsson et al. 2014). However, numerous taxa characterized by life-history traits likely encouraging cryptic diversification have been poorly investigated, and many of them are still analysed with a traditional morphological approach (Bickford et al. 2007;Struck et al. 2018). As a consequence, a clear picture of the distribution of cryptic species across the tree of life is still lacking.
Mimicry complexes are good sources of cryptic species. In particular, cryptic species are expected among closely related species involved in a same Müllerian mimicry ring (Ruxton et al. 2004, cf. Boppré et al. 2017, because phenotypic divergence is constrained by stabilizing selection on shared warning signals (but see Lawrence et al. 2019). Interestingly, it was shown that look-alikeness among co-mimics could also be achieved via introgressive hybridization of patterndetermining genes from related species (Giraldo et al. 2008;Pardo-Diaz et al. 2012;Jiggins 2017). The imbalance between species divergence and phenotypic divergence is a prerequisite for the evolution of cryptic species (Struck et al. 2018). Thus, species involved in mimicry complexes provide a good testing ground to study the contrasting effects of divergence and stasis in the evolution of cryptic species. However, the systematics of these species often results inaccurate or partial, especially if only based on the external morphology. As a matter of fact, the presence of cryptic species within mimicry complexes has predominantly been reported in some well-studied tropical taxa, which turned out as excellent test groups to evaluate how mimicry affects patterns of diversification and speciation (Pfennig 2012, and references therein;Jiggins 2017). Yet mimicry complexes are not limited to the tropical regions.
In this study, we investigated the presence of cryptic species within a group of mimetic moths occurring in temperate regions whose systematics has long been controversial, the Euro-Anatolian Syntomis. These are distasteful, day-flying moths showing aposematic colourations and belonging to a Müllerian mimicry complex that include also several unrelated co-mimics, e.g. yellow melanic morphs of Zygaena ephialtes, Z. transalpina (Zygaenidae) and Callimorpha dominula (Erebidae) (Bullini et al. 1969;Sbordoni et al. 1979;Zilli 1996). Being mainly based on features of the habitus and genital morphology (Obraztsov 1966), the taxonomy of Syntomis moths appears outdated and confused, evidently flawed by the incorrect delimitation of real species boundaries. Disputes also exist on whether the genus Syntomis should be distinguished from the genus Amata-here we treat Syntomis as a distinct genus, in agreement with part of the recent literature (e.g. Schneider et al. 1999;Igniatyev and Zolotuhin 2005;Anikin et al. 2017;de Freina and de Prins 2018), and with the evidence that Amata is polyphyletic (Przybyłowicz et al. 2019: 629-630) whereas Syntomis include a monophyletic set of Palaearctic and Oriental species characterized by synapomorphies and genetic differentiation (see Schneider et al. 1999). In particular, we focused on S. marjana, a taxon ranging from Provence, Sicily and mainland Italy across the Balkan Peninsula, Ukraine and Southern European Russia (de Freina 2008;Fibiger et al. 2011). Based on their disjunct distribution, and slight differences in size, male genitalia, wing shape and spotting, some populations have occasionally been considered distinct species, i.e. albionica (Provence), quercii (Apennines), kruegeri (Sicily), marjana proper (Balkans) and sheljuzhkoi (Caucasus) (Verity 1914;Turati 1917;Stauder 1928Stauder -1929Obraztsov 1966;Dufay 1970;Igniatyev and Zolotuhin 2005). Preliminary genetic studies showed a substantial divergence between two populations of the taxa kruegeri (Mt Pellegrino, Palermo, Sicily) and marjana s. str. (Stari Grad, Croatia), which were proposed to be distinct species (Cianchi et al. 1980;Bullini et al. 1981). Nevertheless, this genetic evidence was not considered anymore in recent taxonomic reviews, which recognized only a single widely distributed species with more subspecies (de Freina and Witt 1987;de Freina 2008;Fibiger et al. 2011). Such discordance on the systematics of the group claims for comprehensive genetic investigations.
Here we provide the first extensive assessment of the genetic variation across some Euro-Anatolian Syntomis taxa. Firstly, we assessed the genetic variation at both nuclear and mitochondrial markers within broad sense 'S. marjana' and in two closely related, sympatric species, S. phegea and S. ragazzii; then, we used ordination-based clustering analyses and species delimitation methods to circumscribe putative species; finally, we compared the estimated time of divergence among sympatric and allopatric taxa. The aim of this paper is thence to clarify the taxonomy and systematics of Western Palaearctic Syntomis and disclose the extent of cryptic diversity within a mimicry complex occurring in temperate areas.

Sampling and laboratory procedures
The sampling strategy was designed to include most of the geographic variation from Syntomis marjana, S. ragazzii and S. phegea. We analysed 1150 specimens from 34 localities representing most of the described subspecies (we did not analyse S. m. odessana, S. m. sheljuzhkoi and S. p. phegea). Sampling localities and sample sizes are given in Table 1 and Fig. 1. Specimens were preliminarily identified following the characters described in Obraztsov (1966), de Freina and Witt (1987) and Fibiger et al. (2011). Biological tissue was stored at − 80°until subsequent analyses. We also obtained two dry specimens from the related species S. nigricornis and S. aequipunta for use in the phylogenetic analysis (see below).
DNA was extracted from the legs of a subsample of available specimens (Table 1) following the standard cetyltrimethylammonium-bromide (CTAB) protocol (Doyle and Doyle 1990). A fragment from the mitochondrial Cytochrome Oxidase I gene (COI) was amplified and sequenced. The polymerase chain reaction (PCR) primers were REVNANCY and PAT2K837 from Simon et al. (1994). Amplifications were performed in a 10-μL reaction volume containing MgCl 2 (2 mM), four dNTPs (0.2 mM each), two primers (0.2 μM each), the enzyme Taq polymerase (0.5 U, Promega), its reaction buffer (1X, Promega) and 10-100 ng of DNA template. PCRs were conducted with an initial step at 95°C for 5 min, 32 cycles at 94°C for 1 min, 45 s at 57°C (COI) and 1 min at 72°C, followed by a single final step at 72°C for 5 min. Purification and sequencing of PCR products were conducted on both strands by Macrogen Inc. (http://www. macrogen.com), using an ABI PRISM® 3730 sequencing system (Applied Biosystems). All sequences have been deposited in GenBank (accession numbers: MW343496-MW343511).

Allozyme data analysis
The effective number of alleles per locus was calculated with GENALEX (Peakall and Smouse 2006). Allele frequencies, mean observed and expected heterozygosity and proportion of polymorphic loci were computed for each sampling site (population) using BIOSYS-2 (Swofford and Selander 1999). BIOSYS-2 was also used to evaluate departures from the expected Hardy-Weinberg equilibrium for each locus at each sampling site, and the linkage equilibrium between each pair of loci after application of the Bonferroni correction for multiple tests. As no departures from Hardy-Weinberg equilibria were observed, indicating panmixia in all populations, we conducted a population-based assessment of genetic affinities.
A pairwise matrix of unbiased Nei's genetic distances (Nei 1978) among all populations was computed by BIOSYS-2. Clustering of populations was assessed via principal coordinate analysis (PCoA) of Nei's distances among populations, as implemented in GENALEX. PCoAdefined clusters were then compared with candidate taxa, and a matrix of fixed differences (i.e. fully diagnostic loci) among them was generated.
We applied three different species delimitation methods traditionally employed to delimit groups of sequences that potentially correspond to distinct species: (i) Automatic Barcode Gap Discovery (ABGD; Puillandre et al. 2012a) that detects gaps in the distribution of pairwise genetic distances, assuming that it corresponds to a threshold between intra-and inter-specific distances; (ii) General Mixed Yule Coalescent model (GMYC; Pons et al. 2006), which starts from an ultrametric tree and tests whether the branching rates fit better with a coalescent model or a speciation model, using the transition point between speciation and coalescence to delimit species; and (iii) Bayesian implementation of the Poisson Tree Processes (bPTP; Zhang et al. 2013), which also compares speciation and coalescent models but relies on the substitution rates calculated for each node instead of the branching rates. These three methods were chosen because they proved to be effective in recognizing cryptic species either in large and small groups of sequences, and are among the most widely used with mtDNA data (Fujita et al. 2012;Puillandre et al. 2012a, b;Zhang et al. 2013;Schwarzfeld and Sperling 2015).
ABGD analysis was performed on the ABGD web server platform (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb. html), using the default parameters (gap width X = 1.5, prior intraspecific divergences from P = 0.001 to P = 0.1, with 20 steps) and the Kimura-2-parameter (K80) model to compute a pairwise genetic distance matrix. The phylogenetic trees used as input for the GMYC and bPTP analyses were generated using the Bayesian method implemented in BEAST 1.8.1 (Bayesian evolutionary analysis by sampling trees; Drummond et al. 2012). We selected the Yule pure-birth speciation model as tree prior, a strict clock model and the best-fit model of molecular evolution estimated by Jmodeltest 2.1.3 (Darriba et al. 2012) under the Bayesian Information Criterion, i.e. the TIM2 transitional model with a proportion of invariant sites (+I). The Markov chain Monte Carlo (MCMC) length was 10 million generations, sampling trees every 1000 generations. The independence of the effective sample size (ESS values »200) for the estimated parameters was evaluated using Tracer 1.6, after removing the first 10% of samples as burn-in. The consensus tree was then generated by TreeAnnotator 1.8.1 (BEAST package), using the maximum clade credibility criterion, after removing the first 1000 sampled trees as burn-in. We included in the analysis also the sequences of S. nigricornis and S. aequipuncta. The GMYC analysis was then run using the splits R package (Ezard et al. 2015), applying both the single threshold and the multiple threshold methods. bPTP analysis was conducted on the web server platform (https://species.h-its.org/ptp/), with a 100,000 MCMC length and a 10% burn-in.
The phylogenetic relationships among the candidate species were inferred using maximum likelihood (ML) and Bayesian methods. The ML-tree was generated using the algorithm implemented in IQTREE (Nguyen et al. 2014) applying the TIM2+I model of substitution; the robustness of the inferred tree topology was assessed using the non-parametric bootstrap method with 1000 pseudo-replicates, and the SH-like approximate likelihood ratio test (SH-aLRT), also with 1000 bootstrap replicates. The Bayesian phylogenetic analysis was performed using the timecalibrated procedure implemented in BEAST 1.8.1 (Drummond and Rambaut 2007), to get a comparative framework of the times of divergence among the candidate species and to get an approximate historical contextualization of the speciation events. In the absence of internal calibration points derived from fossil data or geologic events directly linked to the evolutionary history of Syntomis,totime-calibratethetreeweusedanuncorrelatedrelaxed clock model and a fixed substitution rate of 0.0015 per site per million of years, i.e. the mean COI substitution rate estimated for arthropods by Brower (1994). Although the use of non-fossilbased calibrations is usually discouraged (Schenk 2016), this option is often the only available for most of the soft-bodied invertebrates,anditcontinuestobeextensivelyusedindatingdivergences in insects, including Lepidoptera (e.g. Hinojosa et al. 2018). We chose the substitution rate by Brower (1994) among the several available ones (Papadopoulou et al. 2010) since it falls in the middle of the ranges estimated forthis gene fragment in insects,and it is so far the most used, facilitating comparative inferences. The Yule pure-birth speciation model was chosen as tree prior, and the TIM2+I as substitution model. Two independent runs were performed, each with a Markov chain Monte Carlo (MCMC) length of 10 million generations, with sampling every 1000 generations. The independence of the estimated parameters (ESS values > 200) and the convergence between runs were evaluated using Tracer 1.6, after removing the first 10% of samples as burn-in. The two runswerecombinedusingLogCombiner1.8.1(BEASTpackage), and an annotated maximum clade credibility tree was computed with TreeAnnotator 1.8.1. The pairwise matrix of population genetic distances D (Nei 1978) is given in Supplementary Tab. 2, whereas a heatmap of the D values is shown in Fig. 2a. The inspection of the Table 2 Reduced matrix of the pairwise mean genetic distance (D Nei 1978) between (below the diagonal) and within (on the diagonal) the four genetic clusters identified by allozyme loci, and number of 100% diagnostic loci between them (above the diagonal). Names of taxa as per ; these values agree with those commonly observed among subspecies . A summarized matrix of the mean genetic distances observed within and among these groups is given in Table 2.

Allozyme data analysis
The ordination-based clustering analysis resulting from PCoA also defined four main clusters of populations (Fig.  2b). The first and second PCoA axes explained 40.48% and 34.10% of the total variance, respectively. Populations of S. m. kruegeri + S. m. quercii + S. m. albionica were grouped in a distinct and well-defined cluster with respect to populations of S. m. marjana, which were also combined into a single cluster; the other two clusters coincided with populations of S. phegea and S. ragazzii.

Sequence data
We amplified and sequenced a 781-bp fragment from the final section of the mitochondrial COI gene from 38 Syntomis individuals. No indels, stop codons or nonsense codons were observed. We found 16 different haplotypes defined by 81 (10.3%) variable positions, of which 68 (8.7%) were parsimony informative. The mean haplotype diversity (h) and nucleotide diversity (π) values for this dataset were 0.950 (± 0.021 SD) and 0.032 (± 0.008 SD), respectively. A full list of the haplotype found within each studied population is presented in Table 1.

Discussion
Both mitochondrial and nuclear markers consistently supported the existence of two distinct, highly differentiated, genetic clusters within S. marjana, distributed West and East of the Adriatic Sea, respectively. The extent of the genetic differentiation between these two clusters and the estimated time of their divergence are of the same magnitude of those occurring between other well-recognized species of this complex. Furthermore, all species delimitation methods clearly point to these two genetic clusters as distinct species. These results are consistent with preliminary evidence showing the genetic divergence between the taxa kruegeri and marjana s. str. (Cianchi et al. 1980;Bullini et al. 1981), and suggest considering the populations of the lineage kruegeri + quercii + albionica as a distinct species, under the oldest available name having priority, i.e. S. quercii Verity 1914 (see below). Our data also highlight further genetic sub-structuring at the intraspecific level, suggesting a role for recent Pleistocene events in triggering population-level genetic differentiation. The taxa kruegeri, marjana and quercii were originally described as aberrations or subspecies of other species, the first two of S. phegea (Ragusa 1904;Stauder 1913), the last one of S. mestralii (a species from the Middle East). Subsequently, kruegeri and marjana were considered valid species by Turati (1917), who ranked quercii as a subspecies of marjana. However, the substantial similarity of these taxa in habitus, habitat preferences (dry grasslands), and phenology led most authors to combine them into a unique polytypic species (Obraztsov 1966;Bertaccini et al. 1997;Igniatyev and Zolotuhin 2005;Freina 2008;Fibiger et al. 2011). Despite this, we detected substantial genetic differentiation between S. marjana s. str. and S. quercii (the latter redefined here to comprise also the taxa albionica and kruegeri) that can be attributed to a speciation event occurred at some point in the Early Pleistocene or before. Our estimate of the divergence time between these two lineages is probably biased by the use of a non-fossil-based calibration of the molecular clock (Papadopoulou et al. 2010;Schenk 2016), but it provides an approximate timing of their split. Indeed, the great changes in bioclimatic conditions characterizing the Plio-Pleistocene transition prompted speciation in several Mediterranean taxa (Hewitt 1996(Hewitt , 2011, including Lepidoptera (Schmitt 2007), and likely affected the evolutionary history of European Syntomis. These changes consisted of a generalized increase of drought and a reduction of the forest habitats in Southern Europe (Hewitt 2011), followed by the beginning of glacial cycles, that resulted in alternate expansion and retreat of xeric and forest habitats (Fauquette et al. 1999;Suc and Popescu 2005;Abrantes et al. 2010). Most European Syntomis are thermophilous (with the partial exception of S. phegea), occur in dry habitats, and are distributed in South-Eastern Europe with several relatives in steppe areas from the Middle-East to Central Asia, that likely represents the centre of origin of this group ("phegea-Gruppe" of Obraztsov 1966). The common ancestor of S. marjana and S. quercii might have benefited from the increase of aridity at the end of Pliocene (Malatesta 1985;Fauquette et al. 1999) by spreading across South-Eastern Europe. Then, with the onset of the first glacial cycles, the cooler conditions and the spreading of broad-leaved forests in the Italian and Balkan Peninsulas (Abrantes et al. 2010) could have trapped the ancestor in the south of both peninsulas, where populations would have found relatively hot and dry refugia. This vicariance process has been claimed as the most likely explanation for the existence of several pairs of sibling species with allo-parapatric distribution in the Italian and Balkan peninsulas (Racheli and Zilli 1985). The Isonzo river in NE Italy and Slovenia, which virtually separates the distribution of these species, has been considered a suture zone and biogeographical boundary for several vertebrate and invertebrate species (Taberlet et al. 1998;Hewitt 1999). However, in our case it is hard to define a boundary between these two species because most of the populations west of the Isonzo river have not been found anymore over the last 50 years and could not be genotyped. Further genetic assays on museum specimens are therefore required.
Pleistocene climatic changes also support the origin of intraspecific genetic sub-structuring. Although our data do not allow fine-scale phylogeographic inferences, they reveal substructuring in S. quercii, S. ragazzii and S. marjana s. str., which can be explained by the interaction between climate fluctuations and the physiographic heterogeneity of the Italian and Balkan peninsulas. In particular, we found genetic differentiation between the Sicilian and Apennine populations of S. quercii, between the central and southernmost Apennine populations of S. ragazzii, and between the northern and southern populations of S. marjana (Fig. 3 and Supplementary Tab. 2). Substantial intraspecific differentiation between Sicilian and mainland populations, and between the southernmost and central Apennine populations, is commonly observed in phylogeographic studies on Italian species (e.g. Canestrelli and Nascetti 2008;Canestrelli et al. 2010;Chiocchio et al. 2019;Scalercio et al. 2020). These patterns can be related to the presence of several geographic discontinuities that acted as ecological and physiographic barriers at some point during the Pleistocene. In fact, during Pleistocene, this area has been affected by several glacio-eustatic oscillations of the sea level (Bonfiglio et al. 2002) that induced repeated fragmentations of populations in distinct sub-refugia. The genetic isolation within distinct sub-refugia led populations toward genetic differentiation in more distinct genetic lineages, a process that boosted the evolution of biodiversity hotspots in southern Italy as well as in the other Mediterranean peninsulas (Hewitt 2011). This scenario, commonly named as "refugia-within-refugia" scenario (Gómez and Lunt 2007), fits well with the geographic pattern of genetic structure observed in S. ragazzii and S. quercii. In contrast, no apparent paleogeographic evidence supports the genetic differentiation observed within S. marjana in a relatively restricted area along the NE Adriatic, with no evident orographic discontinuity. However, this area is characterized by geological discontinuities (Cvetkovic et al. 2016) that affect groundwater availability and, as a consequence, vegetation and bioclimatic conditions (DMEER 2017). These bioclimatic discontinuities likely represented ecological barriers that could have trapped S. marjana populations in distinct climatic refugia, at least during the last part of the Pleistocene. A similar pattern has been outlined in some other species (Hewitt 2011;Pabijan et al. 2015). Further investigations will then be required to track the evolutionary history of S. marjana in the Balkan peninsula during the Pleistocene. We found a general low level of intra-population genetic diversity in all of the analysed Syntomis species, a pattern quite uncommon in insect populations (Schmitt and Seitz 2004). Such low values may reflect the effects of recent bottlenecks, perhaps induced by the last glaciation, in combination with rapid postglacial range expansions. This seems to be the case of S. phegea. Although we did not analyse populations from the northern and the eastern part of its range, the observed low level of population genetic diversity and genetic differentiation at both nuclear and mitochondrial markers support for the hypothesis of a (very) recent and quick range expansion by this species. On the other hand, we found a latitudinal gradient of population genetic diversity in S. ragazzii, with the southernmost populations showing the highest levels of diversity. This pattern may be ascribed to repeated cycles of fragmentation within the southern subrefugia (see above) followed by gene admixture, a model used to explain similar patterns of genetic diversity observed in several other taxa with sub-Apennine distribution (Canestrelli et al. 2010). However, as other causes may explain these patterns in Lepidoptera (e.g. Wolbachia infections, Duplouy et al. 2010) we postpone any speculation to further investigations.
In the phylogenetic relationships outlined by mitochondrial DNA, noteworthy is the recovery of S. ragazzii as more closely related to the Anatolian-Caucasian S. nigricornis than S. phegea, despite S. ragazzii and S. phegea being long considered sibling species and can still hybridize (althought limited to F1, Sbordoni et al. 1982). Vicariance between Southern Apennine and Anatolian taxa is common in Lepidoptera and can be ascribed to several biogeographic processes (Racheli and Zilli 1985) other than a direct geological connection between the two regions, the last of which estimated in the late Miocene (Rögl 1999). Interestingly, the estimated divergence between these two species is more recent than that between the other sister pairs. However, more samples of S. nigricornis and other middle-eastern taxa are needed to resolve the biogeographic history of this connection. Finally, the weak support values recovered for the relationship of S. phegea and S. aequipuncta with the lineages of quercii-marjana and ragazzii-nigricornis stress the need to broaden the analysis to more genetic markers and species to fully resolve the evolutionary history of Western Palaearctic Syntomis (see also Przybyłowicz et al. 2019).
Our results highlight that in Euro-Anatolian Syntomis the degree of genetic divergence is not strictly related with that of phenotypic divergence, as expected for species involved into mimicry complexes. These moths happen thus to be at the crossroad between independent evolutionary divergence and constraints limiting their phenotypic diversification (Leimar et al. 2012). Indeed, species belonging to the same Müllerian mimicry ring share the costs of predation, but they have to be perceived very similar by predators. As a consequence, the external appearance of these species is expected to be under strong selection to maintain a common signal. A different set of characters (e.g. morphological, genetic, physiological and behavioural) will reveal their interplay with the various evolutionary forces, which is in turn affected by the allopatric/ sympatric occurrence of populations of these moths with other co-mimics. Many Syntomis species are sympatric and coexist in the same on very close biotopes, at least to the scale of bird predators' home ranges. S. quercii and S. marjana have alloparapatric distribution, but they both overlap with S. phegea (except in Sicily, where the latter is absent, all records hitherto being misidentifications), which is hugely abundant and has already been shown to be a model also for some distantly related co-mimics (Sbordoni et al. 1979). In this context, the departure from a successful aposematic signal would be strongly counter-selected during and after speciation. Strong ecological constraints would then be responsible for morphological stasis, in the face of substantial genetic divergence, thus eventually leading to the evolution of so-called cryptic diversity (Bickford et al. 2007;Struck et al. 2018). In Syntomis, such stasis also extends to the configuration of the genitalia, all fairly homogeneous across the group, which cannot be modelled by predation or other external agents. This circumstance contrasts with the striking diversification in genitalia seen in other almost indistinguishable moths (e.g. the sympatric species pairs Grammodes geometrica-G. occulta, Dysgonia algira-D. torrida, Noctua fimbriata-N. tirrenica and Cilix glaucata-C. hispanica, to name just a few), suggesting that Syntomis genitalia could have maintained stable configurations due of phylogenetic inertia, as they perhaps do not play any species isolation role.
Mechanical compatibility of genitalia between different Syntomis species has been proved by the detection of rare hybrids S. phegea x S. ragazzii (Sbordoni et al. 1982) and heterospecific fertile matings S. phegea x S. marjana reared from the wild (Rupinpiccolo, NE Italy, A. Zilli, unpublished). Mate recognition and species integrity in Syntomis have, therefore, to rely on other systems than mechanical compatibility. As seen above, visual stimuli are unfit candidates, as any departure from the shared pattern would weaken the efficacy of their aposematic signalling. Bickford et al. (2007) recognize that cryptic diversity is particularly common in taxa not using visual mating signals. Sheer evidence that Syntomis are unlikely to use visual mate recognition systems is based on field observations of pseudocopulae with distantly related species belonging to an alternative mimicry ring based on different colours (red and black), i.e. some burnet moth species of the genus Zygaena (Zygaenidae). Such mismatings are not uncommon during the peaks of abundance of these moths, when they cluster together on flower heads and get confused by congestion of pheromonal plumes all around and tactile abdominal stimuli that lead them to clasp almost anything within their reach (Lees and Zilli 2019).
The occurrence of temporal displacement as reproductive isolating mechanism among Syntomis has been suggested by Sbordoni et al. (1982), who reported partial temporal displacement between S. phegea and S. ragazzii in Central Italy, and a complete displacement between S. phegea and S. marjana in the Karst Plateau. However, during the fieldwork campaigns run by us in the Karst in different years, S. phegea and S. marjana were often found flying simultaneously, occasionally indulging in heterospecific copulae. Furthermore, S. phegea, S. ragazzii and S. quercii can be found flying together in a same biotope in the Apennines. Therefore, the issue of temporal isolation between sympatric Syntomis is worth of further investigation, as a trade-off situation is expected between selective pressures promoting synchrony-Müllerian mimicry being most effective under such circumstance-and opposite ones enhancing temporal displacement as a way of ethological reproductive isolation.
Preliminary breeding experiments run by us (not shown) are fully compatible with chemically mediated interactions between sexes in Syntomis moths, with females typically resting in a calling posture and attracting males. Bendib and Minet (1998) demonstrated the presence of sex pheromone glands in female S. phegea. In contrast, Schneider et al. (1999) highlighted that males of all Syntomis species bear androconial organs at the base of forelegs. Thus, all evidence points to chemically mediated interactions as the main driver of mate recognition in Syntomis moths.

Taxonomic remarks and conclusions
Our research confirms that in Syntomis, one of the taxonomically most controversial groups of European moths, one additional cryptic species is involved. This is the case of a group of populations living in the Italian peninsula, Sicily and Provence that has been noted by different authors as quercii, kruegeri and albionica, respectively. Genetic information points out that these populations have to be combined in a distinct species under the name Syntomis quercii Verity 1914 (bona sp., stat. nov.). In fact, when first named by Ragusa (1904), the oldest name kruegeri was infrasubspecific and before it was made available for nomenclature by Turati (1917), Verity's (1914) description intervened.
Further genetic sub-structuring has been detected in S. quercii, S. ragazzii and S. marjana. Automatically defining subspecific rank to genetically differentiated populations is considered a common mistake (see Zilli 1996, and references therein). However, it is worth noting that at least in S. ragazzii, diversification in southern Italy (Calabrian Apennine) had already been pointed out with the description of two subspecies (silaensis Obraztsov 1966 and asperomontana [Stauder & Turati], 1917). Furthermore, partially distinct genetic makeups also characterize the Sicilian populations of S. quercii, which coincide with the original concept of nominotypical kruegeri.
Research on the Euro-Anatolian taxa of Syntomis highlights how the Palaearctic species of this group provide an interesting example to study trade-offs between morphological stasis and genetic divergence in complexes involved into mimicry relationships of Müllerian type. New studies have to be focused on understanding the genomic basis of speciation and origin of cryptic diversity within this group of moths.
Funding Open access funding provided by Lund University. AC was supported by the SYNTHESYS Projects GB-TAF-2574, awarded by the European Community Research Infrastructure Action under the EU FP7 Framework.
Data availability The dataset generated during the current study is available in the GenBank repository (see corresponding GenBank numbers in the "Material and Methods" section)

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.

Code availability Not applicable
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.