We studied about 130 adult specimens for this work. We evaluated 18 species and included 30 terminal taxa in the analyses (see the section on the phylogenetic integrity test), with 11 species comprising the ingroup, including the two new species; all species are listed in Table S2. We used the following species as outgroups: B. amazonum Estévez and Polhemus, 2001, B. denticolle Montandon, 1903, from B. denticolle group Lauck, 1962; B. minor (Pasilot de Beauvois, 1805), from B. minor group Lauck, 1962; B. candidulum Montandon, 1903, B. horvathi Montandon, 1903, B. oxyurum (Dufour, 1863) and B. sanctulum Montandon, 1903, from B. oxyurum group Lauck, 1962. We chose the outgroups based on the criteria of Nixon and Carpenter (1993), more precisely when taking into account the existence of some level of morphological proximity as suggested by Lauck (1962), Nieser (1975), Estévez and Polhemus (2001), Estévez and Polhemus (2007), and Ribeiro and Estévez (2009).

The material studied comes from the collections of the American Museum of Natural History, New York, United States of America (AMNH); epartamento de Parasitologia, Instituto de Ciências Biológicas, Universidade Federal de Minas Gerais, Minas Gerais, Brazil (DPIC); Departamento de Zoologia, Instituto de Biologia, Universidade Federal do Rio de Janeiro, Rio de Janeiro, Brazil (DZRJ); the Hungarian Natural History Museum, Budapest, Hungary (HNHM); Institut Royal des Sciences Naturelles de Belgique, Brussels, Belgium (ISNB); Instituto Nacional de Pesquisas da Amazônia, Manaus, Brazil (INPA); Museo Argentino de Ciências Naturales, Buenos Aires, Argentina (MACN); Museu Nacional, Universidade Federal do Rio de Janeiro, Rio de Janeiro, Brazil (MNRJ); Museum National d’Histoire Naturelle, Paris, France (MNHN); Museu de Zoologia da Universidade de São Paulo, São Paulo, Brazil (MZSP); and the Snow Entomological Museum, University of Kansas, Lawrence, United States of America (SEMC). All acronyms of depositories are based on Arnett-Jr. et al. (1993), except DZRJ and MNRJ. The list of examined material of the two new species is associated with the formal description of each one presented in the ‘Systematics’ section of our results.

To obtain the total body length and the greatest width of the body of adults, we used a manual calliper with an accuracy of 0.1 mm. We took all other measurements with the aid of a 0.01 mm accuracy millimeter ocular micrometer on a stereoscopic microscope. The measures taken are based on the measurement according to Ribeiro (2007) and Ribeiro and Alecrim (2008). We used all available materials of each species in this study to obtain the measurements. The referred measurements and relations are summarised as mean values in Table S3, and are as follows: TBL, total body length (from the apex of the head to the end of the abdomen); T/G, ratio between the total body length and the greatest width of body; R1/R2, ratio between the lengths of first and second visible segments of the rostrum; A/I, ratio between the median lengths of the anteoculus and interoculus; IW/A, ratio between the posterior interocular width and the median length of anteoculus; IW/O, ratio between the posterior interocular width and the width of an eye; W/L, ratio between the width and median length of the ventral diverticulum of the phallosoma.

We prepared all line drawings and sketches with a drawing tube mounted on a ZEISS (AXIOLAB) optic microscope and a ZEISS (SV-6) stereo microscope. We took photographs with a ZEISS stereo microscope with a MC 80 camera (with AXIOVISION 3.0 SP4 software), housed at Laboratório de Captação de Imagens do Instituto de Biologia, Universidade Federal do Rio de Janeiro (UFRJ), Rio de Janeiro, Brazil.

The terminology of the head, thorax, and abdomen follows Hamilton (1981), Lauck (1962), Nieser (1975), Parsons (1964, 1974), Estévez and Polhemus (2007), Ribeiro (2007), Ribeiro and Alecrim (2008), and Ribeiro and Estévez (2009) for specific areas, sclerites, and sutures of the head, thorax, and abdomen; Betts (1986) for basal portions of fore and hind wings; Gorb and Perez-Goodwyn (2003) for the Druckknopf system, which is responsible for coupling of the hemelytra and the body; and Estévez and Polhemus (2007) and Ribeiro and Estévez (2009) for male genitalia.

More specifically, the general terminology of the axillary sclerites used in this study follows Betts (1986) and Yoshizawa and Saigusa (2001). According to Betts (1986), there is a reduction in the number of axillary sclerites through Heteroptera, when compared with the possible ancestral state found in Homoptera s.l. However, the location and number of described sclerites remain constant through the representatives of Heteroptera, except for those of the genus Notonecta Linnaeus, 1758 (Notonectidae). Due to their easy visualisation and handling, we preferred to use the sclerites of the basal portion of the hind wing of the specimens studied after having certified that both wings have the same pattern regarding the distribution and shape of the sclerites. This verification agrees with the pattern found in the arrangement and shape of the axillary sclerites, as reported by Betts (1986) and Yoshizawa and Saigusa (2001). It could be possible to obtain hypotheses of homology (‘primary homologies’ sensu de Pinna 1991) for the sclerites found in the specimens studied, identifying the humeral plate; cubitus; and anal 1, 2, and 3 veins (hup, Cu, A1, A2, and A3, respectively, Fig. 1). By using these sclerites and veins as ‘markers’, we were able to establish all other units of comparison. In this study, we obtained six characters (Table S4, characters 21–26) based on the study of sclerite 1 of the distal median plate (‘DMP1’ for Yoshizawa and Saigusa 2001; dmp, Fig. 1), which was the only one that showed useful variation through the studied specimens.

Figure 1. 

Basal portion of the right hind wing and its axillary sclerites. Six characters were obtained based on the study of sclerite 1 of the distal median plate (characters 21–26), which was the only one to show useful variation through the studied specimens. By using these sclerites and veins as ‘markers’, we established all other units of comparison. These veins (and their branches) are named according to a system devised by John Comstock and George Needham: the Comstock–Needham system. Abbreviations: dmp, sclerite 1 of the distal median plate; hup, humeral plate. Veins: A₁, anal 1; A₂, anal 2; A₃, anal 3; C + Sc, costa + subcosta; Cu, cubitus; M_{1 + 2}, media 1 + media 2; R, radius.

The study of the Druckknopf system components, according to Weber (1930) and Schuh and Slater (1995), or the ‘knob-and-socket’ system, according to Gorb and Perez-Goodwyn (2003), which allows for the complete coupling of hemelytra to the dorsum, follows A. L. Estévez’s ideas (pers. comm.) regarding the purpose of homologies and coding characters. The system is known to occur in several families of aquatic Heteroptera, and the following authors have described it in detail: Cobben (1957) for Saldidae Amyot and Serville, 1843; Gorb and Perez-Goodwyn (2003) for Corixidae Leach, 1815, Naucoridae Fallèn, 1814, Nepidae and Notonectidae Leach, 1815; Larsén (1945) for Nepidae; and Perez-Goodwyn and Gorb (2003) for Belostomatidae. To study the knob-and-socket apparatus, it was necessary to carefully displace the hemelytra, by pulling the anterior portion along their lateral margins near to scutellum after moistening the insertion point of the hemelytra on the thorax. Concerning the procedure, it neither damaged the forewing nor completely decoupled the forewings from the body. Moreover, the involved musculature kept itself intact, allowing us to return the hemelytra to its initial position after the procedure. We have proposed characters by studying only the complementary structure of the apparatus, located in the thorax, and not that on the wing apparatus (Fig. 2).

Figure 2. 

Mesepimeron and part of wing-to-body coupling mechanism. Part of the ‘Druckknopf’ system component, located in the thorax. The system allows for the complete coupling of hemelytra to the dorsum. Character 27 was proposed by studying only the complementary structure of the apparatus, located in mesepimeron, not that on the wing apparatus. A Knobshaped ‘Druckknopf’ system rounded, as long as wide B knobshaped ‘Druckknopf’ system piriform, longer than wide.

The procedures for the preparation of male genitalia follow the same ones described by Ribeiro (2007); we adapted them for small insects. After the examination, we stored the dissected parts in microvials with glycerine, properly sealed, and then kept together with the voucher specimens.

