Analysis of the Holarctic Dictyoptera aurora Complex (Coleoptera, Lycidae) Reveals Hidden Diversity and Geographic Structure in Müllerian Mimicry Ring

Simple Summary We evaluated the red net-winged beetle populations’ morphological and genetic divergence within the Holarctic region. In contrast with relatives, D. aurora occurs in an exceptionally large range and very different ecosystems. In Northern America, we found an earlier undetected cryptic species isolated by the Bering Strait since the mid-Miocene. D. aurora colonized Fennoscandia from at least two refugia after the last glacial maximum. The absence of morphological differentiation is supposedly affected by the selection for similarity in the Müllerian mimicry ring. The results exemplify the phylogeographic history in the Holarctic region and contribute to understanding the morphological stasis in an extensive circumpolar range. Abstract The elateroid family Lycidae is known for limited dispersal propensity and high species-level endemism. The red net-winged beetle, Dictyoptera aurora (Herbst, 1874), differs from all relatives by the range comprising almost the entire Holarctic region. Based on a five-marker phylogeny and 67 barcode entries (cox1-5′ mtDNA) from the whole range, we recovered two genetically distinct species within traditionally defined D. aurora and resurrected the name D. coccinata (Say, 1835) as the oldest available synonym for Nearctic populations. Yet, no reliable morphological trait distinguishes these species except for minute differences in the male genitalia. D. coccinata is a monophylum resulting from a single Miocene dispersal event, ~15.8 million years ago, and genetic divergence implies long-term isolation by the Bering Strait. Far East Asian and west European populations are also genetically distinct, although to a lower extent. Two independent colonization events established the Fennoscandian populations after the last glacial maximum. Besides intrinsic factors, the high morphological similarity might result from stabilizing selection for shared aposematic signals. The rapidly accumulating barcode data provide valuable information on the evolutionary history and the origins of regional faunas.


