Relationship between genetic diversity and morpho-functional characteristics of flight-related traits in Triatoma garciabesi (Hemiptera: Reduviidae)

Background Triatoma garciabesi, a potential vector of the parasitic protozoan Trypanosoma cruzi, which is the causative agent of Chagas disease, is common in peridomestic and wild environments and found throughout northwestern and central Argentina, western Paraguay and the Bolivian Chaco. Genetic differentiation of a species across its range can help to understand dispersal patterns and connectivity between habitats. Dispersal by flight is considered to be the main active dispersal strategy used by triatomines. In particular, the morphological structure of the hemelytra is associated with their function. The aim of this study was to understand how genetic diversity is structured, how morphological variation of dispersal-related traits varies with genetic diversity and how the morphological characteristics of dispersal-related traits may explain the current distribution of genetic lineages in this species. Methods Males from 24 populations of T. garciabesi across its distribution range were examined. The cytochrome c oxidase I gene (coI) was used for genetic diversity analyses. A geometric morphometric method based on landmarks was used for morpho-functional analysis of the hemelytra. Centroid size (CS) and shape of the forewing, and contour of both parts of the forewing, the head and the pronotum were characterised. Length and area of the forewing were measured to estimate the aspect ratio. Results The morphometric and phylogenetic analysis identified two distinct lineages, namely the Eastern and Western lineages, which coincide with different ecological regions. The Eastern lineage is found exclusively in the eastern region of Argentina (Chaco and Formosa provinces), whereas the Western lineage is prevalent in the rest of the geographical range of the species. CS, shape and aspect ratio of the hemelytra differed between lineages. The stiff portion of the forewing was more developed in the Eastern lineage. The shape of both portions of the hemelytra were significantly different between lineages, and the shape of the head and pronotum differed between lineages. Conclusions The results provide preliminary insights into the evolution and diversification of T. garciabesi. Variation in the forewing, pronotum and head is congruent with genetic divergence. Consistent with genetic divergence, morphometry variation was clustered according to lineages, with congruent variation in the size and shape of the forewing, pronotum and head. Graphical Abstract Supplementary Information The online version contains supplementary material available at 10.1186/s13071-024-06211-x.

Triatoma garciabesi has historically been considered to be a single species or synonymised with T. sordida due to morphological similarities and partial overlap in geographical distribution.Its validation as a separate species was published in 1967 and was based on the observation of morphological differences in the head and phallic structures [10].However, T. garciabesi was synonymised with T. sordida by Lent and Wygodzinky in 1979 [3], with these authors concluding that the phenotypic differences were not sufficient to consider two distinct species.Regardless, two decades later, evidence based on microhabitat characters, male genitalia and cytogenetics reaffirmed its taxonomic status as T. garciabesi [for review see 11].Various approaches have been used to justify the taxonomic validity of this species, including geographical distribution [7,10,12], morphology of the head, eggs and genitalia [13,14], morphometry of the forewings, head and pronotum [8,14,15], cuticular hydrocarbons [9], isozymes [17], cytogenetic characteristics [7] and molecular analyses [7,8,18].Triatoma garciabesi is usually found in peridomestic environments, mainly associated with trees where chickens roost [19].The species occasionally invades rural houses [20,21] and has been recorded in wild environments [22,23].Geographical distribution records based on morphological characteristics indicate that this species is present in the semi-arid regions of northwestern and central Argentina and western Paraguay [7,15,16,24].Based on cytogenetic analysis, its distribution has been extended to the Bolivian Chaco (Santa Cruz department), but there are no records of its presence in other areas of Bolivia [7].According to an ecological niche study, the drier areas of northwestern Argentina provide ideal conditions for the occurrence of T. garciabesi [15].
The analysis of nuclear and mitochondrial DNA (mtDNA) sequences has proven to be suitable for determining genetic diversity and population structure of triatomines [25][26][27].In other insects, the mitochondrial cytochrome c oxidase I gene (coI) has been widely used to improved understanding of species diversity and delimitation [28].Assessing the genetic differentiation of a species across its geographical range can provide insights into dispersal patterns and habitat connectivity.Dispersal also plays a significant role in shaping the spatiotemporal distribution of a species' genetic diversity [29].Active dispersal tendencies of a species or its populations evolve according to changes in locomotory and navigation traits [30].The relationship between the shape and function of insect forewings would be subject to strong selective pressures, with not only behavioural activities like oviposition and foraging preferences being taken into consideration, but also flight-related characteristics that enhance mobility within a specific environmental context [30,31].In triatomines, flight is considered to be the most important active dispersal mechanism [32,33].For Triatoma infestans Klug, 1834, flight distances of between 200 m [34] and 1500 m [35] have been reported.Hemipteran forewings, named hemelytra, are composed of a stiff proximal region and a more flexible membranous apex.The morphological structure of the hemelytra is linked to their function: the membranous part can be deformed by aerodynamic and inertial forces, and the stiff part limits deformations and provides support [36].The shape, size and development of the parts of the insect forewing define the characteristics related to flight efficiency.For example, in butterflies and damsflies, elongated and narrow wings determine greater deformation capacity with lower energy cost for long distance flights, whereas short and broad forewings are compatible with short, slow and less stable flight [37,38].Although wing morphology in triatomines is quite different, for T. guasayana, elongated and thin wings were suggested to allow sustained long-duration flights, and short and wide wings would be related to short flights [39].Wing aspect (4 × [wing length 2 /wing area]) is also an indicator of flight capacity, with individuals having high values of this parameter being able to fly faster and farther than individuals with low wing aspect values [40].
In insects, the ability to orientate in space and in relation to the specific environmental cues (seasonal events, daily cycles and specific temporal cues) contributes to the navigational abilities, which are linked to sensory and spatial cognitive abilities.All of these abilities are provided by organs, such as eyes for vision, antennae for olfaction and the brain for celestial navigation and cognitive maps (e.g.[30]).In several insects, changes in head phenotype have been associated with flight capacity [30,41,42].In the triatomines Mepraia spinolai Porter, 1934 and T. guasayana, some characteristics, such as narrow head with highly developed compound eyes, are considered to be associated with flight dispersal [41][42][43].When developed, flight muscles are located under the pronotum.In some triatomines, an increased dispersal capacity was found to coincide with a wide thorax, despite the intermediate size of forewings.A wide thorax may imply highly developed forewing muscles [44].
Flight dispersal is a complex and plastic attribute determined by a suite of morphological, physiological and behavioural characteristics [29], all of which have genetic basis and are influenced by the environment [45].Similar to several other insects, triatomines respond rapidly to environmental changes; therefore, they are considered a suitable model for studying the effects of genetic diversity on changes associated with dispersal.For example, evolutionary changes in flight morphology have been linked to range expansion, colonisation success, habitat fragmentation and migration.The aim of the present study was to understand the structuring of genetic diversity of T. garciabesi, the variation in morphology of dispersal-related traits in relation to genetic diversity and the relationship between morphological characteristics of dispersal-related traits and the current distribution of genetic lineages in this species.

