Males of T. garciabesi from 24 populations across Argentina and Paraguay were included in this study, covering nearly the entire distribution range ( Table 1 ). Owing to morphological similarities between T. garciabesi, Triatoma sordida (Stål, 1859), and Triatoma rosai Alevi, Oliveira, Garcia, Cristal, Delgado, Bittinelli, Reis, Ravazi, Oliveira, Galvăo, Azeredo-Oliveira, Madeira, 2020, male identification was confirmed using molecular techniques for one or two individuals from each population. For this, coI gene was used, and total DNA was extracted from the legs fixed in 70% ethanol using a standard phenol–chloroform technique procedure. An approximately 624-bp coI fragment was amplified by PCR and the 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://technelysium.com.au/wp/chromas/). A phylogenetic tree using the maximum likelihood (ML) method was done and the best-fitting substitution model under the Bayesian Information Criterion (BIC) was used (for more details see 26). The origin (field or colony reared in the laboratory) and the collection site (peridomestic or sylvatic) of the populations included in this study were described elsewhere (26). Only males were included in this study due to the sexual dimorphism observed in morphometric traits for this species (19). T. garciabesi distribution encompasses five ecoregions where the Eastern and Western lineages are found (26). The Eastern lineage is restricted to a portion of the center and eastern part of the Argentine provinces of Chaco and Formosa, predominantly within the Humid Chaco ecoregion, with one population near the boundary between the Humid and Dry Chaco and another in the Paraná Flooded Savanna ( Figure 1 ). The Western lineage is present throughout the remainder of the species distribution range, with most of the populations being located in the Dry Chaco ecoregion with a few populations in the High and Low Monte ecoregions. These ecoregions exhibited a marked precipitation variation and significant differences in phytogeographical characteristics (30).

A total of 198 males from 24 populations were analyzed, each with a sample size of between 5 and 21 male specimens ( Table 1 26). Digital images of the ventral view of the head and the dorsal view of the right forewing were captured to analyze the size and shape variables of flight-related traits, specifically the head and forewing structures ( Figure 2 ). These traits were further characterized by assessing both stiff and membranous portions of the forewing, as outlined elsewhere [ Figure 2 , Verly et al. (26)]. All methodologies and techniques used for this study were previously explained in Verly et al. (26), which provides further details on the analytical framework and the population-specific data.

To obtain the geographic variables of each population, latitude, longitude, and altitude were recorded using Google Earth Pro. Climatic and vegetation data were extracted using ERAS 5 (Copernicus Climate Data Store) (31) at a spatial resolution of 28 × 28 km, covering a 10-year period (2013–2023). For this study, we included only variables considered influential in the spatial distribution of Triatominae or relevant to flight dispersal in insect dispersal (13, 25, 32–35).

The climatic variables considered were (1) annual mean temperature, (2) maximum temperature of the warmest month, (3) minimum temperature of the coldest month, (4) wind speed, (5) wind direction, (6) annual rainfall, and (7) relative humidity. Wind direction and velocity were extracted from the North–South and East–West components of winds. Relative humidity was calculated as: RH = 100 - 5 × (T − T _d), where T is the average monthly temperature and T _d is the average monthly dew point temperature (36). The vegetation indexes included were (1) the Normalized Difference Vegetation Index (NDVI), (2) leaf area index low vegetation (LAILV), and (3) leaf area index high vegetation (LAIHV) (31). The leaf area index (LAI) quantifies the amount of leaf material in a canopy. By definition, it is the ratio of one-sided leaf area per unit ground area (31).

Statistical analyses and results were first conducted by including all 24 populations in the study (i.e., species level of analysis) and then separately for the Eastern and Western lineage (i.e., lineages level of analysis). Each structure or part of the structure was analyzed separately. The landmark configurations were superimposed using a generalized Procrustes analysis (GPA) (37).

