We classified the vegetation of GNP through the analysis of woody species composition. Woody plants represent one of the most distinguishing features of the vegetation in southern Africa (Germishuizen & Meyer 2003; O’Brien, Field & Whittaker 2000), and have been widely used for modelling vegetation structure and physiognomy in the region and elsewhere (Daru et al. 2015; Kubota et al. 2014; O’Brien 1993).

A total of 330 relevés were conducted between March and July 2010 (Figure 2), following the methodology of Timberlake, Nobanda and Mapaure (1993). Sampling was targeted to distinct vegetation communities that were clearly recognisable on satellite imagery (Landsat and Google Earth) and/or in the field. A plotless method was used for each relevé, whereby all woody species present within an area of 0.5 ha – 2 ha were recorded and classified in five height classes: seedlings, less than 0.5 m (mainly regeneration), 0.5 m – 3 m (shrubs and young trees), trees taller than 3 m (tall shrubs and trees); and mature trees (canopy height trees). Each species was allocated to a cover-abundance value in each layer using a six-point modified Braun–Blanquet scale: + = < 2%, 1 = 2% – 10%, 2 = 11% – 25%, 3 = 26% – 50%, 4 = 51% – 75% and 5 = 76% – 100%. Any species that could not be reliably identified in the field was collected and later compared with the collection available at the National Herbarium in Harare.

Environmental parameters considered relevant to the interpretation of plant distribution were also recorded. These included landscape type, lithology, soil type, slope, aspect, presence of termitaria, proportion of bare ground and litter cover, presence of surface capping, gulley and sheet erosion.

A database was created in TURBOVEG (Hennekens & Schaminee 2001). The cover values in different layers for each species were subsequently combined, in order to obtain a single value of cover per species per relevé. Calculations were performed under the assumption of random species distribution within a plot.

A Digital Elevation Model (DEM) with a spatial resolution of 90 m, representing the geographic area of interest, was downloaded from the data set of the CGIAR Consortium for Spatial Information (www.cgiar-csi.org). The DEM cell size was re-scaled to a resolution of 100 m. Eighteen bioclimatic variables were downloaded from WorldClim (www.worldclim.org) at a spatial resolution of 1 km, and subsequently downscaled to a resolution of 100 m according to the methodology of Zimmermann et al. (2007). A linear regression presenting altitude as the response variable and the bioclimatic variables as explanatory covariates was performed in STATISTICA (version 8, StatSoft Inc., Tulsa), after testing for the normality of the covariates using the Kolmogorov–Smirnov test. In case normality was rejected, an appropriate Box–Cox transformation was used. The intercept and angular coefficient values of each bioclimatic variable were then used in ArcGIS (version 10.1, Esri Inc., Redlands, California) to produce bioclimatic raster maps at a spatial resolution of 100 m. Four additional topographic variables were calculated: slope, curvature, aspect and incoming solar radiation.

Two classifications were carried out: a traditional one based on the percentage cover for each species summed across all layers and one based on the individual cover values of each species in each layer, which can be considered similar to the Integrated Synusial approach developed in Switzerland and France (Gillet & Gallandat 1996), where the classification scheme is built based on the classifications of separate vegetation layers.

In order to identify the optimal number of vegetation types for the 330 relevés, a dissimilarity matrix was created using the Kulczynsky distance method (Faith, Minchin & Belbin 1987) to calculate dissimilarities between them. For each matrix, classifications based on 2, 3, 4, and continuing up to 15 types were calculated. The clustering technique used was Partitioning Around Medoids. For each of these classifications, the Gap statistic was obtained, and the optimal number of types was set as the minimum number of types no lower than one standard error of the maximum Gap statistic observed (Tibshirani, Walther & Hastie 2001). The modified Rand Index was used in order to quantify the similarity between the two selected classifications. All statistical analyses were performed in R (www.r-project.org).

For the two classifications, for each type, diagnostic, constant and dominant species were identified in JUICE (Tichy 2002) according to the following selection criteria:

For the identification of diagnostic species, in order to remove dependence from the size of the vegetation types, a standardisation was first performed. The Fisher’s exact test was calculated, and a zero fidelity value was given to species with significance p > 0.001.