Insects and collection sites
Triatoma garciabesi male bugs were collected from 24 populations found in Argentina and Paraguay (Table 1; Fig. 1).We included only males in the study because their sample size was larger than that of females and because sexual dimorphism has been observed in some populations from Argentina (Fiad et al. unpublished data).The collected specimens were first identified as males based on morphological features and then the results were confirmed by molecular analyses for one or two individuals from each population (see section DNA sequence analysis).The localities included in the study cover almost all T. garciabesi's distribution area (Fig. 1).All populations were collected from the field, except for one colony, Balbuena, which consisted of a population reared for one generation in the laboratory (Table 1).Individuals from Balbuena colony were supplied by the "Unidad Operativa de Vectores y Ambiente (UnOVE)", "Centro Nacional de Diagnóstico e Investigación en Endemo-Epidemias (CeNDIE)", which is part of the "Administración Nacional de Laboratorios e Institutos de Salud 'Dr.Carlos Malbrán' (ANLIS), from Santa María de Punilla, Córdoba, Argentina.In the field, T. garciabesi populations were collected from peridomestic structures, except for those from "Reserva Natural Bosques Telteca", Mendoza province, Argentina, which were collected from sylvatic habitats.The distribution of T. garciabesi encompasses five ecoregions (Fig. 1): Humid Chaco, Dry Chaco, High Monte, Low Monte and Paraná Flooded Savanna (the latter region contains 1 population).Ecoregions are considered to be terrestrial units that contain a geographically distinct set of species, natural communities and environmental conditions [46].The Humid Chaco ecoregion is an extremely flat plain, with very gentle slopes oriented from west to east.The climate is humid temperate, and rainfall follows a marked longitudinal gradient, with maximum records (exceeding 1300 mm rainfall) in the east [47].The Dry Chaco ecoregion is characterised by a warm subtropical continental climate, with annual rainfall concentrated in the summer and ranging between 500 and 700 mm [48].The climate in the High Monte ecoregion is temperate-arid, with annual rainfall of up to 200 mm, while the Low Monte ecoregion is also a temperate-arid climate, with annual rainfall of about 100 mm, occasionally reaching up to 200 mm [48].

