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.