DNA Extraction and Sequencing

The bodies of the specimens used in the DNA analyses were stored in 99% ethanol at –20°C and wings were kept separately as vouchers. The DNA analyses were conducted employing sequences retrieved from 191 individuals, including all the species of the M. phoebe group. Sampling sites in the Iberian Peninsula are plotted in Fig. 1; here, the specimens were identified using the wg phylogeny ( Supp Fig. S1 [online only]). In total, we used 160 COI and 180 wg, 111 RPS5, 87 MDH, and 90 EF-1α ( Supp Table S1 [online only]).

Total genomic DNA was extracted using Chelex 100 resin, 100–200 mesh, sodium form (Biorad), under the following protocol: one leg was removed and introduced into 100 µl of Chelex 10% and 5 µl of Proteinase K (20 mg/ml) were added. The samples were incubated overnight at 55°C in a shaker and were subsequently incubated at 100°C for 15 minutes.

Primers and PCR protocols used for the amplification of COI (barcode region), wg, EF-1a (three fragments), RPS5, and MDH are written in  Supp Tables S2 and S3 (online only). Universal tails were included in all primers. PCR products were purified and Sanger sequenced by Macrogen Inc. Europe (Amsterdam, the Netherlands). All sequences have been deposited in GenBank ( Supp Table S1 [online only]).

Phylogenetic Reconstruction