DNA sequence analysis
A fragment of the mitochondrial coI of 34 insects collected from the 24 collection sites was analysed.Due to incongruences and misclassifications within the Sordida subcomplex species in the GenBank nucleotide database, a number of specimens from other Triatoma species were taken into account: T. sordida, T. sordida La Paz and T. rosai individuals; these had previously been determined by Panzera et al. [7] using chromosomal markers.In addition, newly available sequences from T. guasayana Genbank Acc.Nos.OR650776 to OR650782, also determined by Panzera et al [7], and a sequence from Triatoma rubrovaria Blanchard, 1843 Genbank Acc.No. AF021206 were used as outgroups to the Sordida subcomplex.Total DNA was extracted from the legs fixed in 70% ethanol (3 per specimen) using a standard phenol-chloroform technique procedure.A approximately 624-bp coI fragment was amplified by PCR using primers coIb (5′-TAA GCG TCT GGG TAG TCT GAR TAK CG-3′) [49] and ACO (5′-CTT GCA GGA GGA GGA GAY CC-3′) [50].The cycling parameters for the coI gene consisted of 5 min at 95 °C, followed by amplification for 30/35 cycles of 30 s at 95 °C, 30 s at 58 °C/47 °C and 1 min at 72 °C and terminated by 10 min at 72 °C, using the Taq ® 2× Master Mix Kit (New England Biolabs, Beverly, MA, USA).The PCR products were sent to Macrogen Inc. (Seoul, Korea) for DNA purification and subsequent sequencing.Both sequence strands were aligned and manually curated by chromatogram evaluation using Chromas (https:// techn elysi um.com.au/ wp/ chrom as/), and then deposited in the GenBank database (http:// www.ncbi.nlm.nih.gov), under accession numbers OR650712 to OR650782.Sequence alignment was performed using MAFFT v7.310 [51].The R package phangorn [52] was used to estimate a phylogenetic tree using the maximum likelihood (ML) method.The tree was bootstrapped 100 times to assess node support.The phangorn JModeltest function was used to find the best-fitting substitution model under the Bayesian Information Criterion (BIC).The phylogenetic trees were visualised and edited in ggtree [53].The R core base version 4.3.0 ( R Foundation for Statistical Computing, Vienna, Austria), the vegan (https:// CRAN.R-proje ct.org/ packa ge= vegan), ape v5.7.1 [54] and pegas v1.2 [55] packages were used to calculate the within and intergroup pairwise Kimura two parameter (K-2p) distances [55].A table and a heatmap were produced with the distance values.In addition, using the K-2p distance matrix, we conducted a multidimensional scaling (MDS) analysis with the base R package stats.The resulting MDS configuration was plotted in a two-dimensional scatter plot.The x and y coordinates of each point in the plot represented the positions of the individuals in the reduced space, with the distances between the points reflecting their genetic dissimilarities [56].To determine if the dataset exhibited distinct groupings or clusters based on genetic distances, a clustering analysis was performed on the MDS plot using k-means clustering, with the function k-means of the R stats package.
In addition, a minimum spanning network (MSN) haplotype network was constructed using PopART software [57].

Geometric morphometry
Individuals were grouped according to lineages based on their results of the coI genetic analyses.A total of 198 individuals belonging to the 24 populations were included in the morphometry analysis (Table 1).For all individuals, digital images of the dorsal view of the right forewing and pronotum and of the ventral view of head were taken using a digital camera (model S9900; Nikon Corp., Tokyo, Japan) mounted on a stereomicroscope (model Stemi SV-11; Carl Zeiss AG, Jena, Germany) at 6× magnification.All images included a reference scale.The maximum forewing length and the total forewing area, which were used to estimate the wing aspect ratio, were measured using ImageJ software version 1.52 (https:// imagej.nih.gov/ ij/).The wing aspect ratio (4 × [wing length 2 /wing area]; [40]) was compared among lineages using one-way analysis of variance (ANOVA).A quantitative shape analysis was also performed using geometric morphometry based on the statistical analysis of landmark coordinates for the forewing, membranous Filled circles indicate the location of the studied populations, with dark-purple circles corresponding to the populations assigned to the Eastern lineage by coI analysis and lilac circles corresponding to populations assigned to the Western lineage by coI analysis.Numbers correspond to the population code in Table 1.Ecoregions are according to Dinerstein et al. [88] (https:// ecore gions.appsp ot.com/).coI, Cytochrome c oxidase I gene and stiff portion of the forewing, the head and the pronotum.Landmarks are point locations that are biologically homologous across specimens [58].Semilandmarks were used to analyse the lateral contour of both the stiff and membranous parts of the forewing and the contour of the pronotum.These parts do not exhibit sufficient homologous points to be used to characterise the morphology of the specimens using traditional landmarks; therefore, the contour of these parts using semilandmarks were quantified [59].To analyse forewing variation, we included eight type I landmarks of the right forewing (Fig. 2a).Landmarks 1, 2, 4, 5 and 6 are positions between the stiff and membranous portions of the hemelytra.The stiff and membranous parts of the forewing were analysed separately.A combination of landmarks and semilandmarks were used to characterise the contour of these portions of the forewing.For the stiff portion, seven type I landmarks were used to characterise the shape of the boundary between the stiff and membranous parts, and 13 equidistant semilandmarks were used to characterise the lateral contours of this part (Fig. 2b).For the membranous part, six type I landmarks were used to characterise the shape of the boundary between the stiff and membranous parts, and 16 equidistant semilandmarks were used to characterise the lateral contours (Fig. 2b).For the head, seven coplanar type II landmarks on the right side of the ventral view of the head were defined and collected; these landmarks included the anterocular, ocular, postocular and neck regions (Fig. 2c).For the pronotum, one coplanar type II landmark in the right anterolateral angle and 13 equidistant semilandmarks on the contour of the right side of the structure were defined and collected (Fig. 2d).The pronotum and head from individuals belonging to the Aguirre and Balbuena populations were not included in the analysis because only the forewings from these populations were available.
Semilandmarks were transformed to landmarks and analysed together with traditional landmarks [59].We employed TPS DIG version 2.31 [60] to record landmarks and semilandmarks of the described structures.Data on the head, pronotum and forewing shape were extracted with a generalised full Procrustes fit and a projection to shape tangent space [61].Procrustes coordinates as shape variables were used.Centroid size (CS; i.e. the square root of the sum of the squared distances from each landmark to the centroid of the configuration) was computed as a measure of size for each structure or part of the structure.These steps of morphometry analysis were performed using the free software MorphoJ version 1.07a [62].Geometric variables of the forewing, portions of the forewing, the head and the pronotum were analysed separately.Variations across populations for the CS between Eastern and Western lineages for the forewing, part of the forewing, the head and the pronotum were explored, after testing for normality of the CS data with the Shapiro-Wilks, one-way ANOVA or Kruskal-Wallis tests using the software InfoStat version 2016 [63].For shape measurements of the head, pronotum and forewings, we calculated Procrustes distances between lineages and evaluated the significance of these distances via a non-parametric test based on permutations (1000 runs) using MorphoJ [62].We then performed discriminant function analysis (DFA) and calculated the percentage of phenotypic similarity between lineages using the crosscheck test of discriminant analysis in InfoStat software [63].Allometric relationships between shape and CS for the different structures or portions of the structure studied were checked using a multivariate regression of Procrustes coordinates on CS.The studied structures and the stiff portion of the wing showed no significant allometric effects according to the regression analyses.The membranous portion of the wing showed a significant allometric relationship (P = 0.0027), but with small allometric effects (2.30%).

Divergence of coI
The GTR + G + I substitution model was selected to construct the ML phylogenetic tree based on coI (Fig. 3).The tree shows five lineages within the ingroup, plus T. rubrovaria and T. guasayana individuals as outgroups, which are indicated in grey and black font in Fig. 3.The dendrogram drawn from the heatmap (Additional file 1: Figure S1) and the MDS plot (Additional file 2: Figure S2) reflect the same groupings as the ML tree.Mean distances within and between the lineages also indicate that lineages are well separated (Additional file 3: Table S1; Tables 2, 3).Average intraclade and interclade divergence ranged from 0.3% to 2.5% and from 3.2% to 16.3%, respectively.The greatest interspecific divergence was found between T. garciabesi Western lineage and other lineages, with the greatest divergence with T. guasayana (16.2%), followed by that with T. sordida sensu stricto (T.sordida s.s.; 12.6%), T. rosai (8.0%) and T. sordida La Paz (5.8%).
As expected, the shortest pairwise genetic distance was observed between both T. garciabesi lineages (mean 3.2%, range 2.3-3.8%)(Additional file 3: Table S1; Tables 2, 3).Within T. garciabesi, the Western lineage showed the highest genetic distance among individuals, which also had the highest values for haplotype and nucleotide diversity (Table 3; Additional file 3: Table S1).The haplotype network also supports the delimitation of the lineages detected by genetic distance and on the tree topology.All lineages are separated by several mutation steps.Furthermore, the network shows the great diversity of the T. garciabesi Western lineage (Fig. 4).

Forewing size, wing aspect ratio and size of the stiff and membranous portions of the forewing differed between coI-defined lineages
Forewing CS varied significantly between the Eastern and Western lineages (F (2, 197) = 6.929,P = 0.0092), with forewings from the Eastern lineage being larger than those from the Western lineage (mean = 7.741, standard deviation [SD] = 0.475 vs mean = 7.543, SD = 0.468 for Eastern vs Western lineages, respectively).Wing aspect ratio also varied significantly between lineages (F (2, 197) = 140.31,P < 0.0001), being bigger for the Western lineage (mean = 16.712,SD = 1.539 vs mean = 20.537,SD = 2.119 for Eastern vs Western lineages, respectively).The size of the stiff portion of the forewing for the Eastern lineage was significantly larger than that of the Western lineage (F (2, 197) = 8.967, P = 0.0031) (mean = 9.442, SD = 0.708 vs mean = 9.125, SD = 0.633 for Eastern vs Western lineages, respectively).No significant differences between lineages were found for the size of the membranous portion (F (2, 197) = 0.0555, P = 0.8141).

Forewing shape differed between coI-defined lineages
The first discriminant factor of the DFA explained 100% of the total variation.The Eastern and Western lineages were well differentiated along the first DFA axis (Fig. 5).Both Mahalanobis and Procrustes distances were significantly different between lineages (Mahalanobis distance = 2.172, P < 0.0001; Procrustes distance = 0.0279, P < 0.0001,).Forewing shape changes, corresponding to extreme positions of the first axis of the DFA, involved all landmark configurations except for landmark 8. Changes in forewing shape in the Eastern lineage compared with the Western lineage involved a more developed stiff portion and a less developed membranous portion (Fig. 5a).Analysis of the forewings of the Eastern and Western lineages revealed 14.65% misclassified individuals; of all these misclassified individuals, 55.20% corresponded to the Eastern lineage and 44.8% to the Western lineage.When the populations with wrongly assigned individuals based on forewing shape were considered, ten populations (5 from the Eastern lineage and 5 from the Western lineage) had individuals assigned to the other lineage.Most of these populations (7 of 9 populations) were located at the distribution boundaries between the two lineages (Table 4; Fig. 1).

Stiff and membranous portions of the forewing differed between coI-defined lineages
The DFA for the shape of the contour of the stiff and membranous portions of the forewing showed that the Eastern and Western lineages were well differentiated (Fig. 5b, c).The first discriminant factor of both DFAs explained 100% of the total variation.Both Mahalanobis and Procrustes distances were significantly different between lineages for both portions of the forewing (stiff portion: Mahalanobis distance = 3.302, P < 0.0001 and Procrustes distance = 0.0271, P < 0.0001; membraneous portion: Mahalanobis distance = 3.285, P < 0.0001 and Procrustes distance = 0.0169, P < 0.0001).Changes in the shape of the contour of the stiff portion showed that the Eastern lineage had a wider shape than the Western lineage, mainly involving the basal and proximal parts of this portion of the forewing (Fig. 5b).For the shape of the contour of the membranous portion, the Western lineage exhibited a longer and narrower shape of the contour than the Eastern lineage (Fig. 5c).

Pronotum size differed between coI-defined lineages
A Kruskal-Wallis's test of the pronotum CS showed significant differences between the Eastern and Western lineages (H = 5.82, P = 0.016).Pronotum size was larger in individuals from the Eastern lineage than in those from the Western lineage (Dunn post-hoc test, P < 0.001) (mean = 1.508,SD = 0.090 vs mean = 1.488,SD = 0.116 for Eastern vs Western lineages, respectively).

Pronotum shape differed between coI-defined lineages
The first DFA axis accumulated 100% of the total variation for the forewing.Both lineages were defined in the space of the first DFA (Fig. 6a).Mahalanobis distance was significantly different between lineages (Mahalanobis distance = 2.253, P < 0.0001), whereas Procrustes distance was not significantly different between lineages (Procrustes distances = 0.0153, P = 0.0679).Changes in pronotum shape involved mostly landmarks located in the anterolateral angle and those located in the humeral angle (Fig. 6a).Analysis of pronotum shape in the Eastern and Western lineages showed 13.53% misclassified individuals.The analysis including the populations with wrongly assigned individuals showed that 16 populations (6 from the Eastern lineage and 10 from the Western lineage) were misclassified.Regarding pronotum shape, most of the populations that exhibited misclassified individuals were located at the border of the distribution area of both coI lineages (9 of 11 populations) (Table 5; Fig. 1).

Head shape-but not size-differed between coI-defined lineages
No significant differences in head CS between Eastern and Western lineages were observed (ANOVA, F (2, 173) = 0.79, P = 0.3767).For wings, the first discriminant factor of the DFA explained 100% of the total variation.Eastern and Western lineages were differentiated along the first DFA axis (Fig. 6b).There were significant  differences in both the Mahalanobis and Procrustes distances between lineages (Mahalanobis distance = 1.337,P < 0.0001; Procrustes distance = 0.0136, P < 0.0001).Changes in head shape between lineages involved anterocular, ocular and neck regions (Figs. 2, 6b).Heads of individuals from the Western lineages exhibited a widening of the anterior segment and greater eye development than heads of individuals of the Eastern lineages.DFA analysis of head shape in the Eastern and Western lineages showed 24.00% of misclassified individuals.When the populations with wrongly assigned individuals were included in the analysis, 16 populations (all populations from Eastern lineage and 10 from Western lineage) showed misclassified individuals (Table 6).Half of the populations (8 of 16 populations) that had misclassified individuals were located at the border of the distribution area of both coI lineages.

Genetic divergence within T. garciabesi
The results obtained in this study demonstrate the existence of two separate lineages within T. garciabesi, here referred to as the Eastern and Western lineages.The lineages showed a comparatively different non-overlapping geographical distribution, with the Eastern lineage being restricted to and clustered in a part of the centre and east of the Argentine provinces of Chaco and Formosa, which coincides with the Humid Chaco ecoregion (there was only 1 population at the boundary between the Humid and the Dry Chaco and 1 population located in the Paraná Flooded Savanna; Fig. 1).The Western lineage is present throughout the rest of the species distribution range.Most populations are located in the Dry Chaco ecoregion, and two populations are located in the High and Low Monte ecoregions.The Dry Chaco ecoregion exhibits a marked longitudinal precipitation gradient from east to west and marked differences in phytogeographical characteristics [48], suggesting that lineage diversification could be associated with landscape characteristics.Studies using mitochondrial genes (cytb, coI and coII) for phylogenetic separation of closely related species have been carried out in other triatomine species (for review see [4]), such as, for example, Rhodnius prolixus Stål, 1859 and Rhodnius robustus Larrousse, 1927 sensu lato [64], Triatoma rubida Uhler, 1894 and Triatoma recurva Stål, 1868 [65], species of the Triatoma brasiliensis complex [66] and species of the Mepraia genus [67].These genes have also been used in population analyses of Triatoma sanguisuga Leconte, 1855 [68] and T. infestans [69].In the present study, the inference of the phylogenetic tree constructed using the ML method allowed us to classify the samples into seven clades; thus, the populations were separated into: T. garciabesi from the Eastern lineage, T. garciabesi from the Western lineage, T. Fig. 4 Minimum spanning network indicating the same clades as those of the maximum likelihood tree (Fig. 3).Nodes represent the haplotypes (see Fig. 3 for color coding), with their size proportional to the frequency of the haplotype.The numbers of mutational steps separating haplotypes are represented by the short vertical lines along each branch of the network.Individuals within haplotypes are depicted in Additional file 4: table S2.
guasayana, T. rosai, T. rubrovaria, T. sordida La Paz and T. sordida s.s.The calculation of coI gene pairwise distance is generally used to predict cryptic or novel species ( [28] and references cited therein).In general, intraspecific genetic diversity in mitochondrial genes is rarely > 2%, and mostly < 1% [27,70,71].Meier et al. [70] indicate that if the objective is to predict cryptic species based on genetic distance, the smallest interspecific distance has to be used.If the minimum interspecific genetic value distance from congener species is ≥ 2%, it is possible to avoid overestimating species diversity based on empirical thresholds.However, in studies on triatomines, different genes and regions within the genes have been used.Therefore, different consensuses were taken into account to decide whether the variation belonged to different species or to populations within the same species.For example, in one study, the mean pairwise distance between species for the coI sequence was about 7% (e.g.Mepraia gajardoi Frias, Henry & Gonzalez, 1998 and M. spinolai [67]).For the genus Triatoma, the largest mean pairwise nucleotide distance was found to be 5.2% between North and Centre-South Triatoma patagonica Del Ponte, 1929 lineages [72].Similar values between pairs of closely related species (T.infestans-Triatoma platensis Neiva 1913 and T. sordida-T.garciabesi) have been reported for other mitochondrial gene fragments (cytb) [66].For two T. brasiliensis populations, the pairwise distance of > 8.4% found between them was one of the arguments that supported their placement to the species level.Our results for T. garciabesi lineages showed a mean pairwise distance between the Eastern and Western lineages of 3.2%.Hence, the coI genetic distances observed for T. garciabesi lineages would not be sufficient to consider them as distinct species.However, the environmental conditions to which lineages are exposed would play a crucial role in increasing genetic differentiation through divergent selection, a pattern known as isolation-by-environment [73,74].

Congruence between genetic and morphometric variation
In agreement with the diversity that has been observed in coI, our results showed morphometry differentiation between the two T. garciabesi lineages in the forewing, head and pronotum measurements.Patterns of variations Fig. 5 Frequency distribution of the first axis of the discriminant function analysis (DFA) for forewing shape (a) and for the stiff (b) and membranous (c) portions of the forewing for the Eastern and Western lineages of T. garciabesi.The shape changes associated with the first axis of the DFA are visualised as configurations corresponding to extreme positions of the axis.Black, configuration for negative extreme scores; grey, configuration for positive extreme scores in size or shape have been interpreted previously in terms of their evolutionary implications [75].It is known that environmental factors (e.g.temperature and humidity) determine phenotypic characteristics in several triatomine species (reviewed in [76]).Size differentiation of morphological attributes has been suggested to be strongly (but not exclusively) influenced by the ecological characteristics of species and populations [76].If we consider that genetic and historical attributes are often thought to affect shape variation of morphological  attributes [77], we can assume that the consistent size differentiation observed in the studied morphological attributes of T. garciabesi lineages is due to ecological differences between the ecoregions where the T. garciabesi Eastern and Western lineages occur.Individuals from the Eastern lineage showed a greater size, both of heads and pronota, than individuals from the Western lineage.In other triatomine species, forewing size has been found to be a good predictor of body length [78].Moreover, total body length has been found to be closely correlated with total body weight and with female fecundity [79].These results may indicate that genetic variability within populations of T. garciabesi plays a role in size variation in these populations.However, this is a preliminary suggestion, and it may be necessary to carry out additional analyses, experiments or studies to confirm and better understand this association between genetic variability and size variation.Differences in shape variation between lineages suggest genetic variation.Interestingly, the discriminant analysis of shape differentiation showed higher percentages of misclassified individuals in populations located at the distribution boundaries of both clades than in the remaining populations.This pattern was observed mainly in the shape of forewings and pronotum, suggesting the presence of a transition zone or ecotone where the morphological characteristics are in transition between the two lineages.This could be an area of hybridisation or a zone of genetic contact.Phenotypic diversity is often known to be higher at habitat boundaries than in more homogeneous interior regions of the habitats [80].Our results also suggest adaptive divergence associated with habitat characteristics, but the degree of reproductive isolation of the two lineages is still unknown.The mechanisms that limit gene flow between populations can be studied by making experimental crosses and analysing the offspring [81].Further studies to elucidate a possible reproductive isolation between the lineages may help to understand whether genetic differentiation results in some level of reproductive isolation.The association between genetic and morphological variation has been studied previously in triatomines [25,72,82,83], providing data that improve current understanding of the different speciation processes that may be operating in this subfamily.For Rhodnius pallescens Barber 1932 lineages, diversification was found to be mainly influenced by biogeographical events [80], while for T. sordida, the differentiation between populations from eastern and western regions of Paraguay was suggested by the authors of the study to be associated with eco-geographical isolation by distance [83].In T. patagonica, genetic and morphological diversification were explained by an ecological diversification accompanied with cytogenetic differentiation and partial reproductive isolation of the most distant populations [72].

Implications of variation in flight-related traits and mitochondrial lineage distribution
Dispersal abilities are of major importance in establishing and rearranging the geographical distribution of species [29].The coI mitochondrial gene, as all protein coding genes in the mitochondrial genome, is considered to be under purifying selection, since it constantly sweeps away deleterious nonsynonymous mutations [84].On the other hand, morphometry of flight-related traits is more likely to have evolved under adaptative selection.In particular, the forewing parts of the hemelytra are expected to be under strong selective pressure [36].The morphometry characteristics of these flight-related traits and the parts of the hemelytra of this species suggest which forewing morphotype is selected in each lineage in relation to the environment where they develop.These characteristics would determine the type of flight of each lineage, which would improve our understanding of how flight-related traits shape the geographical distribution of genetic diversity in this species.Long and narrow forewings give greater deformation capacity at a lower energy cost for long distance flights.Despite this assertion, no specific tests of flight capacity have been carried out in triatomines.While previous studies on other insects may provide indirect support, the morphology of triatomine hemelytra is different and this type of study was necessary to properly interpret these results.The aspect ratio indicates the flight capacity, with higher values of this ratio being compatible with more efficient and longer flights [38,85].In dragonflies (Order: Odonata), the functional roles of forewing structural components suggest that forewings require a balance between flexibility and stiffness [86].Our results showed that the portions of the forewings of the two lineages of T. garciabesi showed differences in both shape and aspect ratio.The forewing of the Eastern lineage was larger than that of the Western lineage; however, when variation was broken down into the stiff and membranous portions, only the stiff portion was larger in the Eastern lineage while the membranous portion did not differ between lineages.The analysis of shape variation showed that the stiff portion of individuals of the Eastern lineage is wider than that those of the Western lineage, whereas the membranous part of the Western lineage is longer and narrower than that of the Eastern lineage.The aspect ratio varied between lineages, with individuals of the Western lineage exhibiting higher values.The morphological differences observed in the hemelytra of the two lineages suggest different types of flight.The Western lineage would have more aerodynamically stable, less energetically costly and longer flights.The morphological attributes of the Eastern lineage are compatible with short, slow and somewhat unstable flights.Hence, these morphological attributes would determine different types of flight associated with the lineages.Among the characteristics that define the Humid Chaco ecoregion, in addition to the abundant rainfall, which favours vegetation diversity, there is little thermal amplitude, and the local fauna is abundant [47,87].These characteristics would allow a long time window for the development of triatomine populations; similarly, the great availability of local fauna means a reduced need to disperse in search of hosts.Thus, the forewing attributes corresponding to short and agile flights could be an adaptation of individuals from the Humid Chaco region that would allow short and agile flights, and a consequent increased allocation of resources to reproduction.On the other hand, in the arid Dry Chaco and Monte ecoregions, the landscape is less lush and host availability is lower than in the Humid Chaco.These landscape characteristics would require longer and aerodynamically more stable flights, which would favour the search for hosts.It should be emphasised here that this study does not directly measure flight capacity and that these suggestions are based on morphological characteristics and extrapolated to the potential of the flight behaviours of this triatomine species.Dispersal capacity experiments would be necessary to validate and strengthen these statements.
Phenotypic variation in head shape also differed between lineages.Heads from the Western lineages exhibited a widening of the anterior segment and greater eye development than those of the Eastern lineages.These lineage-related changes in head conformation are consistent with the changes in forewing characteristics described above.In M. spinolai, macropteran males were found to have a greater convexity of the compound eyes and a greater interocular distance than micropteran males, which could imply that the former have a greater navigational capacity and sense of orientation during flight [41].The characteristics of the forewings and heads help to understand the possible type of flight observed in each lineage.The types of flight are favoured by selective pressures linked to the landscape in which the lineage evolved.The pronotum does not follow the expected shape changes in relation to these characteristics, since in the Eastern lineage, the pronotum exhibits a greater development of the humeral angle, which suggests a greater development of the alar muscle.
In the present study, we were unable to molecularly characterise all individuals and assumed that all individuals from the same population belonged to the same species as when more than one individual from the same population were molecularly characterised, they were found to belong to the same species.This limitation does not undermine our conclusions since the pronounced genetic divergence is consistent with morphometry variation.

Conclusions
Our results provide preliminary insights into the evolution and diversification of T. garciabesi.There is a pronounced coI mtDNA divergence, with the existence of two lineages with different distribution ranges.Consistent with genetic divergence, morphometry variation was clustered according to lineages, with congruent variation in size and shape of the forewing, pronotum and head.

Fig. 1
Fig.1Distribution area of Triatoma garciabesi based on Ceccarelli et al.[24] and extended to the Paraguayan populations included in this study.Filled circles indicate the location of the studied populations, with dark-purple circles corresponding to the populations assigned to the Eastern lineage by coI analysis and lilac circles corresponding to populations assigned to the Western lineage by coI analysis.Numbers correspond to the population code in Table1.Ecoregions are according to Dinerstein et al.[88] (https:// ecore gions.appsp ot.com/).coI, Cytochrome c oxidase I gene

Fig. 2 a
Fig. 2 a Landmark (open circles) positions for the stiff (brown) and membranous (white) portions of the forewing.b Landmark (open circles) and semilandmark (grey-filled circles) positions for the stiff (brown) and membranous (white) portions of the forewing.c Landmark positions (open circles) for the head.d Landmark (open circles) and semilandmark (grey-filled circles) positions for the pronotum contour, for the studied Triatoma garciabesi individuals

Fig. 3
Fig. 3 Maximum likelihood phylogenetic tree based on coI gene fragments for Triatoma garciabesi.Triatoma rubroviaria, T. sordida s.s, T. rosai and T. sordida from La Paz were included as outgroups.Tree topology reveals that T. garciabesi is composed of two clades, represented by Eastern and Western lineages, respectively.Grey-and black-filled circles above nodes represent statistical support obtained through bootstrap replications, 0.7 and 0.9, respectively.Letters included before the name of the population indicate country of origin: AR, Argentina; BL, Bolivia; BR, Brazil; PY, Paraguay.coI, Cytochrome c oxidase I gene; s.s., sensu stricto

Fig. 6
Fig. 6 Frequency distribution of the first axis of the discriminant function analysis (DFA) for pronotum (a) and head (b) shape for T. garciabesi Eastern and Western lineages.The changes in shape associated with the first axis of the DFA are visualised as configurations corresponding to extreme positions of the axis.Black indicates the configuration for negative extreme scores; grey indicates the configuration for positive extreme scores

Table 1
Genetic lineages of Triatoma garciabesi populations defined using the cytochrome c oxidase I gene

Table 2
Mean Kimura's two-parameter substitution model (K-2P)-based pairwise genetic distances between Triatoma garciabesi Eastern and Western lineages, and for the outgroups Triatoma guasayana, T. rubrovaria, T. rosai and T. sordida sensu stricto and T. sordida La Paz, for the cytochrome c oxidase I gene fragment Values in table are genetic distances calculated using the K-2P method for the cytochrome c oxidase I gene fragment

Table 3
Population genetic diversity indices and demographic history test estimated with cytochrome c oxidase I gene sequences for T. garciabesi Eastern and Western lineages, and for the outgroups T. guasayana, T. rosai, T. rubrovaria and T. sordida sensu stricto and T. sordida La Paz N, Number of individuals; S, segregating sites; Nh, number of haplotypes; Hd, haplotype diversity; pi, nucleotide diversity; Tajima's D, value of Tajima's neutrality test; TajimaPval, P-value of Tajima's neutrality test

Table 4
Number of wrongly assigned individuals derived from discriminant function analysis for variation in T. garciabesi wing shape

Table 5
Number of wrongly assigned individuals derived from discriminant function analysis for variation in T. garciabesi pronotum shape

Table 6
Number of wrongly assigned individuals derived from discriminant function analysis for variation in T. garciabesi head shape