After selection of the best classification, and after verifying the normality of the environmental data, a supervised classification was performed using the observed data (i.e. the classified relevés) as a training set and the bioclimatic and topographic variables as covariates in order to produce a map of the potential natural vegetation. This procedure has its theoretical background in predictive vegetation modelling, which has been defined as ‘predicting the distribution of vegetation across the landscape based on the relationship between the spatial distribution of vegetation and environmental variables’ (Franklin 1995). Based on this assumption, several methods have been used to model the spatial distribution of vegetation types (Attorre et al. 2014; Brzeziecki, Kienast & Wildi 1993; Maggini et al. 2006). In our study, Random Forests (RF) (Breiman 2001) was used as the classifying method, and the probability of each background cell being assigned to one of the types was estimated using their climatic and topographic data. RF was chosen because of its widely recognised efficacy in such vegetation studies (Attorre et al. 2011; Bradter et al. 2011; Mellor et al. 2013; Rodriguez-Valiano et al. 2012).

The optimal type for each cell was decided through a Maximum A-Posteriori criterion. To confirm the validity of the spatialisation procedure, we used two methods: Cramer’s V index was calculated by comparing the percentage of relevés for each vegetation type expressed as a percentage of the total number of relevés with the percentage of background cells attributed to each vegetation type. Cramer’s V index is a measure of association for cross-classification tables. Additionally, we computed the out-of-bag classification error, which is an unbiased measure of misclassification directly obtained as a by-product of RF estimation. The out-of-bag-classification error is based on repeatedly sampling subsets of the available data for estimation; the data not sampled are then used as a test set for evaluating the misclassification probability.

A non-metric multidimensional scaling (NMDS) was performed in R in order to investigate the relationship between the resulting vegetation types and the environmental variables (bioclimatic, topographic and geological features).

A total of 333 woody species were identified. This included three new records for Zimbabwe (Commiphora schlechteri, Croton steenkampianus and Putterlickia verrucosa) and 12 new records for southern Zimbabwe. Fifteen of the recorded species, within Zimbabwe, are confined to the south of the country, with nine of them being recorded only from the GNP. Twenty four of the recorded species are classified as red data species within Zimbabwe, including four species classified as critically endangered (Commiphora neglecta, Periploca nigrescens, Rinorea elliptica and Schotia capitata), two endangered species (Adenium multiflorum and Milicia excelsa) and five vulnerable species (Acacia exuvialis, Adenia fruticosa subsp. simplicifolia, Manilkara concolor, Pachypodium saundersii and Suregada zanzibariensis). Thirteen alien species were recorded comprising four woody species (Lantana camara, Senna occidentalis, Senna septemtrionalis and Withania somnifera) and nine herbaceous species (Acanthospermum hispidum, Alternanthera pungens, Corchorus trilocularis, Datura stramonium, Euphorbia heterophylla, Ocimum americanum var. americanum, Sida cordifolia, Sorghum halepense and Xanthium strumarium). One or more exotic species were recorded from 28 out of the 330 samples (8%), with an overall total of 31 observations.

For both classifications, based on separate or merged cover values, the results of the Gap statistic indicated that 14 was the optimum number of vegetation types. In both cases, the resulting 14 vegetation types were allocated to three major groups corresponding to biogeographical regions: mopane, miombo and azonal alluvial woodlands. While the differences in classification were relatively minor, as confirmed by the relatively high value obtained for the modified Rand Index (58%), through comparing diagnostic, constant and dominant species, and in conjunction with expert knowledge, the classification using ‘separate layers’ was considered more effective and representative than the ‘merged layers’ classification. For instance, in the former classification, only one type was considered intermediate or ecotonal between mopane and miombo woodlands (Type 12), while for the latter classification two types were ecotonal between mopane and miombo woodlands: one comprised a rather fuzzy type of Spirostachys africana woodland and the other had similarities with the surrounding Combretaceae mixed woodlands.