To test the populations for differences in the average size of the forewing, stiff, and membranous portions and head, we performed analysis of variance (ANOVA) followed by Tukey’s post-hoc tests. For these analyses, the software InfoStat version 2016 was used (38). To assess the degree of correlation between size and shape for each structure or part of structure, we fitted a multiple linear regression between the Procrustes coordinates and the centroid size. None of the structures showed an allometric effect, except for the membranous portion of the forewing that showed a significant allometric relationship (p = 0.0027), but with small allometric effects (r = 0.02). To test population differentiation for each structure or part of structure, a canonical variate analysis (CVA) was done. We calculated Procrustes distances between pairs of populations and evaluated their significance via a non-parametric test based on permutations (1,000 runs). These steps of analysis were done with MorphoJ version 1.07a (39). The Procrustes distances were represented in unrooted neighbor-joining (NJ) trees using the free software MEGA X 10.2.6 (40). We calculated bootstrap values for 1,000 replications following Ascarranuz et al. (41). The association of geographic distances with linear and geometric morphometric measurements of forewing, stiff, and membranous portions and head was analyzed using Mantel tests with the free software PASSAGE 2 2.0.11.6 (42). Morphological distance matrices between populations were constructed using Procrustes distances for shape measurements and Euclidean distances for size measurements. We tested the effect of geographic, climate, and vegetation cover on phenotypic changes of forewing, stiff, and membranous portions and head using a two-block partial least squares (PLS) analysis. PLS analysis tests for covariation between blocks of variables (shape or size of each structure within populations) were done using the Procrustes coordinates for shape or CS for size with geographic and climatic variables and vegetation indexes (43). To calculate a p-value and significance test, we generated 9,999 permutations for each comparison. PLS analyses were calculated using the R package geomorph (44) using the two.b.pls function. Additionally, to explore if there are associations between the NDVI, as a measure of habitat heterogeneity (45, 46), and shape variations in flight-related traits, we performed a Procrustes regression with permutations for each structure including the lineage and the populations nested within the lineage in the analysis. These analyses were conducted using the procD.lm function from the geomorph package in R.

At the species distribution range level, the size component analysis revealed no significant differences in forewing size between populations [F _{(197, 23)} = 826.94, p = 0.9907]. However, the membranous and stiff portions of the forewing did show significant differences [F _{(197, 23)} = 8.00, p < 0.0001; F _{(197, 23)} = 9.55, p < 0.0001] ( Supplementary Material 1 ). Head size did not show significant differences between populations [F _{(174, 23)} = 20.93, p = 0.1217].

The ordination of the populations by CVA for forewing showed that the first and second canonical variables jointly explained 52.48% of the total variation (CV1: 34.31%, CV2: 18.17%). CGB was the only population with significantly different Procrustes distances compared to the others ( Supplementary Material 2 ). The forewing NJ tree, based on Procrustes distances, revealed that only four nodes had bootstrap support greater than 70%. One of these nodes grouped together two populations from the Eastern lineage (MA and LE) ( Figure 3A ). For the membranous portion of the forewing, the CVA revealed that the first two axes explained 39.27% of the total variation (CV1: 24.65%, CV2: 14.62%). All pairwise Procrustes distances between populations were significantly different from zero (p < 0.0001, Supplementary Material 2 ). The NJ tree based on Procrustes distances indicated 8 from 21 nodes that were supported by bootstrap (values ≥ 70%). One of these nodes grouped two populations from the Eastern lineage, while five nodes grouped Western populations together ( Figure 3B ). Regarding the stiff portion of the forewing, the first two axes explained 45.03% (CV1: 28.00%, CV2: 17.03%). As with the membranous portion, all pairwise Procrustes distances between populations were significant (p < 0.0001, Supplementary Material 2 ). The NJ tree based on Procrustes distances showed that 7 out of 21 nodes were bootstrap-supported. One of these nodes grouped two populations from the Eastern lineage, while the other nodes grouped Western populations together ( Figure 3C ). For head, the first two discriminant factors from the CVA jointly explained 49.92% of the total variation (CV1: 28.42%, CV2: 21.50%). CGB had the highest number of significant Procrustes distances with 18 out of 21 pairwise comparisons ( Supplementary Material 2 ). The NJ tree based on Procrustes distances indicated that 3 out of 20 clusters were supported by bootstrap values ( Figure 3D ).