Our proposition of homologies is based on Brower and Schawaroch (1996), in which the hypothesis of primary homology (de Pinna 1991) is an outcome of topographical identity between structures. Furthermore, based on Thiele’s (1993) ideas, we considered putative homologies from measurements and metric relations to be a particular distribution of frequency values in a sample of specimens of a taxon in this study. Such a distribution (or particular trend) is considered to be inheritable from its ancestor (Rae 1998). Forty-nine characters used in the phylogenetic analyses are listed in supplementary file 3, seven of which are continuous (Table S3, Table S4). We grouped the characters into the following categories: body size (1–2), head (3–8), thorax (9–27), abdomen (28–49), and male genitalia (30–49). Character coding was contingent for discrete data (Forey and Kitching 2000) and the character statements were formulated following the logical basis of Sereno (2007). We scored character states as dashes (-) if inapplicable and as question marks (?) if ambiguous or missing. We constructed the matrix of discrete characters in Mesquite v.3.61 (Maddison and Maddison 2019) and reconstructed character optimisations in MacClade (Maddison and Maddison 2005) using trees based on parsimony analyses with discrete and continuous characters combined. Only unambiguously optimised characters are discussed in a more detailed way through the text. We performed character optimisations on a morphological strict consensus tree from MP analyses with implicit weightings assigned to characters (analysis six, see below). Apart from MP optimisations, we optimised comparable characters based on male genitalia with Bayesian maximum credibility trees. When any character optimisation on trees was different, we have discussed such conflicts. For each character, we calculated the number of steps and the consistency (CI) and retention (RI) indices (Kluge and Farris 1969). Of the 42 discrete characters, 13 are multistate characters and have been coded as unordered (or ‘maximally connected’ in the terminology of Slowinski 1993) (Fitch 1971). Otherwise, we coded continuous characters as ordered (Farris 1970) and treated them as such, according to Goloboff et al. (2006). To be used as such, however, we ln-transformed the mean values of measurements and metric relations to reduce the supposed overdominance that they may present during implicit weighing analyses under MP (Koch et al. 2014).

We treated specimens with intraspecific variations (‘exemplars’, in the terminology of de Pinna 1999) as terminal taxa, receiving the word ‘var’ followed by a number, and we thus considered intraspecific variations as distinct character states. This procedure allowed us to test the phylogenetic integrity (de Pinna 1999) of some putative species represented by individual organisms. Following the protocol described by de Pinna (1999), we corroborated the phylogenetic integrity of a species containing individual organisms with intraspecific variation if individuals did not themselves constitute para- or polyphyletic groups or even monophyletic groups with subclades in a cladogram. The following taxa had specimens with the intraspecific variation being evaluated: B. amazonum, B. denticolle (with three variations), B. minor, B. candidulum, B. oxyurum, and B. sanctulum (outgroups); B. lariversi, B. micantulum, B. minusculum, B. nicaeum, and B. pygmeum (ingroup). We have only discussed variation found among representatives of the ingroup. Especially regarding B. micantulum, studied specimens may apparently constitute distinct populations, and this seems to be the case here because of the existence of variation in the shape of their arms of phallotheca. Besides, many authors had cited this species as actually a species complex (see Stål 1860; Nieser 1975; Ribeiro 2007;). In this study, we partly adapted the original method by using in some cases groups of individuals rather than only one specimen per putative species. The species treated in this way were B. lariversi, B. micantulum, B. minusculum, B. nicaeum, and B. pygmeum. If we could refute the phylogenetic integrity of some species-candidate assemblage, then we mention the species name inside quotation marks throughout the text.

In B. candidulum var2, the character referring to the apex of the prosternal carina (character 12) could not be coded, as the apex of carina was damaged. In B. sanctulum var2, male genitalia characters (30–49) could not be coded, as unfortunately we only studied females (see below).

Apart from non-ambiguous character optimisation in MP trees, we optimised characters 31, 32, 33, 43, 44, 45, 46, 47, and 48, which are discrete characters of male genitalia considered to be important to the taxonomy of B. plebejum group, by using both MP and ML procedures over the maximum credibility tree obtained by BI. With this additional procedure of optimisation via ML in the topology obtained by BI, we also sought to explore the veracity of the reconstructed ancestral states via MP (Yang 2014). We chose these characters mainly because they can provide a general panorama about the evolution of some male genitalia components likely to be known to define or even characterise the species of B. plebejum group according to Lauck’s (1962), Nieser’s (1975), and Estévez and Polhemus’ (2007) papers. Besides, these characters are considered to be comparable among all taxa in the analysis, therefore not evaluating the presence/absence of the structure and its aspect simultaneously.

The nine discrete characters had previously been modelled using the Mk model (Lewis 2001), assuming that transitions among states follow a Markov process. However, we also tentatively fitted alternative models for multistate characters, taking into account that the transition between all states of these characters did not necessarily occur at the same rate of evolution (extended Mk models according to Harmon 2019). We allowed these characters to evolve according to (1) an equal rate model (ER), the simplest model, which assumes that the rate of transition from one state to another is the same whatever the direction, as well as a stationary distribution of character states (Harmon 2019); (2) a multi-rate symmetric ordered model (ordSYM), similar to ER, but allowing pairs of states to differ from any other pair and assuming no stationary distribution of character states; (3) a model similar to ordSYM, now assuming that the forward transition rate is twice the backward rate (symORD2); (4) a model that assumes transition only occurring between adjacent states (ORD); and (5) an all-rates-different model (Paradis et al. 2004)—the most complex model because it has the highest number of free parameters (p = 3)—assuming every possible type of transition has a different rate. In sum, we also aimed to test for differences between the forward and backward rates of character change.

We assumed that each state of the aforementioned characters could occur at the root of phylogeny with equal probability (flat prior according to Harmon 2019), instead of treating it as an estimated nuisance parameter being fitting together with the other free parameters as suggested by Fitzjohn et al. (2009). Besides, we considered the characters to be independent of each other, so those different models could be fitted to the evolution of each of these characters separately. We used ML to find a better estimate of free rate parameters, fitting the statistical models to the dataset. We used the small sample size correction of the Akaike information criterion (AICc) to compare models and to identify those providing the most efficient way to describe patterns in our data with few parameters. Accordingly, we also estimated ΔAICc values, a relative distance between AICc scores of the proposed models, choosing the model with the smallest value, and the relative support for each model using Akaike weights (ω) (Burnham and Anderson 2002).

A summary of the results obtained by different combinations of weightings and characters is shown in Table 1. The main morphological structures and characters, as well as strict consensus trees, are indicated in Fig. 4C–E and Fig. 4H–J (male genitalia), Fig. 5 (unambiguous optimisation of discrete characters over the tree), Fig. 6 (sclerite 1 of distal median plate), Fig. 7 (prosternal carinae), and Fig. 8 (optimisation of continuous characters over the tree). In all analyses, the B. plebejum group is monophyletic (BS = 95%; BRel = 69%), whereas the B. plebejum group s.s., as conceived by D. Lauck, was recovered as paraphyletic in relation to the B. pygmeum group (Fig. 5). Of the six synapomorphies supporting the B. plebejum group, four are non-homoplastic and one of them refers to the phallosoma with the diverticulum strongly curved ventrally (character 31: 1). It is also worth mentioning a general trend in a decreased body length through its members (character 1, Fig. 8). The B. pygmeum group sensu Lauck (1962) was also recovered as a monophyletic group in almost all analyses, except in the scheme DISC + EW (analysis one), despite its low support (BS = 55%; BRel = 52%). The clade is supported by three synapomorphies—all of them based on the morphology of male genitalia—two of which are non-homoplastic: margins of the diverticulum of the phallosoma clearly thickened (character 32: 1); apical margin of the diverticulum conspicuously acute, in a lateral view (character 48: 1). Considering continuous characters, there was a general enlargement of male genitalia through the members of this group, which is evidenced by the increase in the ratio of the width to the length of the diverticulum (character 30) (Fig. 8). Based on the obtained topologies found by MP analyses of the data matrix consisting of different combinations of matrices of discrete and continuous characters, the clade composed of all species of the B. pygmeum group remained consistent throughout them (Fig. 8). Besides, the internal relation between B. lariversi and B. pygmeum was recovered across almost all schemes combining character matrices with different weights, except in the consensus tree of analysis one (i.e., without continuous characters and under EW; supported only by BS: BS = 81%, BRel = 16%). Our findings support such a clade by the presence of a distinct lobed and robust process, situated along the apical portion of sclerite 1 of the distal median plate (as in Fig. 6A) as a homoplastic synapomorphy.

Figure 5. 