Three major groups were identified consisting of miombo vegetation (covering approximately 46.97% of the GNP), mopane vegetation (38.19%) and alluvial woodlands (6.67%). Approximately 8.18% of the relevés (Type 12) could not be unequivocally assigned to any of these three main groups and have been classified as ‘intermediate’.

The vegetation types identified through the ‘separate layers’ classification, together with the number of relevés for each type, are shown in Table 1. Although there was considerable variation in the number of relevés per type (range = 12–40, Table 1), none of the vegetation types were represented by either a very small or excessively large proportion of relevés (range = 3.6% – 12.1%, Table 2). Descriptions of the 14 vegetation types based on their dominant, constant and diagnostic species, and a synoptic table with the frequencies of all woody species in each of the 14 types, are presented in Online Appendixes 1 and 2, respectively.

In the miombo woodlands group (Types 2, 5, 6, 9, 10 and 14), dominant species included Brachystegia tamarindoides subsp. torrei, Combretum celastroides subsp. celastroides, Combretum collinum subsp. collinum, Guibourtia conjugata, Julbernardia globiflora, Millettia usaramensis subsp. australis, Pteleopsis myrtifolia and Terminalia sericea. Constant species included Xeroderris stuhlmannii (five types), Cassia abbreviata subsp. beareana, Combretum apiculatum, Hugonia orientalis and Strychnos madagascariensis (each four types), and Combretum collinum subsp. collinum, Guibourtia conjugata, Monodora junodii v. junodii, Pteleopsis myrtifolia and Senna petersiana (each three types).

The dominant species in the mopane woodlands group (Types 1, 4, 7, 8, 11 and 13) were Colophospermum mopane and Combretum apiculatum, followed by Acacia exuvialis, Combretum zeyheri and Mimusops zeyheri. Colophospermum mopane and Combretum apiculatum were constant throughout all seven mopane types.

Type 12 was somewhat intermediary in nature: Combretum zeyheri achieved dominance in 4% of the samples, and the presence of Colophospermum mopane and Combretum apiculatum was 70% and 67%, respectively, which was generally lower than for the other mopane types. While the group included a number of typical mopane constant species (Cassia abbreviata subsp. beareana, Cissus cornifolia, etc.), it also shared a number of constant species that were characteristic of the miombo woodlands (Strychnos madagascariensis and Xeroderris stuhlmannii, etc.).

The alluvial woodlands were quite distinctive. The only dominant species was Cordyla africana, but there were a high number of constant species: Acacia tortilis subsp. heteracantha, Berchemia discolor, Boscia mossambicensis, et cetera. Most of these were unique to the alluvial woodlands, other than Combretum imberbe, Combretum mossambicense, Diospyros loureiriana subsp. loureiriana, Drypetes mossambicensis, Flueggea virosa subsp. virosa and Philenoptera violacea, which were shared predominantly with mopane types, but also a few species with the miombo types (Combretum mossambicense and Drypetes mossambicensis).

The results of the NMDS are shown in Figure 3a for the 14 vegetation types and Figure 3b for the three macro-groups. The correlation between the climatic gradient and the mopane and miombo macro-groups is evident (Figure 3b). The mopane macro-group is centred on the left of the graphic; these plots were generally located in ‘warmer, more arid and steeper’ areas, while the miombo macro-group is centred on the right of the graphic representing ‘cooler, moister, and higher altitude’ locations. Geology does not appear to be important for the segregation of vegetation types, except for the alluvial group.

The segregation between the 14 vegetation types is less clear (Figure 3a). Vegetation types are more readily identifiable within the miombo macro-group (e.g. Types 2, 6, 9, 10 and 14). Within the mopane macro-group it is difficult to identify single vegetation types except for Types 11 and 13. Of note is the intermediate position of Type 12 which is located at the centre of the diagram.

Potential natural vegetation maps for the 14 types and the three biogeographical regions are presented in Figure 4. The most potentially widespread types are Types 9 and 10 of the miombo region, comprising mixed Combretaceae woodland, respectively, with Burkea africana (Type 9) or Guiburtia conjugata (Type 10), generally found on sand-loam soils on flat upland areas (Figure 4a). Type 3 mixed woodland on alluvium is mainly localised on the clay-loam soils along the Runde River. The mopane types are mainly restricted to the intervening lowlands, where they occur interspersed with one another (Figure 4b).