The Eastern lineage showed significant differences in the size component of forewing [F _{(52, 5)} = 11.90, p < 0.0001], membranous, and stiff portions [F _{(52, 5)} = 16.73, p < 0.0001; F _{(52, 5)} = 12.37, p < 0.0001 for membranous and stiff portions, respectively] and head [F _{(52, 5)} = 6.44, p < 0.0001]. In all cases, CGB and CO exhibited the smallest structure or portion of structure significantly different from the other populations (Tukey post-hoc tests, p < 0.010) ( Supplementary Material 1 ). The first two discriminant factors of the forewing CVA jointly explained 85.72% of the total variation (CV1: 62.81%, CV2: 22.91%). The forewing NJ tree based on Procrustes distances showed no bootstrap-supported clades ( Figure 4A ). For the membranous portion of the forewing, the first two CVA axes explained 82.68% of the total variance (CV1: 66.47%, CV2: 16.21%). All pairwise Procrustes distances between populations were significantly different from zero (p < 0.0001, Supplementary Material 3 ). The NJ tree based on Procrustes distances showed that three out of four clades were bootstrap-supported. One of the bootstrap-supported clades separated the CGB population from Formosa province from the TI and LC populations in Chaco province. Other bootstrap-supported clades grouped the TI and LC populations from Chaco province together ( Figure 4B ). For the stiff portion of the forewing, the first two axes explained 71.62% of the total variation (CV1: 55.83%, CV2: 15.79%). All pairwise Procrustes distances between populations were significantly different (p < 0.0001, Supplementary Material 3 ). The NJ tree based on Procrustes distances showed one out of four bootstrap-supported clades ( Figure 4C ). For the head, the first two axes of the CVA explained 81.03% of the total variation (CV1: 62.53%, CV2: 18.49%). Procrustes distances for the head from the LE, LC, and MA populations did not differ significantly from any of the other lineage populations ( Supplementary Material 3 ). The NJ tree based on Procrustes distances did not show any bootstrap-supported clades ( Figure 4D ).

The Western lineage showed no significant size differentiation for the forewing and head between populations [F _{(125, 17)} = 0.43, p = 0.9771; F _{(125, 17)} = 0.15, p = 0.9320 for forewing and head size, respectively]. However, significant differences were found for the size component of the membranous portion of the forewing [F _{(125, 17)} = 5.90, p < 0.0001]. The CE and AV populations had the smallest membranous portions, while HI, CAS, and CAA had the largest membranous portion (Tukey post-hoc tests, p < 0.01, Supplementary Material 1 ). The stiff portion of the forewing also showed significant differences across populations [F _{(125, 17)} = 8.25, p < 0.0001], where AG had the smallest stiff portion and PY, CAA, CAS, and YO populations had the largest ones (Tukey post-hoc tests, p < 0.01, Supplementary Material 1 ). For the shape component of forewing variation, the first two discriminant factors of the CVA explained 50.20% of the total variation (CV1: 36.01%, CV2: 14.19%). Most of the populations showed significant Procrustes distances ( Supplementary Material 4 ), except for CAA, where 11 out of 17 Procrustes distance comparisons were not significant ( Supplementary Material 4 ). The NJ forewing tree based on Procrustes distances showed no bootstrap-supported clades ( Figure 5A ). For the membranous portion of the forewing, the first two CVA axes explained 43.77% of the total variance (CV1: 26.80%, CV2: 16.80%). All pairwise Procrustes distances between populations showed significant differences (p < 0.0001, Supplementary Material 3 ). The NJ tree based on Procrustes distances showed that 7 out of 16 clades were bootstrap-supported. One of the bootstrap-supported clades separated the RI population from Salta province from five of the six Paraguayan populations (YO, CAA, SAN, CAS, and CAN). Additionally, clades grouping the Paraguayan populations YO and CAA, as well as CAS and CAN, were also supported by bootstrap values ( Figure 5B ). For the stiff portion of the forewing, the first two axes explained 43.97% of the total variation (CV1: 27.86%, CV2: 15.94%). All pairwise Procrustes distances between populations were significantly different (p < 0.0001, Supplementary Material 3 ). For this portion of the forewing, the NJ tree based on Procrustes distances shows that 7 out of 16 clades were bootstrap-supported. One of the bootstrap-supported clades separated the RVP population from La Rioja province from two of the six Paraguayan populations (YO and SAN). In addition, the clade with YO and SAN was supported by bootstrap values ( Figure 5C ). For the head, the first two CVA axes explained 55.83% of the total variation (CV1: 32.89%, CV2: 22.95%). Approximately 65% of the Procrustes distance comparisons between population pairs showed no significant differences (p > 0.050, Supplementary Material 4 ). The NJ tree based on these Procrustes distances showed 2 out of 14 bootstrap-supported clades ( Figure 5D ).