Maximum parsimony inference from analysis six (all datasets under implied weights, K = 3). Consensus of four most parsimonious trees (L = 142.3, CI = 0.44, RI = 0.74, fit = 14.61) with non-ambiguous discrete characters optimised. Apart from non-ambiguous character optimisation, the discrete characters based on male genitalia and their states are indicated on male genitalia sketches and are coloured accordingly. Green character labels and their optimisations refer to the characters based on male genitalia. Numbers inside circles refer to character numbers and coloured circles indicate the states of each aforementioned character. Clades are coloured according to traditional classification: Belostoma plebejum group (black clades) and B. pygmeum group (light blue clades). The coloured names refer to the taxa whose phylogenetic integrity was not corroborated. Coloured squares refer to consistency index values. Outgroups were omitted. Abbreviation: CI, consistency index.

Figure 6. 

Sclerite 1 of the distal median plate. Different aspects of the sclerite: form, degree of development of process from posterior portion articulated with A₁ and A₂ veins, degree of development and form of apical lobed process, shape of basal margin of process from posterior portion articulated with A₁ and A₂ veins, and form of the basilateral margin of the process from posterior portion articulated with A1 and A2 veins. Six characters were obtained based on the study of sclerite 1 of the distal median plate (characters 21–26). A Belostoma denticolle Montandon var1 B B. minor (Pasilot de Beauvois) C B. estevezae Ribeiro and Alecrim D B. oxyurum (Dufour) E B. sanctulum Montandon F B. micantulum (Stål) var1 and var2 G B. parvum Estévez and Polhemus H B. plebejum (Stål).

Conversely, almost all analyses recovered the clade B. plebejum + B. minusculum with low support (BS < 50%; BRel = 15%), except in analyses one (DISC + EW) and four (DISC + CONT + EW) (Figs 5 and 8). Two homoplastic synapomorphies supported the clade, including the presence of medial and longitudinal depression or sulcus poorly developed on the diverticulum of phallosoma (Fig. 5). Regarding continuous characters, there is either a slight increase in the distance between eyes or a decrease in the size of eyes, which is indicated by these species sharing the increase in the posterior interocular width/width of an eye ratio values (character 7, Fig. 8). Accordingly, such a clade was not consistent across analyses when that condition was combined, in some way, with any of the continuous characters 1, 2, 3, 4, 6, and 30 (Table S8).

Figure 7. 

Lateral view of part of the head, rostrum, prothorax with the prosternal carina, and forelegs. Characters were obtained based on the different aspects of the prosternal carina in lateral view: form and shape, degree of enlargement taking into account its base, positioning relative to the longitudinal axis of the body, form of the apex and anterior margin, presence of tubercles on the apex. A Belostoma amazonum Estévez and Polhemus var2 B B. horvathi Montandon C B. oxyurum (Dufour) D B. minor (Pasilot de Beauvois) E B. nicaeum Estévez and Polhemus F B. pygmeum (Dufour).

Figure 8. 

Maximum parsimony inference from analysis six (all datasets under implied weights, K = 3). Consensus of the four most parsimonious trees (L = 142.3, CI = 0.44, RI = 0.74, fit = 14.61) with continuous characters optimised. A The continuous characters and their states are indicated by coloured rectangles on the left-top corner and are coloured according to their standardised values. Measurements based on male genitalia are indicated on male genitalia sketches B Codes after ln-transformed values referring to the ancestral state reconstructions presented on nodes were used in order to reduce space. Numbers inside rectangles above or below branches refer to bootstrap (if > 50%) with a Poisson correction and relative Bremer supports, respectively. Coloured squares refer to consistency index values. Letters inside grey squares refer to nodes that were also recovered in all bootstrap resamples C A line graph illustrates trends in bootstrap support of those nodes, as long as combinations of continuous characters were used together with discrete characters. Black and red rectangles indicate the existence of a node recovered by a maximum parsimony analysis and with a combination of continuous characters (indicated by the Subtree Pruning Regrafting distance), respectively. Outgroups were omitted. Abbreviations: BS, bootstrap support; CI, consistency index; C1, continuous character 1; C2, continuous character 2; C3, continuous character 3; C4, continuous character 4; C6, continuous character 6; C7, continuous character 7; C30, continuous character 30; lvd, medial length of the diverticulum; SPR differences, the Subtree Pruning Regrafting distance; wvd, largest width of the diverticulum.

Analyses three (DISC + IW_est) and six (DISC + CONT + IW_est) provided low but some support for two topologies (Table S5, Table S6): an overarching clade comprising B. lanemeloi sp. nov., B. nessimiani, and B. nieseri sp. nov. (BS < 50%; BRel = 21%), as well as an internal relation between the last two species (BS < 50%; BRel = 21%) (Figs 5 and 8). The first clade is supported by the presence of slightly longer arms of the phallotheca compared with other species, a non-homoplastic synapomorphy, so this condition renders to such arms the possibility of clearly reaching the posterolateral margins of the diverticulum (character 35: 0). There is also a slight decrease in the ratio of the anteocular to interocular medial lengths (character 4). Accordingly, the monophyly of that clade was not recovered when only that metric condition was combined, in some way, with any of the continuous characters 1, 2, and 7 (Table S8), such that the resulting topologies showed differences of six SPR units from the strict consensus tree generated by the whole dataset with IW_est (i.e., analysis six). Likewise, the second clade recovering B. nessimiani as the sister of B. nieseri sp. nov. is supported by the arms of the phallotheca clearly covering the lateral margins of the diverticulum of the phallosoma (character 41: 1), also a non-homoplastic synapomorphy (Fig. 5). There is also both an increase in the posterior interocular width/the anteocular length ratio values and a decrease in the ratio between this width and the width of an eye (characters 6 and 7, respectively). Finally, the continuous characters 6 and 7 provided further support for the aforementioned clade whenever they were combined in the matrices (Table S8).

The seven continuous characters herein used allowed us to produce 127 different combinations of continuous characters (from one-to-one to seven-to-seven combinations), each of them having been submitted together with the matrix of discrete characters to MP analyses. Of these, MP analyses using 87 combinations of continuous characters together with discrete characters exclusively found strict consensus topologies that did not differ from that topology found using only the 42 discrete characters with K-value estimated from the dataset (analysis three: DISC + IW_est); the SPR distances were equal to zero. On the contrary, 36 combinations produced together with discrete character topologies revealed two units of SPR distance, and only four combinations produced topologies with six units of SPR distance when compared to the strict consensus topology found by MP analyses when using only discrete characters (Table S8).

The measurements that somehow compose the four combinations of characters that produced the worst topologies (i.e., topologies with the greatest distances of SPR) are total body length (character 1), the ratio between the total body length and the greatest width of body (character 2), the ratio between the first and second visible segments of the rostrum (character 3), the ratio between the median lengths of the anteoculus and interoculus (character 4), the ratio between the anterior interocular width and the width of an eye (character 6), and the ratio between the width of the diverticulum and its length in ventral view (character 30). The combination with the highest number of continuous characters concomitantly inducing the production of strict consensus topologies with the highest conflict with that topology found only with the partition of discrete characters (i.e., SPR distance equals to six) is the one that uses all the aforementioned measures, except the ratio between the median lengths of the anteoculus and interoculus.

The progressive increment in the number of continuous characters during MP analyses necessarily increased the SPR distances between strict consensus topology found with exclusive use of the partition of discrete characters and the strict consensus topologies found with the use of continuous and discrete partitions (Ordinary Least Squares [OLS] regression: b = 0.08; p = 0.0187). However, a proxy of node stability—herein measured by nonparametric bootstrap—was not dependent on SPR distance values, when scaling the number of continuous characters added to bootstrap values of the obtained clades through a strict consensus tree. Accordingly, the relationship between these two variables was not significant, so the estimated slope does not suggest a scenario that the increase in SPR distances shows some degree of covariation with the decrease in node supports (OLS regression: b = 0.01; p = 0.7831). Only the decrease in retention index values showed some significant degree of covariation with the increment in number of continuous characters (OLS regression: b = –0.01; p < 0.001).