Assessment procedures highlighted the reliability of the map. In particular, Cramer’s index is very low (< 0.007), confirming a very small (and not significant) difference between the proportion of relevés assigned to each vegetation type based on species composition and that of cells classified in the same type based on the RF model (Table 2). This result is confirmed by the low out-of-bag prediction error values obtained for each type (range = 0.000–0.048, Table 2).

The paper was aimed at providing an updated and comprehensive analysis of the woody vegetation of the GNP, within the framework of the GLTFCA.

From a methodological point of view, we demonstrated that a classification based on layers can be as effective as a traditional one where species abundance measures are merged in one value. We speculate that in averaging some information might be lost; hence, we recommend clustering based on layers measurements, which give a more direct assessment. The analysis of the 330 relevés allowed the identification of 14 vegetation types belonging to three macro-groups: alluvial (one type), mopane (six types) and miombo (six types) woodlands and one intermediate type between the last two (Table 1, Online Appendixes 1 and 2). The first two correspond to well-known biogeographical regions and confirm the findings of Daru et al. (2015) that the GLTFCA belongs to a specific phyloregion ‘Zambesian transition zone’ at the interface of mopane or miombo. We also confirmed that the mopane and miombo, even though phylogenetically similar, differ in respect to environmental parameters, with mopane vegetation being located in the lower lying, drier and warmer parts of the GNP and miombo types in moister and cooler upland areas (Figures 3b and 4b). Specifically, within the GLTFCA, beside in the GNP, true miombo (Brachystegia or Julbernardia) vegetation has been found only in the ZNP (Stalmans & Peel 2010), which is also located in the northern part of the area, while it is absent to the south (BNP, LNP and KNP).

When comparing the results of our analysis with the only previous comprehensive vegetation classification for the GNP (Sherry 1970), which found 15 vegetation types, we recorded only 8 types in common (Table 1). Among those that did not match well are, for instance, Sherry’s Spirostachys africana-Terminalia prunioides woodland that can be generically associated with the mopane macro-group, but not to any one of our types, as Spirostachys africana had significant presence in Types 4 and 12, and Terminalia prunioides in Types 11 and 13 (Online Appendix 2). Moreover, Sherry’s Kirkia, Commiphora, Adansonia open woodland occupies an intermediate position; Kirkia acuminata had a significant presence in Type 6 (miombo macro-group), Commiphora mollis was predominantly recorded from relevés classified under the mopane macro-group but was not significant for any type, and Adansonia digitata was significant for both Type 6 and for Type 3 mixed woodland on alluvium. Similarly, for Sherry’s Acacia nigrescens-Acacia welwitschii tree savanna, Acacia nigrescens was significant for Type 1 while Acacia welwitschii was not a significant component for any of the 14 types here and was mainly recorded from relevés classified under the miombo macro-group.

Other types of Sherry (1970) are difficult to match with any types from the present study. None of their characteristic species, Combretum fragrans, Terminalia stenostachya, Baphia obovata and Millettia stuhlmannii, were identified as constant, dominant or diagnostic species in any of the 14 vegetation types described here.

These differences could be because of differences in terms of data collection and classification methods. Additionally, they may relate to vegetation changes that have occurred during the past 40 years, for example, as described by Tafangenyasha (1997, 2001), Gandiwa and Kativu (2009) and Gandiwa et al. (2011a, 2011b), who investigated the decline of several tree species in the GNP.

Methodological sampling differences also influenced the comparison between the current GNP classification and those of the other parks of the GLTFCA. However, based on our results, we were able to identify several vegetation types in common.