At the species, Eastern, and Western distribution levels, the Mantel tests revealed positive and significant association between shape of all flight-related traits and geographic distances ( Table 2 ).

At the species distribution level, forewing size showed no significant association with geographic distances. However, the Mantel test showed a positive and significant association for the stiff portion, membranous portion, and head ( Table 2 ). For the Eastern lineage, no significant association was found between head size and geographic distances. However, a positive and significant association was observed for forewing, stiff, and membranous portions. In contrast, for the Western lineage, forewing size did not show significant results. Nonetheless, the Mantel test revealed positive and significant associations for the stiff portions and membranous portions and head ( Table 2 ).

At the species distribution level, the results of the PLS analyses for all flight-related traits revealed significant correlations with the PLS1, which accounted for, in all cases, 100% of the covariation ( Supplementary Material 5 ). For the forewing size component, PLS1 was primarily associated with relative humidity (r = −0.53) and wind speed (r = 0.44) ( Figures 6D, E , Supplementary Material 5 ). For the membranous portion size, PLS1 was predominantly correlated with relative humidity (r = −0.51) and wind speed (r = 0.41), while for the stiff portion size, PLS1 was mainly correlated with relative humidity (r = −0.44), wind speed (r = 0.39), and annual rainfall (r = −0.38) ( Figures 6C–E , Supplementary Material 5 ). For head size, PLS1 was mainly correlated with latitude (r = −0.47) and wind speed (r = 0.40) ( Figure 6E , Supplementary Material 5 ). At the Eastern lineage populations, PLS analyses for the size of all flight-related traits revealed significant correlations with the PLS1, which accounted for, in all cases, 100% of the covariation ( Supplementary Material 5 ). For forewing size, PLS was primarily associated with wind direction (r = 0.41) ( Figure 6F ) and the minimum temperature of the coldest month (r = −0.39) ( Figure 6B ). For the size of the membranous portion, PLS1 was primarily correlated with the minimum temperature of the coldest month (r = −0.41) ( Figure 5B ) and wind direction (r = 0.40) ( Figure 6F ). For the size of the stiff portion, PLS1 was mainly associated with latitude (r = −0.48), annual rainfall (r = −0.40) ( Figure 6C ), and wind speed (r = 0.37) ( Figure 6E ). For head size, PLS1 was mainly correlated with wind direction (r = 0.40) ( Figure 6F ), relative humidity (r = 0.37) ( Figure 6D ), and latitude (r = −0.37) ( Supplementary Material 5 ). For the Western lineage populations, the size of all flight-related traits revealed significant correlations with the PLS1 accounting for 100% of the covariation ( Supplementary Material 5 ). For the forewing size, PLS1 was mainly correlated with latitude (r = −0.49) and relative humidity (r = −0.43) ( Figure 6D ). The membranous portion of the forewing was mainly associated with annual rainfall (r = −0.46) ( Figure 6C ), latitude (r = −0.46), and relative humidity (r = −0.42) ( Figure 6D , Supplementary Material 5 ). For the stiff portion, PLS1 was mainly correlated with latitude (r = −0.48) and annual rainfall (r = −0.40) ( Figure 6C ). In the case of head size, PLS1 was primarily correlated with latitude (r = −0.49), longitude (r = −0.38), and minimum temperature of the coldest month (r = −0.38) ( Figure 6B , Supplementary Material 5 ).