Introduction
Molecular data have become essential sources of information in biodiversity studies, and they are especially valuable if diagnostic morphological traits are contentious [1][2][3][4][5]. Then, DNA data can solve the delimitation of cryptic species [2,6,7]. Mitochondrial markers are also widely used for the rapid genetic screening of biodiversity, and the method can be increasingly used with next-generation sequencing [8][9][10][11]. Although the mtDNA resolution to detect subtle intraspecific differentiation is low [7,[12][13][14], even short fragments are sufficient for obtaining the first clue to shallow relationships and identification of discrepancies between traditional morphology-based taxonomy and genetic diversity [2,15,16].
Here, we focus on an exceptionally widespread red net-winged beetle, Dictyoptera aurora (Herbst, 1784) (Erotinae: Dictyopterini) that belongs to the predominantly tropical gene was performed with ModelFinder [45,46] implemented in IQ-TREE2 (-MFP option). The GTR model was applied in all analyses. Ultrafast bootstrap [47] values were calculated in IQ-TREE (options -bb 5000) to assess nodal supports for focal relationships.
Further, we generated a haplotype network based on the cox1-5 fragment using the TCS algorithm [48] in PopART [49]. The output was processed in Adobe Illustrator.
The dated tree of the genus Dictyoptera Latreille was taken from the whole family analysis (the dataset used by [50] and expanded by unpublished data). Beast v.1.8 [51] was used to date splits among the Lycidae family with calibrated separation of Lycidae + Iberobaenia to 176 my. The following settings were employed: the unlinked partitions (genes and codon positions were partitioned in the (1 + 2),3 scheme, the HKY + I+G model, and the birth-death speciation prior. The run was computed for 5 × 10 7 generations with a sampling frequency of 10,000. The ESS values and plateau phases were checked using TRACER v.1.681 [52]). The maximum credibility tree was generated with TREEANNOTA-TOR v.1.8.1 [51].

Morphological Study
About 100 individuals deposited in the senior author's laboratory from the Nearctic and Palearctic regions were investigated to identify putative diagnostic characters and intraspecific variability. We specifically compared the shapes of the antennae, pronotum, and elytral costae. Further, we dissected multiple individuals from the whole range to study the intraspecific variability in the male genitalia. The voucher specimens were relaxed in 50% ethyl alcohol for an hour and then detached abdomens were treated with a hot 10% KOH solution to remove muscles and body fat. The external characters and genital morphology were observed under an Olympus SZX-16 microscope and photographed by an attached digital camera. We used Helicon Focus (www.heliconsoft.com, accessed on 21 January 2022) and Photoshop 13.0 (Adobe Inc., San Jose, Ca, USA) to assemble photo stacks in figures showing the principal morphological traits.
The individuals used for morphological investigation and sequenced specimens bearing the UPOL voucher number are deposited in the collection of the Laboratory of Biodiversity and Molecular Evolution, CATRIN-CRH, Olomouc.

Phylogeny, Genetic Diversity and Distribution
In our five-gene phylogenetic analysis, we find two distant clades represented by the individuals traditionally assigned to the morphospecies D. aurora ( Figure 1). The North American populations are represented by a single individual from the Eastern USA, and the Euro-Asian populations by three individuals from the Russian Far East and four from Central Europe ( Figure 1C). The dating analysis using the cox1-3 , rrnL, nad5 mtDNA dataset suggests that Nearctic and European populations separated already in the mid-Miocene, 15.75 million years ago ( Figure 1D). Their generic divergence is documented by uncorrected pairwise distances in mitochondrial and nuclear markers (complete SSU RNA 0.3%, D2 loop of LSU rRNA 0.3%, rrnL mtDNA 4.3%, cox1-3 13.8%, nad5 9.6%). We suggest the presence of D. aurora in the Western Palearctic shortly before the end of the Miocene, 6.23 mya when the Far East and European populations split ( Figure 1D).
Further, the cox1-5 mtDNA (barcode) fragments were analyzed from the whole range. The uncorrected pairwise distances of the barcode cox1-5 fragments reached 12.5% (Table S1). The phylogenetic network based on the cox1-5 fragment confirms separate clusters of haplotypes corresponding to the main branches recovered by the phylogenetic analysis ( Figure 1B,C,E). Further, we identified European populations' relatively high intraspecific genetic diversities. The genetic divergence between Fennoscandian populations reaches 7.2% in cox1-5 (barcode) fragment. Additionally, we constructed the barcode networks for three lycid species occurring in Central Europe and Fennoscandia ( Figure 1F-H). L. sanguineus and Pl. minuta comprise genetically close populations ( Figure 1F,G). P. nigroruber contains genetically divergent populations in Fennoscandia ( Figure 1H).

Differential Diagnosis
Both D. coccinata and D. aurora are morphologically highly variable species, and we have not found any external traits for identifying these sister species (Figures 2 and 3). We only noted minor, potentially unreliable, differences in the shape of the phallus. The basal process (pointing ventrally to the phallobase) is slenderer in D. coccinata. Still, there are intraspecific differences, and although the trait is very apparent in some individuals ( Figure 3M), it is hardly observable in others ( Figure 3O). Additionally, the ventral edge of the D. coccinata phallus is concave to straight in the basal half, unlike the simply concave edge in the whole length of D. aurora (observed in the ventrolateral view; Figure 3G-Q).

Intraspecific Variability
High intraspecific variability of external characters includes the shape of the male antennae, pronotum, and the arrangement of elytral costae (Figure 2A-J). The males have

Intraspecific Variability
High intraspecific variability of external characters includes the shape of the male antennae, pronotum, and the arrangement of elytral costae (Figure 2A-J). The males have slenderer antennae than females, but we noted differences in the relative length of antennomeres also in males ( Figure 3A-E). Similarly, we observed differences in the shape of the pronotum ( Figure 2E-J). Unlike fully sclerotized beetles, the asymmetry of the pronotum is commonly encountered, and defective structures are common in soft-bodied lycids ( Figure 2F). The net-winged beetles are so weakly sclerotized that dry-mounted individuals have regularly deformed elytra. Then, it is difficult to compare the strength of elytral costae. Therefore, we photographed the specimens immediately after mounting to show the differences (Figure 2A-D). The specimens in Figure 2A,C,D have slightly stronger primary costa 2 that reaches the apex of elytra; the specimen in Figure 2B has merged primary costae 2 and 3. The primary costa 2 is not only shorter but also relatively weaker. Similarly, the male genitalia is variable. There are differences in the basal part of parameres (pointed in Figure 3L,M and straight in Figure 3P,Q), the conspicuousness of thorns at the inner margin of parameres (acute in Figure 3L,M,G, obtuse in Figure

Biology
D. coccinata is a locally common forest species in mid to high-latitude mountains; larva is unknown. Adults are active from March in the south until late June in the hi mountains and the northern part of its range.

Biology
D. coccinata is a locally common forest species in mid to high-latitude mountains; its larva is unknown. Adults are active from March in the south until late June in the high mountains and the northern part of its range.

Distribution
Nearctic Region: USA (from northern Florida and Eastern Texas to the Canadian border; Arizona, New Mexico, eastern slopes of the Rockies, southwestern South Dakota, Alaska) and Canada to the northern limit of spruce forests. Most records are available from the Rockies and Appalachian Mountains [19]. Further information was taken from www.inaturalist.org, accessed on 2 July 2022. The records close to range limits were revised, and misidentifications were excluded). The northernmost record is known from central Alaska (Chena River, 65.9 • N).

Note
The mature larva of D. aurora was reported by Kazantsev and Nikitsky [31]. Here, we describe the second instar for comparison.

Biology
The larvae live in rotten, red-colored wood that remains wet for most of the year; typically in large stumps and trunks. Larvae were found in tree roots up to 20 cm deep in the soil. The species is locally common in mountain forests but occurs in lowlands in high-latitude regions. Larvae pupate in autumn, and adults overwinter under the bark or in wood crevices, usually in aggregations, which sometimes contain up to several dozens of individuals. They leave the wood in the early spring. Adults remain primarily inactive, either sitting on rotten wood or herb leaves under the forest canopy, close to the place where they developed. The flight of adults is slow, restricted to shaded situations and favorable weather conditions (high humidity, no wind). The adults occur in nature for about three weeks, from late April until June, depending on the local climate. The males die shortly after copulation, and the females live for another week. The females lay eggs in rotten wood, usually in shallow holes. The first instar larvae hatch after two weeks and immediately actively search for shelter deep in rotten wood. The larvae do not build tunnels and use crevices to move within the substrate. Larvae were growing rapidly in moist wood, but we did not manage to observe the whole life cycle, and the length of development remains unknown.

Taxonomic Decision
There are seven available species-group names in the D. aurora complex [56][57][58]60,61,63,64]. Kleine [37] broadened the concept of D. aurora and considered other names as junior subjective synonyms. Phylogenetic analyses show that this complex consists of at least two genetically distinct clades geographically separated by the Bering Strait. We consider the Nearctic clade a separate species (see Results: phylogeny). Several available synonyms are listed for D. aurora, and the oldest available name referring to Nearctic populations is Omalisus coccinatus Say, 1835. Therefore, we designate the Nearctic species as Dictyoptera coccinata (Say, 1835). Although we found genetic differentiation between Russian Far East and European populations, we refrain from the resurrection of D. superba (Motschulsky) due to the lack of high-resolution nuclear markers and information on the eventual sympatric occurrence. Similarly, the cox1-5 genetic difference between sympatrically-occurring individuals of D. aurora from Fennoscandia (Table S1) surpasses the widely accepted 2-3% threshold for the species delimitation [1], but see [66] for the genetic differentiation of geographically distant populations). The taxonomic status of European populations can be reconsidered when D. aurora is densely sampled in the whole continent.
The holotype of O. coccinatus is unavailable as Say's collection was destroyed [67]. The name has been listed as a junior subjective synonym of D. aurora since Kleine [37], but due to an error, it is listed as a synonym of Lygistopterus sanguineus in the GBIF database (https: //www.gbif.org/species/161597475; accessed on 1 June 2022). L. sanguineus does not occur in the Nearctic region, and the original description clearly states that the pronotum bears the median areola [56]. We refrain from a designation of the neotype as the present study is not a taxonomic revision, and we consider the identity of the species uncontentious [68].