Sequences were visualized, edited, and aligned in Geneious Prime 2019.0.3 ( https://www.geneious.com). A haplotype network of the COI barcode region was created in POPART v1.7 (Leigh and Bryant 2015) under the TCS method. The best-fitting substitution model was estimated in jModelTest (Darriba et al. 2012) under the Akaike information criterion.

The alignments of the nuclear genes were concatenated and a phylogeny was reconstructed using BEAST v2.5.0 (Bouckaert et al. 2014). Distinct partitions were selected for each gene using the best model indicated by jModelTest and four rate categories if included gamma and base frequencies were estimated. Parameters were estimated using two independent runs of 30 million generations each and convergence was checked with TRACER 1.7.1 (Rambaut 2018). A burn-in of 10% was applied.

A maximum likelihood (ML) inference was obtained with the concatenated (partitioned) alignment of nuclear loci using RAxML v8.2.11 (Stamatakis 2014). RAxML was also used to retrieve a phylogeny for every individual nuclear gene, which was employed in ASTRAL. A thorough bootstrapping was employed and we selected a GTRGAMMA model and 1,000 bootstrap replicates.

Species Delimitation

Joint Bayesian species delimitation was conducted using the program BPP (Yang 2015). The method uses the multispecies coalescent model to compare different models of species delimitation (Yang and Rannala 2010, Rannala and Yang 2013) in a Bayesian framework, accounting for incomplete lineage sorting due to ancestral polymorphism and gene tree-species tree discordance. Considering that, in butterflies, the median of the theta is 0.0160 (ranging from 0.004 to 0.043) mutations/site (Mackintosh et al. 2019), a sensible diffuse theta prior would be IG(3, 0.045).

Fig. 1.

Iberian sampling sites of Melitaea pseudornata stat. nov. and Melitaea phoebe. The relationship between flight time and altitude is depicted. All identifications are based, at least, on the sequencing of the wg gene.

Assuming neutrality, the mutation rate for butterflies has been estimated to be about 2.9 × 10^–9 mutations/site/year (Keightley et al. 2015). A divergence time between 3.3–9.6 Mya (95% credibility interval) has been estimated for the M. phoebe species group (Tóth et al. 2017). By multiplying these values, we considered a divergence between the root of the species tree and the present time (tau) between 0.0096 and 0.0278 mutations/site and we assigned a diffuse tau prior of IG (3, 0.03).

The other divergence time parameters are specified by the uniform Dirichlet distribution (Yang and Rannala 2010: equation 2). We selected 500,000 MCMC and 50,000 burn-in and the analysis was run twice to confirm consistency between runs.

Genitalia Morphometry

Geometric morphometry on the processus posterior (male genitalia) was applied to determine the morphological relationships between M. pseudornata and the other species of the M. phoebe species group. In total, 315 specimens have been measured. These individuals are partly identical with the material used in Tóth and Varga (2011), but it has been completed with 12 specimens from Iberia that were identified based on nuclear DNA data ( Supp Table S1 [online only]).

A standard genital preparation method was followed. The abdomens were removed and heated in 15% KOH solution in 80°C for 30 min. Next, genitalia were cleaned and dehydrated in ethanol and mounted in euparal. Genitalia slides were digitalized using a stereo-microscope and a digital camera.

TpsDig2 was used to record nine fixed landmarks at the tips and origin of the main processus (Tóth and Varga 2011). The raw coordinates were transformed using Procrustes generalized least squares using geomorph R package (Adams et al. 2021).

All further analysis was performed based on the transformed coordinates. We used linear discriminant analysis (LDA) to determine the morphological relationships between the studied taxa using MASS R package (Venables and Ripley 2002). Leave-one-out cross-validation classification was also used to quantify the classification success. The significance of the visible pattern was analyzed by pairwise permutational MANOVA using Bonferroni-corrected significance levels using RVAideMemoire R package (Hervé, 2021). Average landmark coordinates of the processus posterior were calculated for all the studied species, then a PCA was performed on these mean shapes.

Flight Time Data (Navarre, Northern Spain)

Adult records were retrieved from the data of two transects located at Taxoare (Aranguren, Navarre). In this locality, only orange-headed caterpillars were found, and DNA results classified all the sequenced specimens as M. pseudornata. Adults were counted following a standardized methodology that consists of sampling every two weeks, from the first week of April to the last week of September, along a walked transect at distances of 2.5 m on both the sides and 5 m ahead of the recorder (Pollard and Yates 1993).

Nomenclature

This paper has been registered in Zoobank ( www.zoobank.org), the official register of the International Commission on Zoological Nomenclature. The LSID (Life Science Identifier) number of the publication is: urn:lsid:zoobank.org:pub:D1410808-7450-4190-BA1A-C57BA477AD46

Phylogenetic Inference and Species Delimitation

The haplotype network based on the barcode region of the COI (Fig. 2A,B) showed that M. ornata and M. pseudornata shared barcode haplogroups with M. phoebe individuals. However, M. ornata from the Balkan Peninsula and eastwards maintained very distinct barcodes. Some individuals of M. pseudornata displayed exclusive haplotypes within Haplogroup 2, although separated from those of M. phoebe by only one or two substitutions.

Nuclear markers retrieved all the taxa as monophyletic with a posterior probability (PP) of 1 ( Figs 2B, Supp Fig. S2 [online only]) and a bootstrap value (BS) > 70, except for M. phoebe in the ML inference, which had a lower support ( Supp Fig. S3 [online only]). Melitaea pseudornata was recovered sister to M. ornata in the Bayesian and ML inferences. Interestingly, the Bayesian phylogeny retrieved two clades within M. phoebe: a clade formed by all the Iberian specimens and a French individual (PP = 0.82), and a clade with the rest of individuals (PP = 0.77).

Both BPP runs selected the seven-species hypothesis as the most likely (PP = 1). All the taxa of the M. phoebe species group and the Iberian taxon were separated as distinct species.

Fig. 2.

A. Haplotype network of the Melitaea phoebe species group based on the barcode fragment of the COI gene. The area of the circles is proportional to the number of sequences they represent. B. European distribution of the three haplogroups of Melitaea pseudornata stat. nov., M. phoebe and Melitaea ornata. C. Bayesian inference phylogeny obtained using a partitioned alignment of the nuclear genes EF-1a, MDH, wg and RPS5. Posterior probabilities at species-level and for higher relationships are indicated, and scale units are presented in substitutions per site.

Genitalia

Despite interindividual variability was present in this group of taxa, the processus posterior of M. pseudornata showed on average a unique shape, distinct to M. ornata and M. phoebe (Fig. 3B–D). The average shape of the processus posterior was very similar to M. ornata as the inner process have the similar length. On the other hand, it showed some similarities to M. abbas and M. punica as the outer process pointed upward. On the LDA scatterplots, the centroid of M. pseudornata was positioned at a distance from all the other groups but close to M. ornata and M. phoebe considering the first three axes, which explain 74.81 % of the variance between groups (Fig. 3). The first axis explained 38.17%, the second 20.53%, and the third 16.11% of variance between groups. Leave-one-out cross-validation classification assigned 77.78% of the individuals correctly in the case of M. pseudornata ( Supp Table S4 [online only]). This success was very similar to the overall precision of the classification in the group, which was 78.34%.

The pairwise permutational MANOVA indicated statistically significant differences (P < 0.05 in all cases) in the processus posterior shapes among the taxa analyzed ( Supp Table S5 [online only]). M. pseudornata showed a similar morphological differentiation compared to the closely related species ( Supp Table S6 [online only]).

Fig. 3.

A. Landmarks on the processus posterior of a Melitaea phoebe male genitalia. B. Typical processus posterior shapes of Melitaea pseudornata stat. nov., M. phoebe and Melitaea ornata. C. PCA built using the mean shapes of all the studied species. Species mean shapes are shown. D. Linear discriminant analysis (LDA) scatterplot for the studied taxa. Group centroids are shown. All the groups had equal prior probability independently from the sample size of the group.

The Iberian Taxon, a New Species

The results here presented suggest that the Iberian taxon should be elevated to the species status. First, molecular evidence (nuclear markers) retrieved the Iberian individuals as a monophyletic clade, well-diverged from M. ornata, although sister to it (Figs. 2C,  Supp Figs. S2–S3 [online only]); species delimitation analyses supported the specific status for this clade. Second, differences in the genitalia between the Iberian taxon and M. ornata were comparable to those found interspecifically among other species of the group (Fig. 3B–D;  Supp Tables S4–S6 [online only]). Third, their phenology is distinct since M. ornata has apparently only one generation each year—second generations have been obtained only in captivity (Russell and Pateman 2013, Russell et al. 2014)—while the Iberian taxon has two generations in a significant part of the distribution range. Worth mentioning, there is no evidence of a close relationship of the Iberian taxon with the north African M. punica, a hypothesis that was suggested due to similarities present in the caterpillars (Sánchez Mesa and Muñoz Sariot 2017, Muñoz Sariot and Sánchez Mesa 2019a).

The denomination of this novel species, however, is not straightforward. Russell et al. (2020) attributed several taxa to the Iberian M. ornata-like taxon. They based their proposals on the external morphology of the adults, but in these traits the Iberian taxon cannot be reliably differentiated from M. phoebe. Several of our specimens showed typical ornata-like characteristics but, based on nuDNA, they proved to be M. phoebe, and vice versa. The thickness of the tip of the antennae seems not to be a defining trait either because thick tips are found in both taxa (see Sánchez Mesa and Muñoz Sariot 2017, Muñoz Sariot and Sánchez Mesa 2019a), although thicker shapes are more common in M. pseudornata than in the Iberian M. phoebe. Overall, the evidence available points that this taxon can only be reliably distinguished by nuDNA data and by the reddish head of the last (L7) instar caterpillars (Sánchez Mesa and Muñoz Sariot 2017)—distinct to M. phoebe, with invariably black head, and to M. ornata, with reddish head from L4 to the last instar (Russell and Tennent 2016). Given the absence of these data in the taxonomic proposals made by Russell et al. (2020), we think that further analyses of the type specimens are required in order to confirm the identifications. In consequence, here we used the name of the first taxon whose identification was based on the color of the head of the caterpillars, which is pseudornata (Muñoz Sariot and Sánchez Mesa 2019a, b). Thus, we tentatively name the novel species as M. pseudornata Muñoz Sariot & Sánchez Mesa, 2019, stat. nov.

The Distribution Range of M. pseudornata

Identifications based on the wg gene ( Supp Fig. S1 [online only]) confirmed the presence of M. pseudornata across most of the Iberian Peninsula, apparently restricted to areas of oceanic influence and/or mountain ranges. So far, it has been only found in Spain. This distribution is similar to the prediction made by Tóth et al. (2013, 2017), who showed through ecological niche modelling analyses that virtually all the Iberian Peninsula represents a climatically suitable habitat for the sibling species M. ornata. Although further exploration is required, M. pseudornata was the sole species found in a vast area of north-western Spain. In contrast, only M. phoebe was found in the south-west of the Iberian Peninsula (although sampling in this region is low) and in Catalonia (except in the southern mountains of Els Ports, where it is replaced by M. pseudornata). Both species are locally sympatric in some mountain ranges in the Baetic System (south-eastern Iberia), but this seems not to be usual and the pattern documented agrees with a situation of parapatry, in which contact zones may reflect some kind of competition or incompatibility.

The ranges of the species pair M. pseudornata and M. ornata match with a distribution pattern typically produced by glacial cycles, even if initial divergence predates them (Ebdon et al. 2021). Glacial periods caused the isolation of populations in the southern peninsulas, which promoted allopatric differentiation and, sometimes, speciation (Hewitt, 2000). In the Iberian Peninsula, there are about twenty butterfly species that have a sibling widespread through Europe (Dincă et al. 2015). Many of them establish contact zones around the Pyrenees—typically in the Ebro River valley, the Pyrenees themselves, or in S. France—such as the pairs Iphiclides feisthamelii (Duponchel, 1832) (Papilionidae)–Iphiclides podalirius (Linnaeus, 1758) (Gaunet et al. 2019) or Aricia cramera (Eschscholtz, 1821) (Lycaenidae)–Aricia agestis ([Denis & Schiffermüller], 1775) (Vodă et al. 2015). In our case, the existence of a contact zone cannot be determined since it is unknown whether M. pseudornata is present in France; meanwhile, the closest area where M. ornata has been reported is Provence. Hence, as far as we know, M. pseudornata is allopatric with respect to M. ornata.