In general, the conflict between the two partitions of characters seems to be weak because the clades found by using character matrices based on measurements and relationships, including their combinations, were consistent with those found by adopting only the partition of discrete characters. More specifically, when comparing the topologies found in analyses one, two, three, and six (under MP), it is clear that the use of continuous characters in this study does not seem to have rescued different relationships among taxa (Fig. 5). This may indicate the presence of some congruence between these partitions of characters, or could have arisen as a kind of sampling behaviour where such a small number of characters—as herein adopted with the use of only seven continuous characters—had no significant effect on the phylogenetic signal present through another partition with an evident higher number of characters (42 discrete characters). Besides, it should be highlighted that the clades found in the strict consensus tree of analysis three are practically the same as found in the strict consensus tree with the use of discrete and continuous characters (Fig. 8). However, the greatest SPR distances between the strict consensus tree of analysis three and the topologies found with combinations of continuous and discrete characters were not obtained with the use of all continuous characters activated, such as topologies found in analysis six, or with these characters each used separately. Furthermore, there was no general increase in overall support with the progressive increment of continuous characters in characters matrices (Fig. 8), but instead—except for the clade comprising B. plebejum and B. minusculum, as well as the clade comprising B. lariversi and B. pygmeum—there was a slight increase in node supports, especially those independently formed by the variations of B. minusculum, B. nicaeum, and B. lariversi. Moreover, there seems to be a very fast increment with the progressive increase of about up to 30% of characters based on measurements and relationships, and their combinations, which in general refers to approximately three characters and all their combinations. After that, the increase through the node supports seems to happen very slowly (Fig. 8, plot).

A representation of differences and similarities between the topologies found under different approaches adopted here is summarised in a nonmetric multidimensional scaling (NMDS) biplot (Fig. 9). Based on the calculation of RF distances among these topologies found in different reconstruction methods, after rotation procedures during a preliminary principal component analysis (PCA), we plotted into the relationships among these topologies taking into account their ordering onto such an orthogonal space forced to be presented in only two dimensions. This procedure preserved the estimated RF distances among topologies, and the more different topologies found under different analyses, the higher the distance among them on the NMDS biplot. When comparing the topologies found in MP, ML, and BI analyses, there is much more discrepancy between the reconstructions based on MP and ML (RF distances between 16 and 17) than those based on MP and BI (RF distances between 10 and 16) (Table S9). So, the major conflict between the topologies herein observed seems to be related to the adoption of different types of inferences rather than being related to different weighting schemes and/or different combinations of discrete and continuous characters. Without outgroups, the greatest observed conflict occurs between the topology found under ML, with only discrete characters, and those strict consensus trees found under MP, with the matrix of discrete characters and certain combinations of characters based on measurements and relationships (RF distances between 16 and 17). The consistent clades are those comprising the B. pygmeum group (BS = 72%); a monophyletic group composed of B. minusculum + B. plebejum (BS < 50%); the clades including variations of B. lariversi (BS = 53%), B. minusculum (BS = 50%), B. nicaeum (BS = 85%), and B. pygmeum (BS = 56%); and the overarching clade comprising the whole ingroup (BS = 86%) (Fig. S1). For the maximum credibility topologies found under BI, with the matrix of whole dataset, only the clade composed of the variations of B. minusculum + B. plebejum is not consistent throughout topologies found using discrete and continuous characters combined. Accordingly, this clade may have enabled a subtle decrease in the distance between the maximum credibility tree found under BI based only on a matrix of discrete characters and the strict consensus tree found in MP under analysis six. In effect, the RF distance changed from 16 to 14. The other aforementioned consistent clades, however, were found with high levels of Bayesian clade probability (BCP > 85%).

Figure 9. 

Nonmetric multidimensional scaling (NMDS) biplot of a Robinson–Foulds symmetric absolute distance matrix. Ordination of topological differences observed between the topologies obtained under different types of searches and inferences used in this study, after rotation procedures using principal component analysis. The relationships among topologies were forced to be presented into only two dimensions, so that the procedure preserved the distances among topologies after iteratively trying to position the objects in the space in such a way to minimise a stress function (stress = 0.0314). The more different topologies found under different analyses, the higher the distance among the points on NMDS biplot. Abbreviations: A1, analysis 1; A2, analysis 2; A3, analysis 3; A4, analysis 4; A5, analysis 5; A6, analysis 6; SPR0, analysis with a combination of continuous characters producing a Subtree Pruning Regrafting distance equal to zero; SPR2, analysis with a combination of continuous characters producing the Subtree Pruning Regrafting distance equal to two; SPR6, analysis with a combination of continuous characters producing the Subtree Pruning Regrafting distance equal to six; MLDisc, maximum likelihood inference performed on the matrix of only discrete characters; BAYESDisc, Bayesian inference performed on the matrix with only discrete characters; BAYESDiscCont, Bayesian inference performed on the matrix with all datasets combined.

In addition, a monophyletic group constituted by the species of B. plebejum group sensu Lauck, including B. estevezae and B. nessimiani, was recovered in the maximum credibility tree under BI, regardless of the presence of characters based on measurements and relationships (BCP_discrete = 55 %; BCP_{discrete+continuous} = 100%). Using the ACCTRAN criterion when discrete character optimisations were made over BI topology based on the combined discrete and continuous datasets (data not shown), this clade is supported by the anterior portion of the frontogenal suture shorter than its posterior portion (character 5: 2)—a homoplastic synapomorphy—as well as the posterior portion of sclerite 1 of distal median plate articulated with veins 1A and 2A without evident narrowing along its lateral margins (Fig. 10; character 26: 1), another homoplastic synapomorphy. BI also recovered all variations of ‘B. micantulum’ as a monophyletic group either with low support when only discrete characters were used (BCP < 50%) or with slight moderate support when all dataset was used (BCP = 72%).

Figure 10. 

Maximum likelihood and maximum parsimony ancestral states of characters based on male genitalia. Key genitalic characters of the species of the Belostoma plebejum group, particularly those whose aspect of structures was studied rather than their presence or absence and their aspect simultaneously (i.e., only comparable characters), were optimised over the ultrametric maximum credibility tree yielded by Bayesian inference performed on the matrix with all datasets combined. Likelihood of ancestral states are shown as pie charts only in the nodes where evolutionary changes are likely to have happened, and changes in terminal branches were omitted. Characters optimised under maximum parsimony indicated by an asterisk (on the left side) are non-ambiguous. The scoring of the states of each character is indicated by coloured rectangles.

We optimised characters based on male genitalic morphology, particularly those for which we studied the aspect of structures other than their presence or absence and their aspect simultaneously (i.e., only comparable characters), over the maximum credibility tree yielded by BI (Fig. 10). Reconstructed ancestral states of nine characters using the MP procedure revealed only one as non-homoplastic synapomorphy being shared unequivocally by the ingroup: a phallosoma with a strongly curved diverticulum downward, in lateral view (character 31:1). In fact, this synapomorphy is supposed to be diagnostic for the ingroup.

Not only did our ML ancestral state reconstructions determine and reinforce the condition of the aforementioned character, reporting the likelihood of the ancestral conditions at branches where changes were more likely to have happened, but they were also able to report the likelihood of the ancestral condition of other characters based on male genitalia. All species of the B. plebejum group share a diverticulum of phallosoma flat, not longitudinally depressed (character 44: 0). The present MP optimisation gives different evidence for this, having reported otherwise the presence of a medial and longitudinal depression or sulcus weakly developed on diverticulum of phallosoma (character 44: 1), as an ambiguous homoplastic synapomorphy shared by all species of the B. plebejum group. Our MP-reconstructed ancestral state procedures indicate that some characters, which have been considered somehow diagnostic for the B. plebejum group, may have instead occurred more than once within the ingroup, such as the degree of development (character 33: 0) and extension of the arms (character 45: 0). ML ancestral reconstructions over the Bayesian tree otherwise determined arms being shorter than the basal portion of phallotheca (character 45: 2) as a synapomorphy shared by all species of the B. plebejum group.

The narrowing of arms seems to have happened, with equal likelihood, from a more robust condition (character 33: 2) to a narrower one (character 33: 0) (AICc = 38.095; ωAICc = 0.356) with a higher transition rate (2.28 changes) between moderate and strong narrowing change, or as a maximally connected character (i.e., unordered) with equal transition rates (AICc = 38.046; ωAICc = 0.365). A moderate narrowing condition of arms seems to have been the most likely ancestral condition in the clade formed by the species of the B. plebejum s.s., as well as in the clade formed by B. nicaeum + B. pygmeum. Conversely, narrower arms are the most likely condition in the clade formed by B. nessimiani + B. parvum + B. plebejum. MP optimisation provides other evidence, having reported narrower arms as homoplastic synapomorphic condition supporting the clade formed by B. nicaeum and another formed by B. nessimiani + B. parvum + B. plebejum only under ACCTRAN. Unequivocally, the same narrowing condition is an autapomorphy of B. micantulum var2.