Morphological and Genetic Differentiation
Morphology is the traditional, cheap, and rapid source of diagnostic traits, but morphology-based identification has limitations and warrants validation by independent molecular data [2,15,65]. Our study of the very widespread morphospecies D. aurora reveals yet unknown genetic diversity (Figure 1) and the absence of undisputable diagnostic traits that would discriminate the genetically distant populations as candidate species (Figures 2 and 3).
The intraspecific morphological variability surpasses identified differences between both species (Figures 2 and 3). We found only very subtle differences between Nearctic and Palearctic specimens in the shape of the basal part and the ventral edge of the phallus ( Figure 3F-Q). The students of net-winged beetles have never mentioned these. Indeed, they would not be considered sufficient to separate the Nearctic and Palearctic populations as distinct biological species if we do not have molecular data indicating genetic differentiation, reciprocal monophyly, and long-term geographical isolation. In contrast with these, male antennae ( Figure 3A-D), the pronotum ( Figure 2E-J), the relative strength and length of longitudinal elytral costae (Figure 2A-D), and the general shape of the male genitalia ( Figure 3F-Q) are quite variable within populations. D. aurora's genetic divergence supports a hypothesized causality between low vagility and high diversification rate [7,12,69]. Yet, the placement of all populations in a single morphospecies challenges the intuitively expected link between genetic and morphological diversification.
The levels of mtDNA divergence detected in the D. aurora clade are comparable to those often found between morphologically species pairs in other insects [1,15,66] and also between the sister species of net-winged beetles [12,27]. The genetic diversification between Nearctic and Palearctic populations of the D. aurora complex is higher than between many species' pairs identified in previous studies. The split between the Nearctic and Palearctic populations was estimated to be 15.75 mya ( Figure 2B). Therefore, we split the D. aurora complex into two (semi)cryptic species with~16 million years of independent evolution.
The genetic divergence exceeds thresholds used for algorithmic identification based on barcode data. Eleven bins (clusters) are proposed in the BOLD database for the D. aurora complex based on 67 barcodes (www.boldsystem.org, accessed on 1 June 2022). It is an anomaly as most cox1-3 mtDNA clusters approximately correspond with biological species, or the species are split into only a few genetically distinct groups if data originate from disjunctive areas [2,66].
Although we deal with a single morphospecies and morphological stasis is a macroevolutionary and paleontological concept [70,71], it is worth discussing the processes possibly leading to such an undisputable lack of morphological divergence (Figures 2 and 3). All net-winged beetles share very similar biology, yet most closely related species are morphologically distinct [27,32,72]. For example, their distinctiveness is often a result of advergence to various Müllerian co-mimics [28,50,[73][74][75]. Further diagnostic traits are regularly found in the male genitalia that plays a role in sympatric speciation [76,77]. Therefore, intrinsic constraints undoubtedly play a role in morphological stasis, but ecological factors also affect the divergence of the phenotype.
The origin of the D. aurora complex and its striking morphological uniformity stands in contrast with a new Ferreira's hypothesis that stable paleoenvironment is a driver for conserved morphology in a paedomorphic lineage [44]. The hypothesis is based on a single, 15 million years (my) old, Dominican amber neotenic species, its extant relative placed in the same genus, and the stable position of Hispaniola in the tropics. All populations of the D. aurora complex are morphologically uniform (Figures 2 and 3), yet the morphological stasis due to stable environmental conditions is exceptionally improbable as D. aurora occurs from Northern African cedar forests to tundra low shrub ecosystems beyond the polar circle in Fennoscandia and D. coccinata from dry shrub and forest vegetation in Arizona to the polar circle area in Alaska and the humid temperate to subtropical forests of Georgia and northern Florida. Following the argumentation presented by Ferreira et al. [44], we would have to analogically propose that the unstable and diverse environments select for the morphological stasis of the D. aurora complex.
The linkage of morphological stasis and neoteny is another aspect of Ferreira's theory that is contradicted by the conserved morphology of fully metamorphosed D. aurora (Figures 2 and 3). First of all, the neotenic net-winged beetles are morphologically disparate. They represent 2% of net-winged beetle species diversity but account for 25% of the described morphology-based supergeneric taxa (12 of 48 proposed; [18,78,79]. We might also consider that the ontogenetic modifications have a substantial impact on the phenotypes, and neoteny, or more broadly pedogenesis, are important macroevolutionary drivers producing morphological disparity [80]. The net diversification rate of neotenics must be considered, i.e., speciation and extinction [81]. The high extinction rate is possibly the reason why a low net diversification was observed in neotenics after the colonization of new landmasses and an initial short-term diversification burst [82,83]. Another question is the time factor. The morphologically uniform D. aurora complex split from its closest, phenotypically similar relatives at least 30 my before the time of Dominican amber (54 versus 23 my or 15 my) [44,84]. An ancient origin of the Dictyoptera-like phenotype is also supported by Priabonian Erotinae [85][86][87]. The earliest split between the cryptic species, D. aurora and D. coccinata, is estimated at 15.75 mya (7.18-26.06 mya 95% confidence interval; Figure 1D). The Dominican leptolycine species underlying Ferreira's hypothesis is slightly younger than the earliest split within the morphospecies earlier designated as D. aurora. The age of the Dominican fossil species disqualifies it from a long-term stasis evaluation. The dated phylogenies of net-winged beetles have shown that many genera with morphologically uniform species can be traced to the Paleogene [13,28,50]. As the studies were published only recently, Ferreira and collaborators [44] possibly missed them. Still, the authors did not compare the time since the establishment of the conserved morphology of their Dominican net-winged beetle with other beetles known for morphological stasis. Examples include morphological similarity of extant and early Mesozoic elateroids; mid-Cretaceous and Priabonian lycids [88][89][90][91][92][93][94] but not [95] who reported a Cretaceous tenebrionoid beetle as a lycid by error; extant and mid-Cretaceous lymexylids [96,97], jacobsonids [98], or micromalthids [99]. Growing information on Kachin amber provides further examples in various lineages [100].
The ancient morphological stasis was also recovered by dated phylogenies [101,102]. To sum up, there is plenty of evidence that morphological similarity of closely related species (within genera, usually with the earliest common ancestor hypothesized in the Paleogene) is a rule and not evidence for a causal link between conserved morphology, the neotenic development, and stable habitat. The D. aurora complex is one of the welldocumented examples of the long-term persistence of an ancient phenotype under very diverse and unstable environmental conditions.

Genetic Diversity in European Populations
The Palearctic populations of D. aurora do not represent a genetically homogenous group, and they split into the western and eastern clades in the five gene analyses ( Figure 1C). Additionally, two distinct cox1-5 haplotype groups are identified in Fennoscandia ( Figure 1B). As net-winged beetles depend heavily on rotten wood, we assume that the range of D. aurora was pushed to the south along with forest ecosystems during glacial maxima [103]. The high genetic diversity in the Fennoscandian population shows that the area only relatively recently reclaimed by forest habitats after the retreat of the continental ice sheet was populated from two long-term separated refugia, presumably located in southeastern Russia and western Europe. Similar diversification was identified in P. nigroruber ( Figure 1H, cox1-5 divergence 4.9% in Finland). Still, it stands at odds with the genetically homogenous Fennoscandian population of the other two North European net-winged beetles, P. minuta, and L. sanguineus ( Figure 1F,G). The dual colonization of Fennoscandia has been reported in some plants, mosses, and rodents [103][104][105]. Our findings open the questions of a possible reproductive isolation of sympatric Fennoscandian populations (mtDNA divergence up to 6.8%; Table S1) and the frequency with which various dispersal routes were employed during the recolonization of Northern Europe after the LGM and which species were able to retreat only to some refugia. Further research must be based on more extensive data than we currently have at our disposal. We need to densely sample the populations from the whole range including Eastern Europe and Siberia and to employ multiple nuclear markers to robustly recover the colonization routes after the retreat of the continental ice sheet.

Conclusions
The high genetic diversification within the traditionally delimited morphospecies D. aurora exemplifies unlinked genetic and morphological divergence. The Holarctic distribution of D. aurora has been unique as no other net-winged beetle species is simultaneously known from two zoogeographic regions [17]. Regardless of morphological similarity, we have to resurrect the name D. coccinata for Nearctic populations. As we do not have deep insight into the structure of the whole genome, we cannot robustly identify the processes that control morphological uniformity. Provisionally, we assume the role of the stabilizing selection in the Mullerian mimicry ring (besides the possible intrinsic morphological stasis). The broad ecological valence of D. aurora questions the recent hypothesis linking the stable environment and morphological uniformity of net-winged beetles. The present results show the potential of poorly dispersing beetles for studies on the origin of the large regional faunas and the distribution dynamics since the LGM. With growing information content, the barcode database (www.boldsystems.org, accessed on 1 June 2022) becomes a powerful tool for taxonomy. Yet, the barcodes provide only the first clue. Therefore, we separate as full species only Nearctic and Palearctic populations for which we have additional data. The robust reconstruction of dispersal routes and eventually delimitation of further intraspecific taxa will need nuclear markers and a dense, now unattainable, sampling from the whole range.

Data Availability Statement:
The data presented in this study are available in the article.