PLS analyses for the forewing shape components across the species distribution level showed a significant correlation. PLS1 explained 79.31% of the total covariation and was mainly correlated with the NDVI (r = −0.37) and LAILV (r = 0.38) indexes ( Figure 6G , Supplementary Material 6 ). The membranous portion of the forewing shape also showed a significant correlation, with PLS1 explaining 78.86% of the total covariation and was mainly correlated with longitude (r = −0.46) and latitude (r = −0.36) ( Supplementary Material 6 ). A similar trend was observed for the stiff portion of the forewing, where PLS1 explained 61.92% of the total covariation, with significant correlations with wind speed (r = −0.44) ( Figure 6E ) and latitude (r = 0.39), among others ( Supplementary Material 6 ). For head shape, PLS analysis also indicated a significant correlation, with PLS1 explaining 77.04% of the total covariation. This PLS was mainly correlated with latitude (r = 0.50), the minimum temperature of the coldest month (r = 0.46), and annual mean temperature (r = 0.40) ( Figures 6A, B , Supplementary Material 6 ). At the Eastern lineage distribution range, forewing shape showed significant covariation with geographic, climatic variables, and vegetation indexes. PLS1 explained 87.54% of the total covariation and was mainly explained by wind direction (r = −0.39), annual mean temperature (r = −0.36), and latitude (r = 0.35) ( Figures 6A, F , Supplementary Material 6 ). The membranous portion of the forewing shape exhibited a significant correlation, with PLS1 explaining 86.18% of the total covariation, primarily associated with latitude (r = 0.40), annual mean temperature (r = −0.36) ( Figure 6A ), and altitude (r = 0.37). For the stiff portion of the forewing shape, PLS1 explained 77.01% of the total covariation and was mainly correlated with annual mean temperature (r = −0.40) ( Figure 6A ), among others ( Supplementary Material 6 ).

Head shape PLS also indicated significant correlation; PLS1 explained 75.45% of the total covariation and was mainly correlated with wind direction (r = −0.39) ( Figure 6F ) and latitude (r = 0.36). For the Western lineage distribution range, forewing shape showed significant PLS correlation; PLS1 explained 83.36% of the covariation and was mainly correlated with minimum temperature of the coldest month (r = 0.41) ( Figure 6B ), latitude (r = −0.39), and longitude (r = 0.37). The membranous and stiff portion of the forewing shape also showed significant PLS correlations. PLS1 explained 84.07% and 78.64% for membranous and stiff portions, respectively, and was mainly correlated with annual rainfall (r = −0.46) ( Figure 6C ), latitude (r = −0.46), and relative humidity (r = −0.42) ( Figure 6D ). The stiff portion was mainly correlated with wind speed (r = −0.44) ( Figure 6E ) and latitude (r = 0.41). Head shape also showed significant covariation in the PLS analysis. PLS1 retained 84.97% of the total covariation and was mainly correlated with latitude (r = −0.46) and with the minimum temperature of the coldest month (r = 0.45) ( Figure 6B ) among others ( Supplementary Material 6 ).

Results of the Procrustes regression performed with the shape component of each flight-related trait with the NDVI variation ( Figure 6G ), lineage, and population nested within the lineage showed significant effects for all components (p < 0.05). However, the correlation coefficients were low in all cases and no significant effect was found for lineage in the membranous portion of the forewing ( Table 3 ).