2. Materials and methods
2.1 Study sites
The study was carried out in Zambia, between the longitudes 22 and 34 °E, and latitudes 8 and 18 °S. Glossina m. centralis tsetse was collected from four sites, namely Mumbwa South (KNP1) and Kasongo Busanga (KNP2) game management areas, and Kasanka (KSP) and Sumbu (SNP) national parks (Fig 1), while G. m. morsitans was captured in five sites: Mulangu (CMR and VNP) and Luano (LVA) game management areas, and South Luangwa (SLP) and Lower Zambezi (LZP) national parks (Fig 1). The habitat of G. m. centralis is characterised by Miombo woodland interspaced with large dambos (grassy wetlands) with high annual rainfall (above 1000 mm) (Wigg, 1949). Mopane woodland is the dominant vegetation in the G. m. morsitans range with moderate to low annual rainfall (less than 800 mm).
Fig 1. Distribution of G. m. centralis and G. m. morsitans in Zambia . Data on each subspecies in Muyobela et al. (2023). The base map layer was obtained from the Database of Global Administrative Area GADM (https://geodata.ucdavis.edu/gadm/gadm4.1/shp/gadm41_ZMB_shp.zip) and under the license https://gadm.org/license.html. The figure was created using QGISv3.0 (http://qgis.org/en/site/).
2.2 Tsetse samples
The data used in this study form a subset of results of a cross-sectional tsetse survey conducted between September 2021 and August 2022 (Muyobela et al., 2023). The subset consists of flies captured in November 2021, and were chosen because this was the only month that recorded catches in all sample sites. The sampling was done using the vehicle-mounted sticky trap (VST) (Muyobela et al., 2021) baited with butanone and 1-octen-3-ol dispensed at a rate of 150 and 0.5 mg/hr. respectively (Torr et al., 1997). Tsetse captured within a two-kilometre radius of a sample site were amalgamated from which 80 non-teneral (40 males and 40 females) flies with intact wings were selected. A total of 720 (360 G. m. centralis and G. m. morsitans ) were used in the study.
2.3 Wing measurements
The right-wing of each fly was mounted on a glass slide and affixed with transparent sticky tape. The wings were then photographed using a Leica M 165C stereomicroscope attached to a Leica camera (DMC-2900) (Leica Microsystems, Germany). The images were compiled using tpsUtil v1.79 (Rohlf, 2015) and digitised with tpsDig2 v2.32 (Rohlf, 2015). Twelve homologous landmarks defined as junctions of wing veins were identified and digitised (Fig 2A). To avoid individual bias, landmark digitisation was undertaken by the same person. To avoid operational bias during digitization, specimens were selected at random during.
Fig 2. Landmark digitisation and general Procrustes analysis.(A) Image of the 12 landmarks and the order of landmark collection from the right wing of G. morsitans . (B) Scatter plot with wireframe links of landmark configurations of all 720 wings in the dataset after Procrustes superimposition. For each landmark, the white circle indicates the location of the landmark for the average shape and the grey dots indicate the locations for individual wings.
2.4 Environmental data and processing
Elevation, annual temperature, isothermality, annual precipitation, land surface temperature, and vegetation cover are among the most important variables affecting the biology of Glossina spp (Challier, 1982; Muyobela et al., 2023; Nnko et al., 2021). Therefore, these variables were selected to assess the environmental heterogeneity of sample sites and to estimate their effect on phenotypic variation. Annual temperature, isothermality, and annual precipitation data were obtained from WorldClim Global Climate Database version 2.1 (Fick and Hijmans, 2017). Moderate Resolution Imaging Spectroradiometer (MODIS) composite time series land surface temperature day (LST) (MOD11A1) (Didan, 2015) and percent tree cover based on the Vegetation Continuous Fields (VCF) (MOD44B) (DiMiceli et al., 2015) were obtained from NASA’s EOSDIS Land Processes Distributed Active Archive Center (AppEEARS Team, 2022). Elevation data was obtained as Global 30 Arc-Second Elevation (GTOPO30) from the Earth Resources Observation and Science Center (Earth Resources Observation and Science Center/U.S. Geological Survey/U.S. Department of the Interior, 1997).
Harmonic regression was performed on monthly time series LST data using the TSA package (Kung-Sik and Ripley, 2012) in R (R Core Development Team, 2015). The first coefficient in the regression, representing the mean of the variable, was selected for further analysis. Data values for all environmental variables at each sampling site were extracted using the Raster package (Hijmans and van Etten, 2012) in R.
2.5 Data analysis
All statistical analyses were done in R using the R stats version 4.4 and Geomorph version 4.0 (Baken et al., 2021) packages. Shapiro-Wilk normality test showed that environmental data values were not normally distributed (P < 0.0005 for all variables). Therefore, linear permutation models with 2000 iterations were used to test for the difference in elevation, annual temperature, isothermality, annual precipitation, land surface temperature, and percent tree cover betweenG. m. centralis and G. m. morsitans sample sites. Principal components analysis (PCA) was used for the multivariate analysis of these environmental variables to identify the most important variables accounting for environmental variability between sample sites.
Procrustes superimposition of landmark configurations was performed using general Procrustes analysis (GPA). The procedure translated all landmark configurations to a common location, scaled them to unit centroid size, and rotated them into an optimal least-squares alignment with an iteratively estimated mean reference form (Zelditch et al., 2004) so that the sum of squared distances between corresponding landmarks of each configuration and the mean configuration was minimized (Klingenberg, 2013). This analysis produced the Procrustes distances which measure shape dissimilarity as well as the centroid size (CS). A scatter plot of superimposed landmarks for all specimens is shown in Fig 2B.
General Procrustes analysis places landmark configurations in non-Euclidean Kendall’s Shape Space, making statistical methods predicated on Euclidean relationships unusable (Webster and Sheets, 2010). Consequently, most shape analysis studies involve projecting data from this space into a Euclidean space tangent to shape space. Projecting to tangent space, however, distorts the distance between each target configuration and the reference form relative to the Procrustes distance between them in shape space. To minimize this distortion, the standard procedure in GPA is to use the mean (consensus) form of all configurations in the sample as a reference (Zelditch et al., 2004). An alternative approach to reducing the difficulties associated with the non-linearity of shape space is to use statistical procedures based on bootstrap or permutation methods (Webster and Sheets, 2010). These methods not only avoid the problem of distortion but also allow for the use of a non-central shape as a reference thus facilitating the comparison of sample subset shapes. Therefore, permutation procedures as described by Adams and Otárola-Castillo (2013) were used for shape analysis. Shapiro-Wilk normality test showed that both CS and log CS were not normally distributed (P < 0.0005 for both variables). Therefore, permutation procedures were also used to analyse CS.
Linear permutation models with 2000 iterations were used to compare wing CS differences betweenG. morsitans males and females, G. m. centralis, andG. m. morsitans subspecies as well as CS differences between sample locations within each subspecies range. The pairwise function in R was used for multiple group comparisons. Linear permutation models were also used to estimate the effect of elevation, annual temperature, isothermality, annual precipitation, land surface temperature, and percent tree cover on wing CS.
Hypothesis testing in multivariate linear permutation models of shape was accomplished using Goodall’s F-test (Goodall, 1991), a statistical approach that partitions the variance of Procrustes distances rather than landmark coordinates. Goodall’s F-statistic is the ratio of explained (between-group) and unexplained (within-group) components of shape variation (Klingenberg, 2016) and has been demonstrated to have higher statistical power than other approaches (Rohlf, 2000). To test whether there was significant covariance between wing shape and size (allometry), multivariate linear permutation regression of wing shape on CS was conducted. Residues from this regression were then used to construct allometry-free shape variables that are recommended in taxonomic investigations (Klingenberg, 2009) and studies that define geographically constrained situations such as islands (Dujardin, 2011). The multivariate regression approach to remove the allometric component of shape variation offers a logical method as it partitions the variation in the dependent variables into predicted and residual components (Klingenberg, 2016). The predicted component corresponds to allometric variation of shape, whereas the residual component encompasses non-allometric variation as residues are uncorrelated with CS. Multivariate linear permutation models were used to compare allometry-free wing shape differences between G. morsitans males and females, G. m. centralis, and G. m. morsitans subspecies as well as differences between sample locations. Where differences between locations were observed, multiple group comparison was done using the pairwise function. Thin-plate spline deformation grids and vector wireframe plots were generated to illustrate the observed relative shape change. To estimate the amount of shape variation that could be attributed to environmental variability, size-adjusted shape was regressed on elevation, annual temperature, isothermality, annual precipitation, land surface temperature and percent tree cover. A Procrustes distance matrix was used to build a neighbour joining tree to illustrate divergence of wing shape among flies from different locations. Alpha was set at 0.05 for all statistically significant analysis (Pirk et al., 2013).
2.6 Ethical Statement
The protocol and procedures employed in this study were reviewed and approved by the Department of Zoology and Entomology at the University of Pretoria.