The shape of the diverticulum seems to have evolved either from not being truncated at the apex until totally deformed, with a deep V-shaped cut, as well as in the opposite direction, or even as a minimally connected character with the possibility of a transition between 0 and 1 and between 0 and 2 (ORD2: AICc = 35.814; ωAICc = 0.371; SYM: AICc = 36.434; ωAICc = 0.272). Because the ORD2 model provided a slightly better fit, our results suggest that the transition 1→0→1→2 is more likely than 2→1→0, being the transition rate between a non-truncated apex (character 43: 0) and a horizontally truncated apex (character 43: 1) higher than (5.19 changes) from this condition to that one of deep ‘V’-shaped cut deformation (character 43: 2). The most likely ancestral condition of the B. plebejum group s.l., as well as the B. plebejum group s.s., seems to be the presence of a diverticulum with a horizontally truncated apex. MP optimisation unequivocally reported the V-shaped deformed condition of the apex of diverticulum as a non-homoplastic synapomorphy shared by the species of B. pygmeum group. Besides, a reversal seems to have unequivocally occurred in B. minusculum var1, given the untruncated condition of the apex of diverticulum reported to this specimen.

The evolution of the dorsomedial region of the diverticulum seems to have occurred from the absence of a longitudinal groove (flat surface) to the appearance of a very deep groove (a surface strongly depressed) (0→1→2), according to the most likely model (ORD: AICc = 18.580; ωAICc = 0.777). The transition rate between the presence of a poorly developed longitudinal groove and a very deep and evident groove is higher (1.79 changes) than between the non-grooved condition and the slightly grooved diverticulum condition (0.84 changes). The most likely ancestral condition of the ingroup seems to be the absence of a groove in the dorsomedial region of the diverticulum (character 44: 0), a plesiomorphic condition given by MP optimisation. However, the slight development of a groove onto the diverticulum (character 44: 1) seems to be the most likely ancestral condition of the clade formed by the species of B. plebejum s.s. Besides, a very deep groove is the most likely ancestral condition in the clade formed by the B. pygmeum group and in the clade formed by B. estevezae + B. lanemeloi sp. nov. + ‘B. micantulum’ + B. nieseri sp. nov. MP optimisation reported rather similar results: the slight development of a groove onto the diverticulum is a homoplastic synapomorphy shared by the species of B. plebejum s.l. only under ACCTRAN, and the strongly grooved condition of diverticulum (character 44: 2) is a homoplastic synapomorphy (a parallelism) in the clade formed by the B. pygmeum group, as well as in the clade formed by B. estevezae + B. lanemeloi sp. nov. + ‘B. micantulum’ + B. nieseri sp. nov., also only under ACCTRAN.

Considering their length when compared with the longitudinal length of the anterior portion of the phallosoma, the arms seem to have evolved with equal transition rates (0.68 changes) between maximally connected available conditions. Accordingly, arms shorter than the base of the phallosoma (character 45: 2) seem to be the most likely ancestral condition of the clade formed by the ingroup (ER: AICc = 31.044; ωAICc = 0.498). This condition becomes increasingly most likely as the clades become less internal and slightly more likely in the clade formed by the species of B. pygmeum group. MP optimisation unequivocally reported that the length of arms equals to that of the base of the phallosoma (character 45: 0) is a homoplastic synapomorphy (parallelism) that defines both the clade formed by B. nicaeum + B. pygmeum and the clade formed by B. minusculum. Moreover, the condition is unequivocally an autapomorphy of B. plebejum. The condition of the arms being shorter than the base of phallosoma (character 45: 2) is a homoplastic synapomorphy (reversal) in the clade formed by B. lariversi, in the clade formed by B. nessimiani + B. parvum, and in the clade formed by B. estevezae + B. lanemeloi sp. nov. + ‘B. micantulum’ + B. nieseri sp. nov. only under ACCTRAN.

The unique binary character proposed as diagnostic that appeared independently in some instances in the groups of species studied is the absence of a protuberance on the apicoventral portion of the diverticulum (character 46: 1) because it occurs both in B. pygmeum and in the clade formed by B. oxyurum. Besides, the aspect of the apicoventral portion of diverticulum appears to have evolved with equal rates of transition from either 0 to 1 or vice versa, with a rate of 0.59 changes. Despite being more likely, this model (ER: AICc = 17.139; ωAICc = 0.425) is only a little more likely than the others evaluated (ORD₀₁: AICc = 17.693; ARD: AICc = 18.178). According to the ER, the protruding aspect of the apicoventral portion is the most likely condition of the clade formed by the ingroup, which only changes in B. pygmeum.

This is the first phylogenetic study to infer the relationship among members of a group of small species in Belostoma. The present analyses based on morphological characters (seven continuous and 42 discrete) (see supplementary file 3: List of characters) seem to provide support for the B. plebejum group as monophyletic, as proposed by Nieser (1975), now with two new species included. The results, therefore, corroborate the ideas of monophyly previously suggested by Estévez and Polhemus (2007) and Ribeiro and Alecrim (2008). Accordingly, in the strict consensus tree of analysis six (Fig. 5) with combined datasets and Mirande’s (2009) K estimated from data, the following unambiguous non-homoplastic synapomorphies provide evidence for the monophyly of the ingroup: (1) prosternal carina usually not elevated and much broader at its base than the apex, in lateral view (character 9: 1); (2) posterior pronotal width (usually not the greatest one) twice its median length (character 20: 2); (3) pilosity covering about one third of the connexivum (character 29: 2); and (4) phallosoma with the diverticulum strongly curved downward, in lateral view (character 31: 1).

The enlargement of the base of the prosternal carina is probably independent of its elongation, which seems to be evident through the analyses. Previous clues provided by some essays about the development of prosternal carina and body size in Belostomatidae (Ribeiro et al. 2014, 2017, 2018) also support the idea that these conditions have to be considered separate components. Estévez and Polhemus (2007) established the condition of the prosternal carina being both wide and projected to the anterior part of body as diagnostic for B. nicaeum and B. plebejum. One can still highlight that the poorly elevated condition of the prosternal carina, in lateral view, is shared by the members of B. plebejum group s.l. and has been thought to be characteristic of its members (Nieser 1975; Nieser and Melo 1997; Ribeiro 2007; Ribeiro and Alecrim 2008; Stefanello 2021). Besides, based on the sketches prepared by Nieser (1975), the prosternal carina of the studied specimens of B. amazonum—a species we used as an outgroup—apparently shows most of the plesiomorphic conditions highlighted by our results. In their taxonomic revision of the B. denticolle group, Estévez and Polhemus (2001) suggested that the condition of the prosternal carina projected anterad is diagnostic for B. amazonum and B. denticolle. However, it is in fact homoplastic because that condition is also shared with other species of the B. plebejum group. Likewise, the presence of tubercles on the apex of the prosternal carina with a broader base (character 15: 1) is diagnostic for B. oxyurum, and this has easily allowed its recognition. However, some exemplars of B. lariversi also show this condition on the prosternal carina, and our results provide further evidence that the aforementioned condition is also homoplastic (parallelism), favouring the hypothesis of independent origins of this condition in B. pygmeum group and in some outgroups.

Conversely, our results based on MP inference do not corroborate the restriction of the B. plebejum group to two other more inclusive groups, but rather suggest the B. plebejum group s.s. to be paraphyletic with respect to the B. pygmeum group. The existence of an apicoventral protuberance on the diverticulum of phallus (character 46: 0) was considered by Lauck (1962) as a diagnostic feature of B. plebejum group. One could have plausibly accepted that the species that the author used to create that group of species had indeed borne such a characteristic. However, the aforementioned diagnostic condition is not exclusive (Ribeiro and Alecrim 2008). Besides, the absence of this protuberance (character 46: 1) is unequivocally a homoplastic synapomorphy shared by the B. pygmeum clade, whereas the former condition proposed by Lauck (1962) is plesiomorphic in the present study. Our results still provide evidence that the apomorphic condition arose independently in both the B. pygmeum and B. oxyurum clades (parallelism)