Colophospermum mopane is the dominant species in three types that differentiate structurally as woodland (Type 7), open woodland (Type 4) and bushland (Type 13). Together with similar vegetation types recognised for the three Mozambican Parks (BNP, ZNP and LNP), they can all be compared to the wide and diversified Cissus cornifolia–Colophospermum mopane community, which was described with a large number of relevés mainly for the KNP in the revision of mopaneveld vegetation by Siebert, Bredenkamp and Siebert. (2003). In the GNP, the structural differences among the three types is hypothesised to be because of the combined effect of environmental factors, with Type 7 mopane woodland found mainly on clay soils and the other two types on various loamy soils, from clay to sandy loams (Online Appendix 1), and the impacts of elephants and fire (Gandiwa & Kativu 2009).

Colophospermum mopane–Acacia nigrescens woodland (Type 1) is mainly found on clay soils and similar communities have been found in the LNP, as Colophospermum mopane–Combretum imberbe woodland (Stalmans et al. 2004).

Colophospermum mopane–Combretum apiculatum woodland (Type 8) occupies hills and ridges similarly to the vegetation found on the central-western part of the KNP (Gertenbach 1983) as well as in protected area adjacent to its western boundaries (Peel, Kruger & MacFadyen 2007) and on the rugged veld of the LNP (Stalmans et al. 2004).

The last vegetation type belonging to the mopane group seems not to have any corresponding type in the other GLTFCA parks. Type 11 (Colophospermum mopane woodland with Combretum hereroense and Terminalia prunioides) occupies the drier areas of the Park on shallow soils and the only similar vegetation can be the Terminalia prunioides–Grewia bicolor thicket of the LNP.

Type 12 (mixed woodland with Colophospermum mopane, Combretum apiculatum and Combretum zeyheri) is an intermediate type between mopane and miombo as is indicated by its central position in the NMDS diagram of Figure 3a. In particular, Type 12 shares species from the two macro-groups: Cassia abbreviata subsp. beareana, Cissus cornifolia, Combretum mossambicense, Dichrostachys cinerea subsp. africana, Flueggea virosa subsp. virosa, Maerua parvifolia and Markhamia zanzibarica from the mopane group, and Strychnos madagascariensis, Combretum zeyheri, Lannea schweinfurthii v. stuhlmannii, Pseudolachnostylis maprouneifolia and Xeroderris stuhlmannii from the miombo group. However, considering the significant frequency of Colophospermum mopane, we decided to include this type in the corresponding group when elaborating the potential natural vegetation maps (Figure 4). Similar mixed communities are recorded from the LNP (Stalmans et al. 2004) and the BNP (Stalmans & Wishart 2005).

Mopane vegetation types identified in the GNP are very similar in term of species composition, as can be observed in the NMDS diagram (Figure 3a) where several types are strongly intermixed (Types 1, 4, 7 and 8), and in the map of the potential natural vegetation (Figure 4a), which is characterised by a marked salt-and-pepper effect. In this last case, the use of a soil map, still lacking for the entire GNP, would have increased the predictive capacity of the RF model.

Of the six miombo types identified in the present study, four correspond to types recognised by Sherry (1970), one of which is Brachystegia tamarindoides woodland (Table 1). The continued presence of this type is noteworthy as Brachystegia tamarindoides is one of the tree species considered to be in marked decline, according to Tafangenyasha (1997, 2001) being almost completely eliminated from some parts of the Park.

Other common types with Sherry (1970) include Guibourtia conjugata woodland (Type 1) and dry deciduous sandveld woodland and scrub, which corresponds to Types 9 and10 of the present study (mixed Combretaceae woodland with Burkea africana or Guibourtia conjugata, respectively). The latter two types are the most widespread ones in the GNP, particularly in the central part of the Park on cretaceous sedimentary upland areas (Figure 4a). Guibourtia-characterised vegetation types have also been identified in the other four GLTFCA parks.

Within the miombo group, we also included two vegetation types (Types 5 and 6) characterised by the presence of Androstachys johnsonii. Type 5 is similar to the Androstachys johnsonii–Croton pseudopulchellus woodland that was described by Sherry (1970) and also Van Rooyen, Theron and Grobbelaar (1981) in the northern part of the KNP, and subsequently recognised also for the LNP (Stalmans et al. 2004) and BNP (Stalmans & Wishart 2005). In the case of the GNP, the significant presence of miombo associated species such as Brachystegia tamarindoides subsp. torrei, Combretum celastroides subsp. celastroides, Combretum collinum subsp. collinum, Guibourtia conjugata, Hugonia orientalis, Hymenocardia ulmoides, Monodora junodii var. junodii, Pseudolachnostylis maprouneifolia, Pteleopsis myrtifolia, Senna petersiana, Strophanthus kombe, Strychnos madagascariensis, Terminalia sericea, Vitex mombassae and Xeroderris stuhlmannii, allowed its inclusion in this group, together with Type 6, which is very similar, as shown in Figure 3a.

