The specimens from the historical collection of the Zoological Institute of the Russian Academy of Sciences, St Petersburg, Russia (ZISP) and recently collected material were examined. Type specimens of Stenodema spp. retained in the Finnish Museum of Natural History (MZH) were also studied. The specimens were initially identified using the keys published in Kerzhner and Jaczewski (1964), Wagner (1974), Vinokurov and Kanyukova (1995), and Yasunaga (2019). The following number of specimens were examined for this study: Stenodema calcarata (71), S. holsata (46), S. laevigata (52), S. rubrinervis (12), S. pilosa (13), S. sibirica (50), S. trispinosa (64), S. turanica (41) and S. virens (39). The collection event data for all of them were entered to the Arthropod Easy Capture Database (https://research.amnh.org/pbi/locality/index.php) and available through the Heteroptera Species Pages (https://research.amnh.org/pbi/heteropteraspeciespage/speciesdetails.php). All specimens were examined externally, and at least 10 males and 10 females from different series for each species were dissected for examination of the genitalia. The list of non-type specimens examined for this study is provided in Suppl. material 1.

For the molecular studies, the specimens from the following species were used: S. calcarata (13 specimens). S. holsata (4 specimens), S. laevigata (11 specimens), S. trispinosa (3 specimens), S. turanica (3 specimens), S. virens (3 specimens), Leptopterna dolobrata (Linneaus, 1758) (1 specimen) and Trigonotylus sp. (1 specimen). The genitalia structures were examined for all Stenodema vouchers.

To examine the male and female genitalia structures, abdomens were removed and boiled in 10% KOH for up to five minutes and dissected in water. Afterward, the abdomens were stored in glycerol. In some cases, aedeagi were inflated after this procedure. Aedeagi were also inflated using 40% lactic acid, following the detailed procedure described in Namyatova et al. (2021). The drawings were completed using Leica DM2500 microscope with the drawing device attached. The terminology of genitalia follows Konstantinov (2000, 2003) for males and Schwartz (2008) for females.

The digital images were taken in stacks using the Canon EOS 5D Mark IV camera equipped with a Canon MP-E 65 mm f/2.8 1–5× Macro lens and a Twin-Lite MT-26EX-RT flash. Partially focused images were combined using the Helicon Focus software. The SEM images were taken from uncoated specimens using the Hitachi TM1000 tabletop microscope.

Measurements were completed using Micromed MS-5 microscope using a graticule and ×10 eyepiece. Measurements statistics is provided in Table 1. Scale bars for habitus images equal 1 mm, the scale bars for genitalia structures equal 0.1 mm. Measurements provided in the diagnoses and descriptions are in mm.

The DNA was extracted from abdomens of ethanol-stored and dry specimens using the Evrogen Extract DNA Blood and Cells kit. The standard protocol was used with two modifications. First, the abdomens were kept overnight in the lysis solution with proteinase K in the water bath. Second, 50 or 25 μl of elution buffer was added at the final stage to increase the DNA concentration. After lysis, the abdomens were kept in glycerol for further examination. To obtain the barcoding region of cytochrome c oxidase subunit I (COI) the primers from Vishnevskaya et al (2016) were used with the annealing temperature equaling 45 °C or 42 °C. To obtain 16S rRNA region, the primers from Menard et al. (2014) were used with the annealing temperature 48 °C. For both markers, temperature of the initial denaturation and denaturation was 94 °C (3 mins and 30 secs, respectively), and extension and final extension temperature was 68 °C (1 min and 10 mins, respectively). The PCR products were cleaned using Evrogen Clean-up S-Cap kits or with Exonuclease I Thermofisher and sequenced in Evrogen (https://evrogen.ru/). The products were between 647 to 847 for COI and between 361 to 403 for 16s rRNA. The base pairs were trimmed at both ends if they were absent in more than half of the sequences in the alignment. The sequences were uploaded to GenBank, the accession numbers are listed in the Suppl. material 2.

The sequence diversity was calculated using P-distance and Kimura-2-parameter (K2P) in MEGA-X (Tamura et al. 2021) within each species, between species and between the clades within species.

Alignments were completed using Geneious algorithm in Geneious v. 11 software for each marker separately. Alignments included 36 original COI and 16S rRNA each. The COI alignment also included 84 sequences downloaded from Genbank: S. calcarata (15), S. holsata (17), S. laevigata (17), S. pilosipes (2), S. rubrinervis (5), S. sericans (3), S. sibirica (4), S. trispinosa (15), S. vicina (5), S. virens (1). Alignment for 16s rRNA additionally included four sequences of S. rubrinervis (2) and S. sibirica (2) from GenBank. Both alignments included original sequences of two outgroup taxa, Leptopterna dolobrata and Trigonotylus sp. All GenBank accession numbers are listed in the Suppl. material 1. Two alignments were concatenated using Geneious. The alignment lengths for COI and 16s rRNA were 787 bp and 399 bp, respectively. Phylogenetic analyses were run on each marker separately and for the combined datasets. Both combined datasets were 1186 bp length. First of them included all sequences available and included 124 terminals (full dataset). The second dataset included 34 specimens for which both markers were obtained (reduced dataset). In all the cases, Trigonotylus sp. was chosen as a root.

Maximum Likelihood approach implemented in RAxML v. 8.2.12 (Stamatakis 2014) with 10000 bootstrap replicates (BS) was performed. The phylogenetic trees were also calculated using Bayesian inference with MrBayes v. 3.2.7 (Ronquist et al. 2012). The main settings for MrBayes included 20 million generations, four chains, and the burn-in was set at 25%. Posterior probabilities were used for the node support (PP). Log files were checked to ensure that the standard deviation of split frequencies reached 0.01. All analyses were run using the server Dell PowerEdge R7525 (Dell Inc., USA).

Automatic barcode gap discovery approach (ABGD) was used via the online tool (https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html) on the alignment of each marker separately. This algorithm searches for a gap, which can be observed whenever the divergence among organisms belonging to the same species is smaller than the divergence among organisms from different species (Puillandre et al. 2012). The P range was set at 0.001–0.01, and Kimura (K80) model was used to estimate the matrix of pairwise distances.

Poisson tree process model (PTP and bPTP) and Generalized Mixed Yule Coalescent approach (bGMYC) were applied to the phylogenies built on a single marker and on combined datasets. Both approaches model the transition in branch length between species in contrast to within species (e.g., Blair and Bryson 2017) as another indication of speciation events. GMYC is a model-based likelihood approach that combines phylogenetics and coalescence theory, was proposed to estimate species boundaries from DNA sequence data. This algorithm identifies the transition points between inter- and intra-species branching rates on a time-calibrated ultrametric tree by maximizing the likelihood score of the model (Pons et al. 2006; Reid and Carstens 2012; Fujisawa and Barraclough 2013). PTP approach does not need an ultrametric tree and model speciation rate by directly using the number of substitutions (Zhang et al. 2013).

For all analyses, bGMYC, PTP, and bPTP, only unique sequences were left in the datasets, because zero-length branches can affect the results (Reid and Carstens 2012). The duplicates were removed using the online tool sRNAtoolbox (Aparicio-Puerta et al. 2022) (https://arn.ugr.es/srnatoolbox/helper/removedup/). It is recommended to run the species delimitation analysis based on several trees, which helps to overcome the problems with the phylogenetic uncertainty, occurring when the species delimitation is applied for the single tree (Reid and Carstens 2012; da Silva et al. 2018). The trees were calculated using BEAST2 v. 2.6.3 software (Bouckaert et al. 2014) using GTR+G+I nucleotide substitution model with 50 mln chain length. The results were checked in Tracer v. 1.7.1.(Rambaut et al. 2018) to make sure that all parameters had effective sampling size exceeded 200, which is considered adequate for convergence (https://beast.community/analysing_beast_output). The LogCombiner application from the BEAST package was used to obtain the .tre file with ~ 100 trees for each case.

Species delimitation using GMYC was run in R with the bGMYC package with the parameters recommended in the instructions (http://nreid.github.io/assets/bGMYC_instructions_14.03.12.txt), the multiple thresholds was used, MCMC equaled 50000, and thinning equaled 40000. This analysis provides the list of all possible species, and we have chosen the set of species with the highest mean supports.

Bayesian and Maximum Likelihood implementations of the Poisson tree process model (PTP and bPTP) (Zhang et al. 2013) using the scripts in Python (https://github.com/zhangjiajie/PTP accessed in 31/10/2021) were used. The number of iterations equaled 100000. All analyses were run using the server Dell PowerEdge R7525 (Dell Inc., USA).

Bayesian Phylogenetics and Phylogeography (BPP) method tests species using the multispecies coalescent model (Yang 2015), and it was applied to the combined datasets, which includes both, COI and 16s rRNA. It tests whether the separated species has higher supports than the clade comprising combination of species. The specimens should be preliminary assigned to a putative species for this analysis. For each dataset, the specimens were assigned to species based on the phylogenetic results and the bGMYC, PTP, and bPTP analyses ran on the corresponding dataset. The root was removed from the datasets. The analysis was run through the interface version for Windows (https://abacus.gene.ucl.ac.uk/software.html). The A11 (species delimitation and species tree) analysis with nsamples = 50000, sampfreq = 2, burnin = 25000 was applied. All other settings were default.

Our study showed that most of the widely distributed Palearctic species can be separated from each other using external characters, as well as male and female genitalia. The diagnoses for those species are provided in this section.

Below we provide the key to species, where we included all widely distributed Palearctic species. We also added S. algoviensis Schmidt, 1934 (Central Europe), S. alpestris Reuter, 1904 (China), S. chinensis Reuter, 1904 (China), S. crassipes Kiritshenko, 1931 (Central Asia), S. khenteica Muminov, 1989 (Mongolia), S. plebeja Reuter, 1904 (China), S. rubrinervis Horváth, 1905 (China, Korea, and Japan), and S. sericans (Fieber, 1861) (Europe) to this key, because we had an opportunity to examine them. Stenodema nippon Yasunaga, 2019 was included, as Yasunaga (2019) provided a detailed illustrated description for this species. Thus, the key is designed to discriminate all Stenodema spp. of the Western Palearctic, Siberia, and the Far East. However, it does not include 16 of the 19 species originally described and currently known only from China. For a taxonomic account of Chinese species of Stenodema, refer to Zheng et al. (2004). Species comparisons are provided following the diagnoses.

Based on the descriptions and material examined, we could delimit five morphological groups within Stenodema.

1. S. calcarata - pilosa group (subgenus Brachystira). This group has the frons not protruding above the clypeus (Fig. 1I) and hind femur possessing ventroapical spines and not tapering towards apex (Fig. 2A, D). The information on S. falki is very scarce, but this species might also belong to this group. According to Kelton (1961), the Nearctic species S. falki is very similar to S. pilosa, but differs in body ratios and male genitalia, although Schwartz (1987) suspected that these species might be synonymous. Among the species with the female genitalia examined, S. calcarata and S. pilosa are similar in the absence of the membranous swelling on the dorsal labiate plate and the absence of the dorsal structure between interramal lobes (Fig. 4A, C, E, H).
2. S. holsata group includes species with the frons not protruding above the clypeus (Fig. 1H) and hind femur without spines, and non-tapered apical region (Fig. 2E, H). Stenodema algoviensis, S. chinensis, S. holsata, S. plebeja, and S. sericans possess this set of characters. Stenodema chinensis differs from the other four species with the presence of the flattened dorsal setae. Stenodema plebeja is longer and differs from other species in the body length/pronotum width ratio equaling 4.9–5.0 in females, while this ratio is < 4.4 in other species. In contrast to other species, Stenodema sericans is pale, without dark stripes on pronotum and hemelytron, and has parameres different from S. algoviensis and S. holsata, with apical half of the right paramere as wide as basal part, and the left paramere without outgrowth or swelling near the apical process (Wagner 1974: fig. 9d–f). Refer to the notes section after the diagnosis of S. holsata for the differences between S. algoviensis and S. holsata.
3. S. laevigata group includes species with frons not protruding above clypeus (Fig. 1H) and hind femora without spines and tapered apical region (Fig. 2E). According to Zheng (1981a), S. antennata Zheng, 1981 is close to S. laevigata, but much larger, with a female body length of 11. Stenodema longula Zheng, 1981 might be close to S. laevigata as well (Zheng, 1981a), although this requires further verification.
4. S. turanica - virens group includes species with frons protruding above the clypeus base (Fig. 2C) and hind femur lacking spines and apical region tapered (Fig. 2B, C). Stenodema crassipes is close to S. virens and S. turanica. However, it differs from them in the widened hind femur, which is 4–5× as long as wide, and the antennal segment II in female is widened basally with long and dense setae. Based on the drawings of the head and hind tibia in Zheng (1981a), S. tibeta Zheng, 1981 also belongs to this group. Nonnaizab and Jorigtoo (1994) compared S. deserta Nonnaizab & Jorigtoo, 1994 with S. virens, and noted that the former was different in the body structure and the paramere shape. However, those differences can be intraspecific variability, because according to our examinations, the parameres and vesica of S. virens are very similar to those depicted in Nonnaizab and Jorigtoo (1994: figs 4–6, 7G–I), and those two species could be conspecific. According to Zheng (1981b), S. hsiaoi is similar to S. virens and S. turanica in the habitus and male genitalia. The same is true for Stenodema mongolica Nonnaizab & Jorigtoo, 1994 (Nonnaizab and Jorigtoo 1994: figs 7–12); however, according to the original description it has flattened setae on the antennal segment I. The Nearctic species S. vicina (Provancher, 1872), S. imperii Bliven, 1858, S. sequoia Bliven, 1955, and S. pilosipes Kelton, 1961 are allied to the species of virens- turanica group (Bliven 1955, 1958; Kelton 1961).
5. S. sibirica group includes species with the frons protruding above the clypeus base (as in Fig. 1C), its hind femur does not have spines, and it is not tapering towards apex (Fig. 2F). Stenodema nippon is very similar to S. sibirica and S. rubrinervis, although distinctly differs from them in the salient features and genital structures (Yasunaga 2019). Stenodema khenteica is also within this group and differs from S. nippon, S. sibirica, and S. rubrinervis in the antennal segment I shorter than pronotum and distinctly narrower than the eye diameter. Many other species described from China, most probably, belong to this group, and some of them might be conspecific with the species listed above. We had an opportunity to examine the paralectotype of S. alpestris Reuter, 1904 preserved at ZISP. In salient features and measurements this species is identical with S. rubrinervis. Horváth (1905) did not compare this species with S. alpestris, and possibly he was not aware of it. The lectotypes of both species should be examined to draw conclusion on their status. Stenodema gridellii Hoberlandt, 1960 has similar parameres to S. sibirica, but it has a smaller body (Hoberlandt 1960). Although in the drawings of Hoberlandt (1960)S. gridellii is shorter than S. sibirica, the provided measurements for the former fit those for S. sibiricaZheng (1981a) compared S. alticola Zheng, 1981 with S. gridellii, but wrote that the former is longer (males 6.7–6.8, female 7.5–7.6), having the erect setae on antennal segment I and its antennal segment II/I length ratio was 2.5. All those characters correspond to S. rubrinervis. Therefore, S. alticola and S. rubrinervis can be closely related or even synonymous. According to Zheng (1981a), S. nigricalla Zheng, 1981 is similar to S. chinensis (from the turanica- virens group). Judging from the drawings of the male genitalia, the shape of vesica of S. nigricalla is more similar to specimens from Siberia with short ridge on vesica (see Notes for the diagnosis of S. sibirica; Zheng 1981a: figs 13D–F, 19). However, the right paramere of this species has a longer apical process, than in those specimens and it is more similar to S. rubrinervis (Zheng 1981a: fig. 18; Yasunaga 2019: fig. 8a, e). Zheng (1981a) compared S. angustata Zheng, 1981 with S. nigricalla. However, the latter is more similar to S. nippon in the shape of the right paramere with elongate apical process and the presence of long and narrow vesica lobe at the left hand side (Yasunaga 2019: fig. 7C, F, G; Zheng 1981a: figs 20, 22). Tang (1994) compared S. qulininginensis Tang, 1994 with S. nigricalla, and, most probably, it also belongs to the sibirica group. Zheng (1992) noted that S. daliensis Zheng, 1992 is similar to S. alticola and S. gridellii and differed from them in the body shape and coloration. According to Reuter (1904)S. elegans has the hind femur without spines and not tapering apically and its frons is protruded above clypeus, which also corresponds to the sibirica group.

We could not place S dorsalis (Say, 1832) and S. parvula Zheng, 1981 to any group listed above. Kelton (1961) proposed to treat S. dorsalis described from the Eastern USA as nomen nudum, because there were no records for it since the original description. Stenodema parvula could be close to S. holsata or S. laevigata because its frons does not protrude above the clypeus; however, the information on the hind femur shape or genitalia structures for this species were not provided (Zheng 1981a).

The resulted trees from the Bayesian analyses are provided in Figs 14–19, and those resulted from the RaxML analyses are provided in Suppl. material 3.

All analyses show that widely distributed Palearctic species are monophyletic with high supports, as well as S. rubrinervis and S. sericans (P = 100, BS > 92). The COI sequences of two species from Nearctic, S. pilosipes and S. vicina, were included in the analyses. Stenodema vicina forms a clade (PP = 96 and 94, BS = 89 and 87 for COI and full datasets, respectively). However, S. pilosipes forms a clade only in the Bayesian analysis based on the full dataset (PP = 53), and in other cases one of the specimens is closer to S. vicina rather than to the second specimen of its species. Stenodema calcarata and S. pilosa always form sister group relationships (PP = 94–100, BS = 66–83).

The topologies built on COI only and the full dataset comprise the greatest number of specimens and species, and they are very similar. They show that the clade formed by S. calcarata and S. pilosa (subgenus Brachystira) forms sister group relationships with the clade comprising all other Stenodema species (nominative subgenus), and the latter has the following supports: PP = 100 and 87, BS = 67 and 70 for COI and full datasets, respectively. Within this clade, S. holsata, S. laevigata, and S. sericans form a clade (PP = 100 and 89, BS = 82 and 75 for COI and full datasets, respectively). In the analyses based on COI only, the relationships between those three species are unresolved. However, in the phylogeny based on the full dataset, S. sericans forms a clade with S. laevigata although with low supports (PP = 63, BS = 55). Stenodema pilosipes, S. sibirica, S. rubrinervis, S. turanica, S. vicina, and S. virens form a clade (PP = 97 and 100, BS = 77 and 91 for COI and full datasets, respectively). Among those species, S. turanica and S. virens are sister groups (PP = 100 and 99, BS = 80 and 87 for COI and full datasets, respectively), and S. pilosipes and S. vicina also form a clade (PP = 100, BS = 98 in both analyses). Those two pairs show reciprocal monophyly (PP = 100 and 99, BS = 77 and 89 for COI and full datasets, respectively). Stenodema sibirica and S. rubrinervis form a clade in Bayesian analysis (PP = 100 and BS = 62 full dataset), and in the RaxML analysis with COI and 16S rRNA (BS = 90).

The phylogeny based on the reduced dataset with COI and 16S rRNA has the topology corresponding to those obtained based on COI and full datasets.

The results obtained with 16S rRNA have a different topology. In this case, S. turanica forms sister group relationships with the clade comprising other species, although with low support (PP = 67). Stenodema virens forms sister group relationships with the rest of Stenodema species (PP = 89). Stenodema sibirica and S. rubrinervis form a clade (PP = 97, BS = 85), which is a sister group to the clade, formed by S. calcarata, S. holsata, S. laevigata, and S. pilosa (PP = 86). Stenodema laevigata is a sister group to a clade formed by other three species (PP = 85, BS = 68). Stenodema holsata is a sister group to a S. calcarata+S. pilosa clade (PP = 94, BS = 83).

At least some analyses show genetic structure within S. calcarata, S. pilosa, S. holsata, and S. laevigata. Analyses based on 16S rRNA and reduced dataset do not show the structure within S. pilosa and S. laevigata.

The phylogenetic structure within Stenodema pilosa is present only in the results of analyses based on COI and full datasets because Nearctic species are included there. The specimens of this species are split into three main clades: two Nearctic and one Palearctic. One of the Nearctic clades (PP = 100 for both, BS = 88 and 84 or COI and full datasets, respectively) is a sister group to the rest of the specimens. The clade comprising some Nearctic and all Palearctic specimens has low to average supports (PP = 88 and 82, BS = 55 and 62 for COI and full datasets, respectively). This clade splits into two groups: one of them Nearctic (PP = 100 for both, BS = 100 and 99 for COI and full datasets, respectively), and the second one is Palearctic (PP = 100 and 92, BS = 99 and 94 for COI and full datasets, respectively).

In the analyses based on COI and full dataset, representatives of S. calcarata from the southern side of Caucasus (Iran, Georgia, Turkey) and a single specimen from Germany form a clade with the highest support, and it is a sister group to the rest of the specimens of this species (PP = 100 and 97, BS = 90 and 71 for COI and full datasets, respectively). Specimens from East Asia (South Korea and Primorsky Territory) form a clade with the highest support, which is a sister group to the clade formed by the rest of the specimens (PP = 100 and 92, BS = 100 and 90 for COI and full datasets, respectively). Only 16S rRNA was obtained for the specimen from Stavropol Territory, and in the phylogeny based on the full dataset it is a sister group to the rest of the specimens (PP = 97, BS = 95). Two specimens from Germany form a clade (PP = 100 and 99 for COI and full datasets, respectively, and BS = 100 for both), and they are the sister group to the clade comprising most of the specimens from the Western Palearctic and a specimen from Altay Republic (PP = 55 and 51 COI and full datasets, respectively, BS = 82 for COI).

In the phylogenies based on 16S rRNA and the reduced dataset, specimens of S. calcarata from Georgia and Turkey form a clade with the highest supports. In the phylogeny based on 16S rRNA and the reduced dataset, single specimen from the East Asia (Primorsky Territory) included in those analyses has many substitutions. In the analysis based on the 16S rRNA it forms unresolved relationships with the clade, comprising the specimens from Georgia and Turkey (PP = 100, BS = 100) and the clade comprising the rest of the specimens (PP = 95, BS = 83). In the phylogeny based on the reduced dataset, the clade comprising species from Georgia and Turkey forms a reciprocal monophyly with the clade comprising the rest of the specimens including the one from the Primorsky Territory (PP = 99, BS = 58). In the phylogeny based on 16S rRNA specimen from Stavropol Province forms a clade with the clade comprising most of the specimens from the Western Palearctic and Altay Republic (PP = 100, BS = 95 in both datasets).

In the phylogenies based on COI and full dataset all specimens of S. holsata from France form a clade (PP = 93 and 94, BS = 87 and 88 for COI and full datasets, respectively), and it is a sister group to the clade formed by the rest of the specimens in the results of the Bayesian analysis (PP = 56 and 51 for COI and full datasets, respectively). Specimen from Karachay-Cherkessia forms a clade with the clade formed by the specimens from Northern and Central Europe (PP = 53 and 64 for COI and full datasets, respectively, BS = 99 for COI dataset). Only four specimens of S. holsata are included to the analyses based on 16S rRNA and reduced dataset. The specimen from Karachay-Cherkessia is a sister group to a clade formed by three specimens from northern Europe (PP = 96 and 91, BS = 82 and 76 for COI and full datasets, respectively).

In the phylogenies based on COI and full dataset, there is a clade within S. laevigata comprising specimens from Greece, Iran, and Voronezh Province (PP = 100 and 82, BS = 98 and 59 for COI and full datasets, respectively). Within this clade, the specimens from Voronezh Province and Germany are more closely related (PP = 100 and 99 for COI and full datasets, respectively, BS = 100 for both datasets). The results of the analysis based on the full dataset does not show any other clades within this species. The Bayesian inference analysis based on COI dataset also show, that the rest of the specimens except for the three specimens mentioned above and one from Iran, also form a clade (PP = 85). Within this clade, a specimen from Crimea forms sister group relationships with the rest of the specimens (PP = 95).

All analyses show identical results for the phylogeny built based on 16S rRNA. In the case of COI, ABGD delimits the smallest number of species, followed by GMYC. PTP and bPTP show identical results for this marker. In the analyses based on the combined datasets, GMYC results in the smallest number of species. For the reduced dataset, PTP, bPTP, and BPP show identical results. For the full dataset, BPP results in the largest number of species, and PTP and bPTP showed the identical number of species. All species delimitation analyses do not mix the specimens belonging to different widespread species. Stenodema sibirica, S. turanica, and S. virens each form a single species in all the cases.

All analyses suggested that S. calcarata can be a complex of at least three species: (1) Far Eastern clade (2) West Asian clade and a single specimen from Germany, (3) Euro-Siberian clade. Additionally, specimen from Stavropol Province, a clade with two specimens from Germany and specimen from Germany in the West Asian clade form separate clades in some analyses.

Stenodema pilosa also can be a species complex. In the analyses with Nearctic specimens (COI dataset and full dataset) the Palearctic representatives of this species are placed in a single species, and Nearctic sequences are grouped in two or three species.

Stenodema laevigata was subdivided into different number of species depending on the analysis. All analyses based on 16S rRNA, ABGD analysis based on COI and GMYC analysis based on the reduced dataset with both markers, and GMYC, PTP and bPTP analyses for the full dataset place all representatives of this species together. Specimen from Crimea, specimen from Iran and the clade formed by the specimens from Voronezh Province, Greece and Germany are assigned in separated species each by some analyses. Additionally, the specimens from Greece also appeared as a separate species in few cases.

Analyses based on COI and full dataset result in three species within S. holsata: (1) all specimens from France, (2) specimen from Karachay-Cherkessia, (3) specimens from northern and Central European areas. Only four specimens (one from Karachay-Cherkessia and three from northern European areas) are included in the analyses based on 16S rRNA and reduced dataset. The analysis based on the reduced dataset shows that the specimen from Karachay-Cherkessia forms a separate species, the analysis based on 16S rRNA places all specimens of S. holsata into a single species.

Interspecific distances are 6–17% for COI and 5–12% for 16S rRNA, and they are provided in Suppl. material 4. The lowest distances for COI are between S. turanica and S. virens (6–7%). Those two species also have relatively low genetic differences with Nearctic species S. pilosipes and S. vicina (~ 7–8%). The highest distances for COI are between S. pilosa and S. rubrinervis (15–17%). For 16S rRNA, the distances between the two species pairs S. sibirica — S. ribrinerve and S. turanica — S. virens are the lowest (~ 5%), and the highest distances are between S. pilosa and S. virens (11–12%).

For the COI analysis, seven, three and four specimens are included, respectively for S. sibirica, S. turanica and S. virens (Suppl. material 4). Although the specimens of S. sibirica and S. virens were collected in different regions (S. sibirica: from Altay to South Korea, S. virens from Finland, Caucasus, and Irkutsk Province), the diversity of their sequences is very low for both markers (< 0.12%). The COI sequences for S. turanica collected in Iran and Tyva Republic are identical. For 16S rRNA, a single specimen of S. turanica is included. There are five specimens of S. sibirica and two sequences of S. virens, and genetic distances within these species are < 0.1%.

Stenodema holsata and S. laevigata have within species mean distance corresponding to 0.8–1.1% for COI and ~ 0.4% for 16S rRNA. The species delimitation analyses resulted in three groups within S. holsata for COI, and the distances between them are 1–4%. The largest number of groups delimited within S. laevigata is five for COI, and the distances between them are 1–3%.

The interspecific distances within S. calcarata and S. pilosa are the largest, ~4% for COI for both species, ~ 2% for 16S rRNA of S. calcarata and 0.1% for 16S rRNA of S. pilosa. The largest number of species resulted from the species delimitation analyses for S. calcarata and S. pilosa using COI are five and four, respectively. The distances between the groups within S. calcarata are 7–9%, and between groups of S. pilosa are 2–6%. The species delimitation analysis based on the 16S rRNA dataset showed four groups within S. calcarata, and the distances between them are 3–4%.

The authors have declared that no competing interests exist.

No ethical statement was reported.

The study is supported by the Russian Foundation for Basic Research, project number 20-54-56011. We thank Heidi Viljanen (FMNH) for providing access to the collection at FMNH.

A.A. Namyatova and F.V. Konstantinov designed the study. A.A. Namyatova contributed to the laboratory work, completed all the analyses, morphological studies, SEM imaging, drawings and created the first draft of the manuscript. P.A. Dzhelali completed the measurements and tables. A.A. Namyatova and P.A. Dzhelali completed the plates and specimen databasing. F.V. Konstantinov and P.A. Dzhelali completed the photographs. All authors contributed to the manuscript editing.

Anna A. Namyatova https://orcid.org/0000-0001-9678-3430

Polina A. Dzhelali https://orcid.org/0000-0002-0741-3655

Fedor V. Konstantinov https://orcid.org/0000-0002-7013-5686

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