Bayesian inference of combined datasets instead indicates B. plebejum group s.s. is a monophyletic group (BCP_{combined datasets} = 100%), which justifies a future study of this group, comprising not only morphological but also molecular datasets, in order to verify its monophyly as proposed by Lauck (1962). Lastly, the putative synapomorphies that support the monophyly of that clade, reported by MP optimisation over the maximum credibility tree, are ambiguous and have not involved any diagnostic characteristic regarding male genitalic morphology (Fig. 10). Accordingly, an anterior portion of the frontogenal suture that is shorter than the posterior portion (character 5: 2) seems to provide some phylogenetic signal to such a group, which had been highlighted by Estévez and Polhemus (2007) and Ribeiro and Alecrim (2008). However, we have suggested for the first time that the posterior portion of sclerite 1 of the distal median plate articulated with veins 1A and 2A without evident narrowing along its lateral margins (character 26: 1) is a putative synapomorphy. ML and BI also recovered the clade formed by B. estevezae + B. lanemeloi sp. nov. as a monophyletic group, despite having low support (Fig. S1). Accordingly, this clade is unequivocally supported by the lateral margin of the diverticulum of the phallossoma that is strongly curved mesally, in dorsal view (character 47: 1), a non-homoplastic synapomorphy. Such a condition is instead a homoplastic autapomorphy (parallelism) of the aforementioned species under MP inference. We have indicated this feature for the first time as important in the taxonomy of these species, and this may be an outcome of the use of more adequate treatment of measures and relations during the analyses, as well as the achievement of fewer errors. Indeed, scenarios with high rates of homoplastic or superimposed changes make the MP inconsistent. It is plausible to suggest that the discrepancy observed between the existence or absence of these topologies, which seem to be dependent on the method used, could be due to the presence of different evolutionary rates through the constructed characters, because we clearly obtained the phylogenetic signal of this clade and the increased accuracy of the resulting topologies. Therefore, that scenario better handles parametric models (Wright 2019), especially when matrices with a low number of characters are implemented (Wright and Hillis 2014).

Our results also provide evidence to support another hypothesis of monophyly: the clade formed by all species of the B. pygmeum group sensu Lauck (1962), considered by some authors to be a group of specialised species (Nieser 1975; Estévez and Polhemus 2007; Ribeiro and Alecrim 2008). Four synapomorphies support the monophyly of the clade, two of which are unequivocally non-homoplastic and based on characters of male genitalia: (1) margins of the diverticulum of the phallosoma clearly thickened, in lateral view (character 32: 1), and (2) an acute apex of the diverticulum, in lateral view (character 48:1) (Fig. 10). Furthermore, the wider condition of the diverticulum, evidenced by a width that is much greater than its length, is a continuous character that helps to support the monophyly of this group. Nieser and Melo (1997) used this relation for the first time in their book about the taxonomy of Belostoma species, allowing the differentiation of some species from other groups of species, such as B. anurum (Herrich-Schäffer, 1848) and B. dallasi De Carlo, 1932. This relation has been used efficiently by some other authors (e.g., Ribeiro 2000, 2007; Ribeiro et al. 2017) and has allowed a more precise understanding of the taxonomy of other species, given its common use throughout some identification keys (e.g., Ribeiro and Alecrim 2008).

The thickened condition of the laterodorsal margins of the diverticulum also allows easy recognition of the representatives of the B. pygmeum group, despite not being clearly visible in B. pygmeum. As mentioned before, Nieser (1975) established such a peculiar condition of the diverticulum of the phallosoma as an extreme case of specialisation occurring throughout the species of the B. plebejum group. The V-shaped deformed condition of the apex of the diverticulum (character 43: 2) is also something quite conspicuous in B. lariversi, and it was recovered as an unambiguous homoplastic synapomorphy for this group in the present study either under MP or BI. Furthermore, MP inference under the implicit weighting scheme recovered such a clade with low support, but model-based inferences corroborate the monophyly of that clade with high posterior probability (BCP_discrete = 98%; BCP_{discrete + continuous} = 100%) and with slight moderate bootstrap support under ML inference (BS = 72%) (Fig. S1). In addition to the already mentioned putative synapomorphies that support that clade, optimised over the maximum credibility tree, the presence of a longitudinal sulcus strongly developed in the diverticulum of the phallosoma, much evident in dorsal view (character 44: 2), supports the clade only under ACCTRAN. However, our results suggest the hypothesis of independent origins (parallelism) of that condition in this clade and in the clade formed by B. estevezae, B. lanemeloi sp. nov., ‘B. micantulum’, and B. nieseri sp. nov. (BPS_discrete = 52%; BPS_{discrete + continuous} = 86%) also only under ACCTRAN. Curiously, many authors have considered the aforementioned condition of the mesal portion of the diverticulum to be diagnostic in B. micantulum (e.g., Nieser and Melo 1997; Estévez and Polhemus 2007; Ribeiro 2007) and, thus, such a characteristic of this species must be re-evaluated in the future.

In the present study, ML and BI assigned B. lariversi as a sister-group of B. nicaeum + B. pygmeum (BPS_discrete = 65%; BPS_{discrete + continuous} = 98%). However, analysis six of MP inference supported a clade composed of (B. nicaeum, (B. lariversi + B. pygmeum)) in parenthetical form (Fig. S1). According to the MP optimisation of this character over the maximum credibility tree, the presence of arms shorter than the length of the base of the phallosoma (character 45: 2) is an ambiguous synapomorphy of the clade formed by B. estevezae, B. lanemeloi sp. nov., ‘B. micantulum’, and B. nieseri sp. nov. under the ACCTRAN criterion, whereas the arms becoming as long as the base of the phallosoma (character 45: 0) unambiguously supports the clade formed by B. nicaeum + B. pygmeum and also appears independently in B. plebejum and B. minusculum (parallelism). Conversely, arms shorter than the length of the phallobasis is plesiomorphic in the strict consensus tree of analysis six recovered by MP inference, whereas arms becoming as long as the base of the phallosoma unambiguously supports the clade formed by the species of B. pygmeum group + B. plebejum + B. minusculum, with a reversal to state 2 in B. lariversi (Fig. 5). The different dispositions of arms, as well as their degrees of development, have long been used to recognise many Belostoma species, such as representatives of B. denticolle, B. oxyurum, and the B. plebejum group, over the history of the taxonomy of small species of Belostoma (Lauck 1962; Ribeiro and Estévez 2007; Estévez and Polhemus 2007; Ribeiro and Alecrim 2008). Moreover, Ribeiro and Estévez (2009) established the condition of smooth convergence towards the apex of them, taking into account their position over the diverticulum (character 42: 2), as a diagnostic characteristic of members of B. sanctulum (an outgroup), one species of B. oxyurum group. Likewise, within the ingroup, arms of the phallotheca diverging towards their apexes (character 42: 3) is a plesiomorphic condition in the strict consensus tree recovered by analysis six of MP inference; it is also shared by the clade formed by the species of the B. pygmeum group. Accordingly, arms converge abruptly at their apex (character 42: 1) in B. micantulum var1, converge smoothly at their apex in B. nessimiani, and are parallelly disposed (character 42: 0) in B. plebejum.

Based on the previous paragraphs, some increment in morphological complexity—especially through male genitalia—seems to occur within the B. plebejum group, and they could somewhat be associated with a process of body size decreasing in Belostoma. Although B. parvum is the smallest species of genus, this process seems to culminate in the species of B. pygmeum group. The most likely process of evolution in body dimensions, according to the models evaluated in this study, is the one supported by a small force of selection channelling the characteristics based on measures and relation to an optimal global value under constant evolution rate (OU model: ss = -378.03; ps = -377.49; ln Bayes factor = 3.3 - Yule + homogeneous OU model vs. Yule + relaxed BM model, or ln Bayes factor = 5.2 - Yule + homogeneous OU model vs. Yule + homogeneous BM model) (Hansen and Martins 1996) (Table S10). Many aforementioned novelties seem to support the hypothesis—especially when male genitalic characteristics are optimised over the maximum credibility tree, obtained under BI with combined datasets—because they showed a certain tendency to arise in groups of species with smaller body sizes. Shorter arms arising of phallotheca (or phallobasis) and not reaching the lateral margins of diverticulum, for example, has been suggested as diagnostic for the B. pygmeum group (Lauck 1962). Besides, Nieser (1975) established that the plesiomorphic condition of this character (i.e., arms reaching the lateral margins of the diverticulum) would also represent a case of specialisation found in the B. plebejum group, consequently allowing the emergence of the lineage that originated the species of the B. pygmeum group. As suggested by our results, however, the clade formed by B. lanemeloi sp. nov., B. nessimiani, and B. nieseri sp. nov. under MP inference is that one that was supported unambiguously by the non-homoplastic synapomorphic condition of this character. Unfortunately, this clade was not recovered under BI, which may be a result of the presence of many homoplasies underlying the evolution of those species, the mutual presence of characters with slower or faster evolutionary rates (Wright and Hillis 2014), somehow highlighted by our investigation about evolutionary rates of some male genitalic characters, or simply an increase in accuracy that occurs as long as the size of character matrices also increases (Parins-Fukuchi 2018).