Vegetation on alluvium (Type 3) was principally found along the Runde River in the north and the Mwenezi River in the south (Figure 4a). It occurred in fringing ribbons or pockets along water courses, on alluvial terraces and on flood plains. At the Save-Runde junction, it constitutes one of the largest expanses of this vegetation type in Zimbabwe. Although riparian woodlands occur widely throughout the country and adjacent regions, the portion around the Save-Runde junction comprises one of the largest remaining and relatively intact portions and is already formally recognised as comprising one of 20 Important Bird Areas in Zimbabwe, on the basis of regularly supporting a significant number of range-restricted bird species (Childes & Mundy 1997). Although alluvial vegetation was found in all the GLTFCA parks, the most similar one is that from ZNP along the banks of the Save River, with many species in common (e.g. Cordyla africana, Philenoptera violacea and Trichilia emetica).

Direct human disturbances to the vegetation have resulted through clearing for agriculture and bush clearing for the control of tsetse fly. Prior to declaration as a park, resident communities cultivated scattered fields on alluvial deposits along parts of the lower Runde River and in the vicinity of the Save-Runde junction. More recently, there has been extensive settlement and cultivation within the northern basalt plain in association with the Gulugi drainage (Dunham et al. 2010). Moreover, the bulk of the vegetation throughout the Park is highly degraded. There has been a massive reduction of trees, with upper canopy trees, in particular, having been virtually eliminated or greatly reduced over large areas, as documented in specific studies focused on Brachystegia tamarindoides (Tafangenyasha 1997, 2001), Colophospermum mopane and Combretum apiculatum in the northern part of the Park (Gandiwa & Kativu 2009), and Androstachys johnsonii and Acacia tortilis (Gandiwa et al. 2011a, 2011b). Many of the remaining large trees have serious bark damage and/or have been reduced to the level of the sub-canopy or even the shrub layer, while many of the younger trees are multi-trunked from having been knocked or burned down to ground level, probably repeatedly, and then resprouted. In many stands, relatively high numbers of dead trees were observed. Despite high levels of degradation of vegetation resources, the physical environment has remained relatively intact, with only localised occurrences of sheet and gully erosion, presumably because most of the terrain is virtually flat or gently undulating.

As for the neighbouring KNP (Eckhardt, Van Wilgen & Biggs 2000), the cause of the tree reduction appears to be because of elephants acting in concert with fire (Gandiwa & Kativu 2009; Magadza, Coulson & Tafangenyasha 1993; Tafangenyasha 2001). During the present survey, evidence of elephant damage to woody vegetation and of burning was widespread throughout the Park. Elephant densities have been sustained at relatively high levels since about 1970, although prior to 1992/1993, elephant populations were kept in check largely through culling. Since then, there has been steady growth, to the extent that the total elephant population is now higher than ever recorded, with the overall density approaching two animals per km², according to results of the most recent aerial survey carried out in September 2009 (Dunham et al. 2010).

All the four woody alien plant species, recorded during the field campaign, occur in the Type 3 mixed woodland on alluvium (Online Appendix 2). This confirms the findings of Foxcroft and Richardson (2003), who for the KNP indicated riverbanks and floodplains as ideal habitats for the establishment of invasive plants, because of disturbances caused by floods and the competitive capacity of alien species being able to germinate and to establish quicker than indigenous species after floods and to produce allelopathic agents to inhibit the growth of other species.

Beside alien species, the principal river systems of the GNP have been seriously impacted through upstream dam construction leading to reduced flow regimes and upstream agricultural activities leading to marked siltation and water pollution in the form of agricultural chemicals.