We tested the phylogenetic integrity (de Pinna 1999) of four species of the B. plebejum group. Belostoma lariversi (BS = 68%; BRel = 45%), B. minusculum (BS < 50%; BRel = 34%), B. nicaeum (BS = 80%; BRel = 44%), and B. pygmeum (BS = 74%; BRel = 42%) with their exemplars formed their respective monophyletic clades, except for ‘B. micantulum’, in all MP analyses, when combining or not combining continuous with discrete characters (Figs 5 and 8).

The variations of B. lariversi and B. nicaeum formed their respective monophyletic groups, suggesting phylogenetic integrity of these species either with MP or BI. Four unambiguous homoplastic synapomorphies supported the clade formed by the variations of B. lariversi, but only one referred to the degree of development of the arms of the phallosoma (Fig. 5, character 45: 2), a unique recovered character based on the male genitalia, while three other unambiguous homoplastic synapomorphies supported B. nicaeum. In this case, the key characteristics based on male genitalia are as follows: narrower arms of the phallotheca (character 33: 0) and the diverticulum of the phallosoma with a weak and poorly conspicuous dilated portion along its apex, in dorsal view (character 49: 0). Besides, both species were supported by a conspicuous enlargement of diverticulum in ventral view compared with its length (character 30), especially in B. lariversi (Fig. 8). Therefore, despite the variations of the prosternal carina pattern morphology, it is not possible to refute the hypothesis that the terminals assigned to B. lariversi and B. nicaeum form respective single species. Such a variation in the prosternal carina has been documented by some papers (Ribeiro 2007; Ribeiro et al. 2017), which corroborates the idea that there is much intraspecific variation in this structure.

The variations of B. minusculum also formed a monophyletic group under the MP criterion (and also under BI), suggesting the phylogenetic integrity of the species despite the lowest support recovered—probably due to the presence of a broad prosternal carina with its anterior margin straight (character 16: 0)—the unique unambiguous homoplastic synapomorphy supporting that clade. The variations are among the ‘Druckknopf’ system and male genitalia, especially in the disposition of arms (characters 17, 27, 42, and 43), suggesting some degree of polymorphism among these features. Against our expectation, the complementary structure of the apparatus, located in the thorax, also seems to be polymorphic, although Ribeiro and Alecrim (2008) and A. L. Estévez (pers. comm.) indicated such a characteristic is important to identify species of B. plebejum group.

Under the MP criterion (and also under the Bayesian criterion), B. pygmeum is supported by the absence of a median apicoventral protuberance on the diverticulum, a characteristic based on the morphology of male genitalia (Fig. 10; Fig. S1), which Lauck (1962) considered being diagnostic to the species of B. plebejum group. Again, the complementary structure of ‘Druckknopf’ system apparatus, as well as the shape of the prosternal carina, are herein indicated as polymorphisms.

On the other hand, the two variations of ‘B. micantulum’ did not form a monophyletic group, suggesting no phylogenetic integrity of this species with the use of the MP criterion (Fig. 5). The exemplars were instead paraphyletic in all analyses, except for analysis one. Variation 1 of B. micantulum, occurring in the states of Pará, Amazonas, Paraíba, Mato Grosso, Espírito Santo, and Rio de Janeiro, formed a monophyletic group with all species of the B. pygmeum group + B. micantulum var2 + B. estevezae + B. plebejum + B. minusculum, despite having low support (BS < 50%; BRel = 33%). The group was supported by the presence of a rectangular prosternal carina, obtuse at the apex, an unambiguous non-homoplastic synapomorphy (character 18: 0). Although all combinations of continuous characters corroborate the clade (Fig. 8), variation 2 of B. micantulum constituted a more internal clade, a finding that indicates segregation of characters and differentiation (de Pinna 1999). Such an internal clade (BS < 50%; BRel = 15%) was not supported by some combinations of continuous characters together with discrete ones, despite being supported by the presence of a wide prosternal carina with anterior margin curved (character 16: 1), an unambiguous homoplastic synapomorphy. Although the variations co-occur partially in Brazil, because both are found in Pará, Amazonas, and Espírito Santo States, some differences in arms of the phallotheca seem not to be just intraspecific variations as long as arms are abruptly convergent at their apexes in representatives of variation 1, which authors have reported to be more widely distributed. However, in representatives of variation 2 arms are evidently narrower, which only occurs in Pará, Amazonas, and Espírito Santo States. Lastly, variation 1 of B. micantulum is slightly narrower than variation 2 (Fig. 8; Table S3).

Under BI, variations of ‘B. micantulum’ formed a monophyletic group with moderate support (BCP = 84%) (Fig. S1). As mentioned before, analysis one of the MP criterion also recovered these variations as a monophyletic group. It might be that the difference between topologies found under the MP criterion and model-based methods is an outcome of the interaction of discrete and continuous characters, as long as most likely evolutionary models were fitted to different evolutionary rates likely to exist among discrete and continuous datasets (Felsenstein 1985; Parins-Fukuchi 2018). Parametric methods, such as BI, are also thought to produce better results than MP in situations with small datasets and the presence of missing data (Wright and Hillis 2014; O’Reilly et al. 2018; Puttick et al. 2019). Moreover, BI usually generates results that are slightly more similar to MP with equal weighing schemes (Silveira et al. 2021). Thus, taking into account these discrepant results in the present study, we are much inclined to not deliberately refute the hypothesis that the variations assigned to ‘B. micantulum’ form a single species. Nevertheless, historically the specific delimitation of B. micantulum has been a problem. Much has been said about the existence of putative species complexes of small Belostoma species. Accordingly, some important small cryptic species, in terms of general body appearance, are as follows: B. micantulum, B. plebejum, B. pygmeum (by some authors), and B. denticolle (A. Bachmann pers. comm.; Nieser 1975; Nieser and Melo 1997; Ribeiro 2007), some of them being differentiated from one another only by body size (Nieser and Melo 1997). Furthermore, the use of a diversified diet to monitor the life cycle of B. micantulum, under laboratory conditions, has suggested the existence of remarkable phenotypic plasticity regarding its size (Pereira et al. 1991, see Pereira 1992 for B. plebejum). Therefore, the idea of these variations found in B. micantulum being instead an outcome from differentiated interactions between specimens and environment through their lives (Armbruster et al. 2004 with examples in plants; Klingenberg 2015) is not unplausible.

Stefanello et al. (2018) highlighted that taxonomy strongly or exclusively based on a study of male genitalic morphology in Belostoma may not always be recommended because some intraspecific variation in male genitalia has been found inside some Belostoma species with the use of geometric morphometric tools. The monophyly of ‘B. micantulum’ we have obtained with the use of model-based methods should also be interpreted with caution. It might still be coherent because the combination of what is supposed to be either uninformative or informative in MP provides only an overall idea of a plethora of different evolutionary rates likely to underlie their respective datasets (Wright and Hillis 2014). Besides, combined with such monophyly recovered, the existence of differences found among the disposition of arms of the phallotheca described in the previous paragraph reinforces the suggestion of the existence of polymorphism in some Belostoma species, including ‘B. micantulum’. Recent results with the use of molecular dataset recovered some Belostoma lineages inside two or three small species with some level of polymorphism or phenotypic plasticity in male genitalia (Gauterio et al. 2021). This would also explain the moderate level of homoplasies found in those characters based on male genitalia that we have studied in more detail.

Contrarily to Thiele’s (1993) approach, we treated continuous characters as continuous series that we did not divide into an arbitrary set of states (Goloboff et al. 2006). In addition, we have used different implicit character weighting schemes and different amounts and combinations of continuous characters to investigate the effect of these practices on the character matrix, considering the evenness in the production of repeated relationships and character congruence. So, we have developed a new procedure for combining continuous characters with a discrete dataset, under MP specifically for this study. Besides, we compared recovered topologies of these analyses under the MP criterion in the light of those obtained with explicitly statistical methods (i.e., ML and BI). In the present study, comparing the topologies obtained by different approaches clearly showed the presence of different scenarios in terms of heterogeneity of evolutionary rates among characters, because we noted a certain disparity between the results obtained with different reconstructions in the biplot produced from the NMDS approach. Based on a triangular matrix of RF distances (Robinson and Foulds 1981) estimated among topologies (Table S9), the disposition of points on the biplot (Fig. 9) indicated that the discrepancy between the topologies reconstructed both under explicitly statistical methods and MP, in general, seems to be much more important than the supposed discrepancy related to the use of different weighing schemes and/or different combinations of discrete and continuous characters. However, it is worth mentioning that the disproportionate number of discrete characters compared with the continuous dataset could be responsible for the almost nonexistence of differences among our topologies recovered with different combinations of measures and relations.

Issues related to the controversial use and treatment of continuous characters have received more attention in parsimony analyses (see Hauser and Presch 1991; Slowinski 1993; Kitching et al. 1998; Wiens 2001; Goloboff et al. 2006; Baur and Leuenberger 2011; Kock et al. 2014). In analyses carried out under different weighing schemes, the influence in using and weighting characters based on the implicit weighting approach (Goloboff 1993) onto strict consensus topologies was similar, whether or not we used continuous characters. However, the influence of adopting any weighting was greater in those topologies where weight equal to one was assigned to all characters than in those where characters received different implicit weightings, according to the number of homoplasies that a referred character could present and topological stability (as explored in Barão et al. 2020). Only when there was a weight equal to one assigned to all discrete characters (analyses one and four) did we find a greater number of alternative hypotheses, which resulted in a greater number of taxa integrating polytomies in strict consensus trees. Nevertheless, some of the clades are still in agreement with the clades recovered in analyses two, three, five, and six, whose matrices were submitted to implicit weighting procedures. Similar clades recovered through different weighting and ordering schemes constitute consistent and stable hypotheses of phylogenetic relationships (Gutberlet 1998; Barão et al. 2020). The B. plebejum group sensu Nieser (1975) and the phylogenetic integrity of some species are clear examples of this condition.

An overlooked fact instead is that a multistate character (which could be ordered or not ordered) exerts greater influence in discriminating among alternative tree topologies than a binary character (Donoghue and Ree 2000; Wiens 2001). A character with four states will have three times as much influence as a binary character, while one with 10 states will have three times as much influence as one with four states and nine times as much influence as a binary character, which is particularly relevant for continuous characters broken into an arbitrary number of discrete states (Swofford and Begle 1993). Furthermore, discrete binary characters have smaller state space and develop much more slowly at low evolutionary rates than multistate characters (Donoghue and Ree 2000; Wright and Hillis 2014; Parins-Fukuchi 2018). Therefore, treating measures and relations as such is also an advantage in scenarios with high evolutionary rates because these characters are less subject to the saturation phenomenon, as they present a high magnitude of variation (Goloboff et al. 2006) and, in turn, the loss of phylogenetic signal is lower (Wright and Hillis 2014).

Continuous characters, therefore, tend to retain more information regarding transformation and segregation, especially when submitted to explicitly statistical methods of tree reconstruction, such as BI (Wright and Hillis 2014; Parins-Fukuchi 2018; Puttick et al. 2019). This seems to be the case in the present study because there was a slight increase in the stability of some nodes, despite a slight decrease in retention index values through MP analyses with the use of different combinations of continuous characters.

Nevertheless, there are known negative effects of high magnitude of variation imposed by characters based on measures and relations under implicit weighting procedures during MP reconstructions (‘an asymmetric influence’ according to Barão et al. 2020). The effects are usually reduced when those characters are standardised or normalised so that they can be weighted more accurately according to their degree of homoplasy (Koch et al. 2014). In the present study, we adopted weighing schemes associated with scaling or standardisation of metric and meristic data during MP procedures. However, due to the large amount of homoplasy normally attributed to the continuous characters (Wiens 1995) and a reduction likely due to the dominance of characters with a large range of values (Goloboff et al. 2006; Baur and Leuenberger 2011; Koch et al. 2014), the phylogenetic signal of discrete character partition still produced some discrepancies. These issues become evident with the extent of conflict or concordance found in those reconstructions produced by probabilistic methods. Hence, representatives of the B. plebejum group underwent a scenario with heterogeneity among the evolutionary rates of their studied morphological traits, more specifically when we estimated the evolutionary rates of some key characters based on male genitalia.

The procedure we have adopted involves the gradual combination of continuous characters to the matrix of discrete characters during the MP procedure. It showed that the progressive increase in the number of continuous characters necessarily increased SPR distances between the consensus tree produced only by the partition of discrete characters with implicit weighting (analysis three) and consensus trees recovered with different combinations of the two partitions, although accuracy and stability were not dependent on the SPR distance values. However, there was a slight overall progressive increase in stability with the inclusion of up to approximately 30% of traits based on measures and relations, as well as their combinations (Fig. 8). Some authors have mentioned that errors during phylogenetic reconstructions are generally smaller when continuous characters are used in conjunction with discrete character matrices (Mounce et al. 2016; Parins-Fukuchi 2018), a strategy that agrees with a ‘total evidence’ reasoning (de Queiroz and Donoghue 1995). Continuous characters are assigned more objectively, while they bring much accuracy because they considerably increase the size of character matrices (Gatesy et al. 1999). On the other hand, accuracy might be reduced if such included characters lead to greater covariance among them, even though they behave relatively well in terms of consistency because of lower levels of saturation so that their phylogenetic signal is not compromised (Donoghue and Ree 2000; de Bivort et al. 2010; Escapa and Catalano 2013; Wright and Hillis 2014). This phenomenon would cause a decrease in their informational content, probably caused by their residual covariance, even under an explicitly statistical model of continuous data evolution (Parins-Fukuchi 2018). Another possibility of incongruence might be related to the existence of different selective pressures upon the modules being represented by these partitions (Clarke and Middleton 2008). As already mentioned, under MP the procedures adopted in the present study showed that an increase in node supports is not so evident after the inclusion of about 30% of continuous characters. We suggest that this seems to be caused by the existence of some covariance component that was gradually added into datasets.

Using model-based methods, our recovered phylogenetic reconstructions indicated that accuracies provided by the two dataset partitions appear to be the same, even though they behave slightly differently when used in different reconstruction methods. IQ-TREE, for example, seems to perform worse when compared with other ML software, especially when character matrices have missing data and are small (Vernygora et al. 2020). In addition, topologies recovered with BI are generally more accurate than those recovered with ML inference when both are evaluated by estimating RF distances (Brown et al. 2017). Our continuous dataset is probably under selective constraint force that seems to act similarly upon each of the seven characters and in a constant way over time. The use of a homogeneous OU model (Hansen and Martins 1996) as the most likely process of evolution to our continuous dataset highlighted such a constraint, therefore, improving the performance of our inferences. In turn, the accuracy of our reconstructions was also improved, regardless of the covariance likely to be found among these characters. This indicates that the evolutionary process of those continuous traits slowly reached stationarity towards an optimum, due to the low estimated overall evolutionary rate (α = 0.202) (Parins-Fukuchi 2018; Harmon 2019).

Polytomies may consistently show low RF distances with any other topology in a given tree space (Rannala et al. 1998). Accordingly, this could also explain why the consensus trees recovered in analyses one and four of MP inference were closer to other topologies recovered by probabilistic methods in the NMDS biplot. Metrics based on distances focussed on topologies, such as RF distance, usually evaluate differences and similarities among trees well (Kuhner and Yamato 2015). However, such a metric has been known to benefit methods of tree reconstruction—and datasets—that recover consensus trees with polytomies (Vernygora et al. 2020), therefore estimating low distances among them when compared with fully resolved trees. Likewise, trees without polytomies but with only a few very stable nodes—that is, with high support values—tend to be a problem. Unlike trees without polytomies but with few stable nodes, polytomies can still be treated as partially correct and, unfortunately, there is currently no metric that avoids such a bias created by its lacking of resolution (Rannala et al. 1998; Holder et al. 2008; O’Reilly and Donoghue 2018). In the present study, this reasoning allowed us to adopt RF distance—a metric that strongly penalises fully resolved trees with low-supported nodes—thus reducing high estimates of RF distance among two trees when they are not quite different from each other (i.e., type-I error) (Vernygora et al. 2020). Regarding MP inferences, the NMDS biplot shows that the additional use of continuous dataset did not improve the phylogenetic information supposed to be present in our discrete dataset, probably because of the disproportional number of discrete characters compared with continuous characters in our matrix.