Livestock grazing and biodiversity: Effects on CO 2 exchange in semi-arid Karoo ecosystems,

2 ). We present four years of measurements from twinned eddy-covariance towers in Nama-Karoo, South Africa, to investigate the carbon fluxes and the impact of grazing intensity on NEE. The design contrasted NEE at a long-term site grazed at recommended levels (LG) with a long-term heavily grazed (EG) site that had been rested for 10 years, and was monitored for two years after which intensive grazing was reintroduced for this experiment. This allowed for the quantification of long-term NEE trends on “ recovering ” vegetations (years I, II) and short-term responses to an intensified land use (years III, IV). The results showed that the net release of CO 2 was slightly higher at LG than on “ recovering ” vegetation at the EG site, where near-neutral exchange was observed during years I and II. However, after grazing was reintroduced to the EG site, differences between sites was reduced but not eliminated. These findings suggest that there

• Effects of grazing on CO 2 exchange were studied at two sites with different management but similar climatic conditions.• The ecosystem that was less diverse and overgrazed in the past showed slightly higher carbon sequestration rates.• The ecosystems can quickly switch from a net carbon source to sink depending on water availability and distribution.).We present four years of measurements from twinned eddycovariance towers in Nama-Karoo, South Africa, to investigate the carbon fluxes and the impact of grazing intensity on NEE.The design contrasted NEE at a long-term site grazed at recommended levels (LG) with a longterm heavily grazed (EG) site that had been rested for 10 years, and was monitored for two years after which intensive grazing was reintroduced for this experiment.This allowed for the quantification of long-term NEE trends on "recovering" vegetations (years I, II) and short-term responses to an intensified land use (years III, IV).
The results showed that the net release of CO 2 was slightly higher at LG than on "recovering" vegetation at the EG site, where near-neutral exchange was observed during years I and II.However, after grazing was reintroduced to the EG site, differences between sites was reduced but not eliminated.These findings suggest that there

Introduction
While terrestrial ecosystems currently absorb roughly 30 % of annual anthropogenic CO 2 , the location and long-term sustainability of these carbon sinks is poorly known (Ciais et al., 2014).Diverse methodologies provide a range of estimates of the sizes of these sinks, but lack of measurements to permit adequate spatial modeling and credible extrapolation (Friedlingstein et al., 2020).In particular, information is needed in southern hemisphere countries to evaluate the magnitude of terrestrial carbon sinks and their vulnerability to climate change and local drivers such as land degradation.In the semi-arid subtropics, extensive rangeland systems represent a vast area that has shown remote-sensed evidence of greening, and potential net sink activity for some decades (Ruijsch et al., 2023).These systems play an important role in global carbon dynamics: by partitioning terrestrial CO 2 among land cover classes, Ahlström et al. (2015) showed that the interannual variability of the net biome production is dominated (39 %) by semi-arid ecosystems.
Due to the global warming in the last decades, assessment of net carbon balance has become a critical issue in environmental science (Ardö et al., 2008).Even though numerous studies on ecosystem carbon budget have been conducted, most of them have concentrated on European, Asian, or North American mid-latitude ecosystems (Anthoni et al., 2004;Bao et al., 2019;Delgado-Balbuena et al., 2019), while a lack of attention has been given to the African continent.
Previous studies determined water availability as the main driver controlling gross primary production (GPP) through increasing photosynthesis and extending the length of the growing season in global semiarid systems (Kutsch et al., 2008).South African semi-arid ecosystems are highly seasonal with large parts of total annual precipitation occurring in only a few months of the year, and generally reveal pulsedriven evapotranspiration events in the rainy season.Phase and amplitude of photosynthesis and respiration are also highly dependent on the precipitation pattern (Kutsch et al., 2008;Merbold et al., 2009).The strength of the ecosystem response to rainfall is not only determined by the amount of an individual event, but also on the timing of preceding events (Huxman et al., 2004;Ivans et al., 2006).
Land use change is one of the key drivers affecting species richness, vegetation and ecosystem structure, thus altering ecosystem-atmosphere carbon exchange (Canadell, 2002).In semi-arid regions, browsing and grazing of rangelands is a dominant land use, but shifts in this use are subject to national policies and economic drivers.While many studies suggest that overgrazing leads to biodiversity losses and soil degradation (Kairis et al., 2015), there are some studies in southern African shrubland rangelands that reveal resilience through shifts in biodiversity rather than mere declines (Rutherford and Powrie, 2013).
With vast areas of the semi-arid sub-tropics subject to differential patterns of intensifying and reducing browsing pressure, quantifying impacts on carbon exchange is an important component of constraining impacts on the global carbon budget, and informing local policy makers of the implications for their mitigation plans and reporting needs.This region has seen substantive shifts in land use pressure over the past decades, from a period of high grazing pressure, to a period post when stock numbers have been reduced country-wide.Reyers et al. (2009) reported that about 52 % of the semi-arid Karoo has been moderately or severely degraded due to historical livestock overgrazing, but Hoffman et al. (2018) showed that the situation has shifted through the reduction in livestock numbers (for example, the number of sheep has reduced from 11 million (1939) to 4 million ( 2007)), and by extension of protected areas (from 0.03 % to 1.6 %).Most studies providing comparisons of carbon exchange under different grazing intensities are of short duration, hence they do not detect seasonal, nor annual changes in the carbon dynamics associated with grazing, and are not based on South African semi-arid ecosystems (Chimner and Welker, 2011;Wang et al., 2016).There are no such studies in southern African rangelands, despite their importance nationally, and potentially in international terms.This region therefore offers important opportunities for quantifying long-term changes in function due to changing land use patterns.
In addition to the impact of changing cover and structure under land use change, the importance of biodiversity in supporting productivity and resilience is increasingly recognized in the context of climate and anthropogenic change.In particular, it is suspected that biodiversity enhances ecosystem functioning (Cadotte et al., 2008;Cardinale et al., 2007), and thus may increase resilience important for establishing sustainable management practices.Numerous observational studies have reported a positive relationship between biodiversity and productivity: more diverse ecosystems are more productive and better adapt to climate change (Isbell et al., 2015;Kreyling et al., 2017).However, some experimental studies have yielded mixed results (Adler et al., 2011;Grace et al., 2007).Observational studies can provide a general picture of the relationships between biodiversity and productivity, but so far existing knowledge is derived from the relatively crude assumption that higher biodiversity can be represented simply by increasing species numbers (richness) and is based on artificial community modifications (Brun et al., 2019;Chalmandrier et al., 2017).This simplification may constrain our understanding of the relationship between biodiversity as expressed in complex natural communities, and measures of productivity that are needed to inform successful management responses.
In this study, we explored the links between NEE, vegetation cover, and composition (biodiversity) in response to historical and current livestock grazing impacts.The study is based on four years of continuous eddy covariance (EC) measurements of CO 2 at the semi-arid Karoo near Middelburg, Eastern Cape, South Africa.This dataset provides an appreciable range of rainfall conditions, from dry to wet, to illustrate the interannual climatic variability characteristic of these understudied ecosystems.Our paired-tower approach allowed for an analysis of the effects of different grazing regimes on CO 2 fluxes under identical climatic conditions.The two study sites represent lenient temporary (LG) vs. experimental grazing (EG).Our expectations, based on broad O. Rybchak et al. consensus in the literature are that 1) carbon sequestration potential would be lower in the EG site than the LG site, due to lower biodiversity of the EG site; 2) reintroducing grazing to the EG site should result in a higher amount of carbon emissions 3) carbon uptake should be increased under higher-than-average rainfall conditions.The experimental design allowed us 1) to quantify the NEE budgets and analyze the diurnal, seasonal, and interannual variability as well as differences between the sites, which is crucial for evaluating the magnitude and susceptibility of terrestrial carbon sinks to climate change and local drivers in South Africa through the use of spatial modeling; 2) to analyze the impact of reintroducing grazing on the EG site (latter two measurement years, under heavy grazing regime) to determine the resilience of the system functioning in relation to the LG site; and 3) to investigate the sensitivity of carbon fluxes to rainfall variability.

Site description
The two study sites were located at an altitude of 1310 m.a.s.l., in the Eastern Upper Karoo vegetation type in the Nama-Karoo, a biome that occupies much of the central to the western inland region of South Africa (Mucina et al., 2006) (Fig. 1).The vegetation is co-dominated by dwarf shrubs (perennial, both succulent, and non-succulent) and grasses (short-lived and perennial), with shrubs, geophytes, sedges, and herbs also present.In the case of high precipitation, grasses may overcome shrubs, while under conditions of low precipitation, grasses and dwarf shrubs are of minor importance (du Toit et al., 2018;O'Connor and Roux, 1995).The grass growth is favored during the warm season, while the cool season promotes the growth of dwarf shrubs (du Toit and O'Connor, 2014;Vorster and Roux, 1983).The soils are loamy, calciumrich Solonchaks overlying solid rock at both sides (Roux, 1993).Four basic seasons can be distinguished throughout the year: cold and dry winter (June-August), warm and relatively dry spring (September-November), hot and wet summer (December-February), and warm and relatively wet autumn (March-May) (du Toit and O'Connor, 2014).During the summer months, days are generally hot (30-40 • C) and nights are moderately warm (10-16 • C), while the winter days are moderate to warm (14-25 • C) and nights are cold (− 4-4 • C).The longterm mean annual temperature was 15 • C. Precipitation (annual variation up to 40 %) and droughts are unpredictable and changeable (Booysen and Rowswell, 1983).Mean annual precipitation was 374 mm from 1889 to 2013 in the range from 163 mm to 749 mm (du Toit and O'Connor, 2014).Precipitation mainly occurs in the mid-summer to midautumn months with March as the wettest month (Venter and Mebrhatu, 2005).The seasonality and amount of precipitation, including droughts and prolonged wet periods, is the main driver of semi-arid ecosystem processes, especially for vegetation dynamics, composition, structure, and functioning (Gentry, 1988;du Toit and O'Connor, 2020).Due to the long dry season, the vegetation needs more time to start growing and is not obviously in phase with seasonal air temperature (T air ) and shortwave incoming radiation (R g ) development but follows the SWC trends (Grossiord et al., 2017;Mucina et al., 2006).Drought is a common phenomenon in the Karoo and severe droughts have occurred with a frequency of about 20 years (du Toit, 2017).
In the 18th century, grazing intensity in the Nama Karoo biome was usually low and determined by the distribution of the surface water available to animals (Milton, 1993;du Toit and O'Connor, 2020).In the 1800s with the arrival of European farmers, much higher livestock stocking rates were introduced which had a detrimental effect on the vegetation (Archer, 2000;Shaw, 1873).Following the drought period (1900)(1901)(1902)(1903)(1904)(1905)(1906)(1907)(1908)(1909)(1910)(1911)(1912)(1913)(1914)(1915) with a mean annual precipitation of 288 mm, a commission of inquiry issued recommendations on the control of livestock grazing regimes, dividing it into paddocks (du Toit et al., 1923;du Toit and O'Connor, 2020).This division allowed an optimized distribution of livestock on the farm and provided rest for paddocks where required.Based on those recommendations and grazing trials (Tidmarsh, 1951), models of vegetation response to grazing and rainfall were developed with guidance for farmers to best manage the veld (Moll and Gubb, 1989;Roux and Vorster, 1983).The main discovery was that repeated livestock grazing in the same area in the same season can strongly affect the vegetation composition (O'Connor and Roux, 1995;Tidmarsh, 1951).
The Lenient Grazing (LG) (31 • 25′20.97″S, 25 • 1′46.38″E)site has been grazed by sheep and cattle (Fig. A1) using a rotational grazing system (approximately 2 weeks grazing followed by 24-26 weeks rest) at recommended stocking rates of 1/16 animal unit per hectare (ALU ha − 1 ) since the 1970s, and is considered to be an excellent condition 'benchmark' site in terms of botanical composition, with a wide diversity of species and co-dominance of grasses and dwarf-shrubs (du Toit and O. Rybchak et al. Nengwenani, 2016).Dominant species are Digitaria eriantha (palatable perennial grass), Pentzia globosa (palatable non-succulent dwarf-shrub), Eriocephalus ericoides (palatable shrub), and Sporobolus fimbriatus (palatable perennial grass) (Table A1) (du Toit and Nengwenani, 2019).The Experimental Grazing (EG) (31 • 25′48.69″S, 25 • 0′57.70″E)site was grazed by Dorper sheep using a 2-paddock rotational system (120 days grazing followed by 120 days rest) at stocking rates approximately double the recommended rate (2/16 ALU ha − 1 ) as part of an experimental trial from 1988 to 2007.The Dorper breed is described as a hardy sheep that prefers shrubs to grasses (Brand, 2000).This severe treatment extirpated nearly all palatable species and nearly all dwarf shrubs, and as a result, the system is dominated by Aristida diffusa (unpalatable perennial grass), Aristida congesta (short-lived unpalatable grass), and Tragus koelerioides (creeping unproductive grass) (van Lingen, 2018).Thus, the EG site is degraded from an agricultural point of view, having shifted from a diverse grassy shrubland to unpalatable arid grassland.The site was ungrazed 2008-2017 but did not recover, after which Dorpers were reintroduced at a somewhat higher stocking rate (1/5 ALU ha − 1 ) in July 2017, and the veld was grazed continuously unless (for short periods) food availability was too low.Vegetation biomass has almost never been absent (du Toit et al., 2018;du Toit and O'Connor, 2020).Animals were moved to the next paddock long before all vegetation was removed.Besides, non-growing vegetation retains its qualities well (almost like standing hay) and remains palatable to animals.The species richness of the LG site is higher containing 32 species compared to 10 species of the EG site (Table A1).At the LG site, the paddock is 550 × 200 m in size, while at the EG site, the paddock is 200 × 250 m with the tower in the middle (Fig. A1).Climatic conditions of the two sites are similar.There is no information on the difference between the sites in soil bulk density, nitrogen, and organic carbon content.

Instrumentation and measurements
Eddy Covariance (EC) towers were installed at the LG and EG sites to measure ecosystem-atmosphere exchange of CO 2 , water and sensible heat fluxes in October 2015.The two towers are placed approximately 1.5 km apart.Wind components were measured with a threedimensional sonic anemometer (CSAT3, Campbell Scientific Inc., Logan, UT, USA) placed 3 m above the ground.The sonic anemometer was coupled with an enclosed path fast-response Infra-Red Gas Analyzer (IRGA) Li-7200 (IRGA, Li-Cor, Lincoln, NE, USA) for CO 2 and H 2 O concentration measurements.Air from close proximity (~5 cm) to the sonic anemometer's array was pulled through a sampling tube into the IRGA through a stainless steel tube of 70 cm length and 6.0 mm inner diameter using the pump from Licor's LI-7200-101 Flow Module with a flow rate of 15 L min − 1 .The IRGA acquisition system collected raw data at 20 Hz sampling frequency and compressed files (meta, flux, and biomet data) for later post-processing.The gas analyzers were manually calibrated every six months using standard air calibration tanks with mixing ratios of zero and span (431.23 ppm).Drift of the sensors was always below 1 %.
Relative humidity and air temperature were recorded with Temperature/Humidity Probe (HMP155, Vaisala, Helsinki, Finland), precipitation with Tipping Bucket Rain Gauge (TR 525, Texas Electronics, Texas, USA), and radiation components with Net Radiometer (CNR4, Kipp&Zonen, Delft, Netherlands).Soil measurements consisted of heat flux plates (HFP01 + HFP01 SC, Campbell scientific, Utah, USA) at two different depths (10 and 20 cm) with three repetitions, soil temperature probes (UMS TH3 sdi12, Meter AG, Munich, Germany) at six different depths (5,6,10,13,20 and 37 cm at LG site and at 10,12,17,22,32, and 57 cm at EG site), soil moisture probe (ML3x Delta T, EcoTech, Bonn, Germany) at two different depths (10 and 16 cm at the LG site and at 8 and 14 cm at the EG site).The two sites differed in soil heterogeneity and characteristics particularly in stone content (soil skeleton >2 mm for the EG site), which were the reasons for different installation depths of the soil temperature and moisture probes.

Data processing 2.3.1. Eddy covariance data processing
Fluxes were calculated on the basis of the EC technique from raw high-frequency data with a half-hourly averaging period, discussed, and explained in detail in previous studies (Aubinet et al., 2000;Baldocchi, 2003;Burba, 2013;Falge et al., 2001;Foken and Wichura, 1996;Moncrieff et al., 1997;Webb et al., 1980) strictly following the recently published processing chain of the FLUXNET methodology by Pastorello et al. (2020).Processing was performed with the EddyPro software package version 7.0.6 to quantify ecosystem-atmosphere exchange of CO 2 , water vapor, and sensible heat fluxes.Flux data were available for the four hydrological years from November 2015 to October 2019.Complete datasets consist of 70,128 30-min flux measurements with 8 % for LG site and 12 % for EG site missing data due to instrument failure.After processing, the gaps increased to 12 % and 26 % of the missing data for LG and EG sites, correspondingly.
The absolute limits assessment was performed to filter out implausible data by setting upper and lower thresholds (Aubinet et al., 2000;Baldocchi, 2003;Foken and Wichura, 1996).Spikes were removed by applying the median absolute deviation (MAD) filter (Mauder et al., 2013).The amplitude resolution test was performed according to Vickers and Mahrt (1997).The vertical and horizontal velocity components were transformed to zero to adjust wind statistics for the misalignment of the sonic anemometer by performing a twodimensional coordinate system rotation (Aubinet et al., 2000;Baldocchi et al., 1988;Wilczak et al., 2001).The time lag compensation between the anemometric and gas analyzer variables was detected using the maximum covariance method (Fan et al., 1990).Linear detrending was used to extract turbulent fluctuations from time-series data (Gash and Culf, 1996).Also, fluxes were corrected following systematic losses in the low and high-frequency range due to the data processing choices, system characteristics, and sampling ranges of the instruments (Moncrieff et al., 1997(Moncrieff et al., , 2006)).Quality control flags were calculated for all fluxes by the flagging methodology after Foken et al. (2006) and Mauder and Foken (2011).After application of all the correction procedures, half-hourly values of CO 2 , latent and sensible heat fluxes were calculated.

Flux filtering and quality control
Based on the criteria after Mauder and Foken (2011), bad quality flux data that obtained the value "2" were discarded from further analyses.In total, after the quality checking and quality flags filtering procedures, the missing CO 2 flux measurements data increased to 27 % for the LG site and 34 % for the EG site.
Numerous studies indicate reduced reliability of the nighttime flux measurements due to periods with low turbulent mixing of the planetary boundary layer (Aubinet et al., 2000;Gu et al., 2005;Twine et al., 2000).Many authors describe a threshold value for the friction velocity (u * ) filtering criterion to remove periods of intermittent turbulence as the EC assumptions fail during low turbulence periods (Aubinet et al., 2000(Aubinet et al., , 2012;;Baldocchi, 2003;Goulden et al., 1996).To identify insufficient turbulent mixing, we estimated the minimum friction velocity u * according to the method described in Papale et al. (2006).Threshold values of u * were 0.10 m s − 1 for both sites.The corrected fluxes were filtered by removing data below these thresholds from the dataset.After u * filtering application, the missing CO 2 fluxes data increased to 44 % and 49 % for the LG and EG site, respectively.After the standard quality assessment and quality control, flux data sets were again checked and filtered for implausible values.

Data gap-filling and flux partitioning
Gaps were filled to obtain continuous dataset for the calculations of the net CO 2 balance and its component fluxes GPP and R eco .
CO 2 fluxes were gap-filled and GPP and R eco were derived using the gap-filling and flux partitioning REddyProcWeb online tool.The detailed explanation of the gap-filling and partitioning methods are described by Wutzler et al. (2018).The marginal distribution sampling (MDS) method after Reichstein et al. (2005) is implemented into the tool, which integrates the covariation of the fluxes with the micrometeorological conditions within a moving time window and the temporal fluxes' autocorrelation based on look-up tables (LUT) and the mean diurnal course (MDC) method (Falge et al., 2001;Wutzler et al., 2018).
The main components of the carbon fluxes (GPP and R eco ) can be derived by two methods: the so-called nighttime approach by Reichstein et al. (2005) and the daytime approach by Lasslop et al. (2010).Both approaches were tested extensively for ecosystems mainly in the midlatitudes but less for highly seasonal ecosystems.As the temperaturerespiration relationship was found to be weak in semi-arid ecosystems (Räsänen et al., 2017), we used the daytime approach from Lasslop et al. (2010) based on the light-response curves.The footprint estimation was performed according to the "simple footprint parameterization" described in Kljun et al. (2004) using TOVI software (Tovi™ Software, 2020).

Uncertainty estimation
Flux measurements by the EC method are subject to errors, and estimating them is an important criterion for analyzing annual net ecosystem exchange (NEE) balances.The annual NEE uncertainty is assumed to take into account the most significant possible error sources (Finkelstein and Sims, 2001;Lucas-Moffat et al., 2018;Massman and Lee, 2002;Moffat et al., 2007;Richardson et al., 2012).These are (1) systematic errors associated with advection, flux divergence and tilt correction, (2) random errors for example associated with inadequate sample size, and (3) bias errors due to the gap-filling of the unavoidable gaps in the EC data sets.As systematic errors are difficult to estimate without additional observations, we follow the approach of Moffat et al. (2007) and Lucas -Moffat et al. (2018) for NEE uncertainty estimation.
The bias errors resulting from gap-filling of the flux measurement data is provided by the formula: where δASum is the offset on the annual sum, Np is the number of gapfilled days and BE is the bias error.Moffat et al. (2007) demonstrate that the effect of different gap-filling techniques on the bias error of the annual balance of NEE for their selection of sites was <0.25 g C m − 2 day − 1 .The random error (RE) was calculated for each 30 minute flux measurement using the Finkelstein and Sims (2001) method.The random error of the annual sum (εASum) is calculated by the formula: Both equations were summed up to estimate total NEE uncertainty (E_NEE) for the year: (3)

Remote sensing-based NDVI index
Remote sensing-based vegetation indices (VIs) were used to put the observed flux tower measurements into the context of plant functioning and responses to environmental conditions and management.To this end, Moderate Resolution Imaging Spectroradiometer (MODIS) time series provided by the Terra satellites were acquired for the period under investigation.More specifically, the MODIS VI product MOD13Q1 (Terra) was downloaded from the National Aeronautics and Space Administration (NASA) Distributed Active Archive Centre (DAAC, https://ladsweb.modaps.eosdis.nasa.gov/search/).Subsequently, the VI values contained in these products were extracted for the two flux tower locations and compared to the EC measurements.
MOD13Q1 is a global data set featuring a spatial resolution of 250 m and a temporal resolution of 16 days.One of the VIs that comes with the product is the normalized difference vegetation index (NDVI).It allows for consistent spatio-temporal comparisons of vegetation canopy greenness and canopy structure.The most recent version (collection 6) of MOD13Q1 was employed in this study.For more information, the reader is referred to the products' user guide (Didan et al., 2015) and their algorithm theoretical basis document (ATBD) (Huete et al., 1999).

Statistical analysis
Statistical analyses were performed using the 'SciPy.Stats' package in Python (Virtanen et al., 2020).The Student's t-test was used to verify the statistical difference with a significance level set at p < 0.05.We compared: -Daily NEE fluxes between the studied sites to investigate the impact of past overgrazing (years I and II) and current heavy grazing (years III and IV), -Daily NEE fluxes between the years to investigate the influence of rainfall on the variability of CO 2 fluxes, -NDVI values between the studied sites, again separately looking at years I + II and III + IV to be able to distinguish the effects of past and current overgrazing.
Polynomial regression analysis was conducted to investigate the relationship between the response variable (NEE) and the main environmental driver (SWC).

Diversity indexes
The Shannon-Wiener diversity index is a quantitative indicator that represents the species abundance and how evenly they are distributed among the population.The Shannon-Wiener's diversity (H ′ ) was calculated by the formula (Morris et al., 2014;Shannon, 1948;Spellerberg and Fedor, 2003): where p i is the proportion of species i, snumber of species collected in a sample.The effective number of species, so-called the true diversity of the community, refers to the number of equally abundant species in a community needed for the average proportional abundance of the species to equal that observed in the dataset.The effective number of species (ENS) was calculated by the formula (MacArthur, 1965;Tuomisto, 2010bTuomisto, , 2010a)): (5)

Artificial neural networks (ANNs) analysis framework
ANNs have demonstrated remarkable effectiveness in identifying patterns within large and noisy datasets, surpassing the performance of classical semi-empirical methods.ANNs based on the backpropagation algorithm can also be regarded as multivariate nonlinear regression models (Bishop, 1995;Rojas, 1996).These can be used in an inductive approach to analyze primary and secondary drivers of ecosystem fluxes, for a detailed description see (Moffat et al., 2010;Moffat, 2012).
We applied this methodology here to identify if there are any noticeable differences in climatic drivers of the daytime NEE between the EG and LG site.For this analysis, the dataset has been filtered to include only non-gapfilled, daytime (QF = 0, R g ≥ 5 W m − 2 ).

Meteorological conditions
The measured daily and seasonal variability in the main climatic conditions (air temperature, radiation, relative humidity) were typical for the Nama Karoo biome (Fig. 2).We analyzed growing (January-May) and dry (June-December) seasons and defined the following periods as hydrological years (HY): - Mean annual temperature (T air ) during the study period (November 2015-October 2019) was 15.7 • C for both sites (Fig. 2a, Table B1).Daily mean relative humidity ranged from a minimum of 6 % in dry seasons up to 99 % during the growing seasons (Fig. 2b).Relative humidity (RH) was lowest in November and December, while the highest values were measured in February-May.Short-wave incoming radiation (R g ) showed typical seasonal fluctuations with higher values during the summer months (December-January) (Fig. 2c, Table B1).The seasonal changes in R g showed that the EG site had slightly higher values.The mean daily R g including nighttime during the four years measurement period were 240 and 250 W m − 2 for the LG and EG sites, correspondingly.The peak values reached 1000-1200 W m − 2 in summer and 400-600 W m − 2 in winter.
The LG and EG sites had mean annual precipitation of 378 mm and 365 mm, respectively.Approximately 55 % of the total annual precipitation happened during the summer months (January-February), with autumn (March-May) as the second most humid season (approx.30 %), winter (June-August) as the driest season with precipitation <10 mm per month and spring (September-November) characterized by short rain events (Fig. 2d, Table B1).Precipitation was highly variable during the measurement period.Years I and II were nearly the same as the longterm mean annual precipitation for the Karoo (376 mm and 365 mm in the LG site and 372 mm and 374 mm in the EG site, respectively).The last year of measurements was the driest year with annual precipitation rates of 295 and 267 mm for LG and EG sites, respectively, while Year III had the highest precipitation of 487 and 458 mm at LG and EG, respectively.
Seasonal variation of soil water content (SWC), expressed as volumetric water content, changed with the highest values in the summer and autumn months (growing season) and the lowest values during the winter (Fig. 2d, Table B1).The SWC values for the EG site were similar, however consistently slightly lower than those for the LG site.The highest values of SWC were in February-April 2018 (30 %) in the period with the highest precipitation, and the lowest in June-November 2019 (12 %), in the driest year over the measurement period.
The NDVI was below 0.3 during the dry season and began to increase rapidly at the onset of the growing season with the maximum values at the peak of the growing season (Fig. 2e).With the onset of the dry season, the NDVI was back at a level below 0.3.The LG site is more diverse with the Shannon index of 2.12 and 8 effective number of species compared to the EG site with the Shannon index of 1.15 and 3 effective number of species.
The wind roses showed that the prevailing wind directions at the studied sites were NNW, N, NNE, and SSE with the wind speed mainly between 1 and 6 m s − 1 (Fig. 1).The footprint indicates that around 90 % of the carbon fluxes measured by the EC method were within 70 m of the flux tower (Fig. A1).

Diurnal patterns of carbon fluxes
The mean diurnal patterns of carbon fluxes are shown in Fig. 2 for each year and lumped together for the months of the dry (June-December) and growing (January-May) seasons.During the dry seasons (June-December), both ecosystems were inactive (Fig. 3a, b, c, d).The carbon flux values were close to zero for the most part as net carbon uptake was constrained due to water unavailability and vegetation dormancy.Mean diurnal carbon exchange rates in the dry seasons were the lowest in the year I and the highest in the year II (Table C1).Despite the similar mean daytime carbon exchange rates, slightly higher carbon During the growing seasons, the diurnal cycles of carbon fluxes for both sites demonstrated a classical behavior with the net carbon uptake predominating during the daytime and only respiration occurring during the night time (Fig. 3e, f, g, h).The daytime mean carbon uptake rates were higher for the EG site in the years I and II compared to the LG site (Table C1), while no significant difference was observed in the years III and IV.Overall, the maximum daytime carbon uptake means and nighttime release were both observed in the year III (Table C1).
The ecosystems turned from a carbon source to a carbon sink mostly from 6:00 to 18:00 LT, with the highest carbon sequestration values near midday (Fig. C1).Midday uptake rates ranged from 48 mg C m − 2 h − (year I) to 206 mg C m − 2 h − 1 (year III) for LG site and from 81 mg C m − h − 1 (year I) to 219 mg C m − 2 h − 1 (year III) for the EG site.

Seasonal and annual NEE, GPP and R eco variations and its relationship with water availability
The NEE, R eco , and GPP at the studied sites varied seasonally with the peak values in the middle of the growing season (February-March) (Figs. 4, C2).Their seasonal patterns were in accordance with changes in precipitation and SWC (Fig. 4).The duration of the total number of continuous NEE < 0 days lasted 54 and 66 days (year I), 96 and 98 days (year II), 154 and 148 days (year III), and 95 and 96 days (year IV) for the LG and EG sites.Throughout the measurement period, the studied ecosystems acted as a net carbon sink for 440 and 511 days (16 months) for the LG and EG sites, respectively.
At the LG site, the magnitude of the daily NEE during the measurement period varied from − 3.2 g C m − 2 d − 1 to 4 g C m − 2 d − 1 , while at the EG site, it ranged from − 3.7 g C m − 2 d − 1 to 3.4 g C m − 2 d − 1 (Fig. 4, Table C2).Year III had the highest precipitation rates and the highest carbon uptake rates (Fig. 4).In comparison, year IV had the lowest precipitation rates and the highest carbon release rates.The difference between the sites decreased after livestock grazing was reintroduced in the EG site.Throughout the measurement period, the EG site had slightly higher carbon uptake and lower carbon release rates compared to the LG site (Figs.4, C2).
The dry periods were mainly characterized by inactive ecosystems with low precipitation (<10 mm month − 1 ), low soil water content (<15 %) and low GPP, R eco , and NEE (close to zero).Under these conditions, both sites acted with similar maximum and minimum daily NEE (Figs. 4, Table C2).Mean values of carbon fluxes throughout the measurement period for the growing seasons were − 70 mg C m − 2 day − 1 for LG site, − 150 mg C m − 2 day − 1 for EG site, while during dry seasons both sites had similar mean values (~100 mg C m − 2 day − 1 ).
In contrast to the dry season, the growing season was characterized by high precipitation (>65 % of annual precipitation) and high soil water content (15 %-35 %), resulting in highest NEE, R eco and GPP values (Fig. 4).The highest daily sums of R eco and GPP were in the year III, whereas the lowest daily sums were in the year I for both sites (Table C2).The maximum GPP during the measurement period reached 7.5 g C m − 2 day − 1 for the LG site compared to 6.3 g C m − 2 day − 1 for the EG site.The maximum values of R eco were 7.5 g C m − 2 day − 1 and 4.4 g C m − 2 day − 1 for the LG and EG sites, respectively.Increased release of CO 2 occurred during the transition period (November-December) from the dry to the growing season.GPP exhibited delayed response to precipitation (~1-3 weeks) compared to R eco , which increased immediately after an initial precipitation event and turned the ecosystem into a carbon source (Fig. 5).Fig. 5 represents year I with the driest spring months to facilitate understanding of GPP and R eco responses to precipitation at the beginning of the growing season after a long drought period.Both R eco and GPP were higher at the LG site as compared to the EG site (Fig. 4) with the absolute values of the GPP/R eco ratios revealing a net carbon uptake (GPP/R eco > 1) during most of the growing season (Fig. C3).
The evolution of NDVI followed the SWC pattern with the lowest NDVI values occurring during spring (September-November) and the maximum NDVI values were measured in mid-February (0.63 at the LG site and 0.72 at the EG site) (Fig. 4).Mean annual NDVI values were slightly higher in the EG site compared to the LG site in the growing seasons, although the difference was not statistically significant (p < 0.05).The NEE of the studied ecosystems is clearly correlated with the P and NDVI (Fig. C5).The EG site demonstrated stronger correlations (R = 0.84 (NDVI); R 2 = 0.95 (P)), whereas the LG site had slightly weaker relationships (R 2 = 0.68 (NDVI); R 2 = 0.73 (P)).The relationships are highly linear at both sites.
The years I and II had similar annual precipitation but different distribution throughout the year (53 % of precipitation from November to February in the year I compared to 78 % in the year II).It resulted in higher NDVI and carbon sequestration in the Year II (Fig. 4, Table C3).Year III was the wettest year (longest NDVI peaks and highest carbon sequestration rates), while year IV represents the driest year with the lowest carbon uptake.
The relationship between daytime means of carbon fluxes (g C m − h − 1 ) during the growing season and the corresponding means of SWC (%) (Fig. C4) from all years of measurement showed stronger correlations (R 2 ranged from 0.53 to 0.73 for individual years) for the LG site, while weaker relationships (R 2 ranged from 0.40 to 0.45) were detected for the EG site.At both sites, the relationship appeared to be highly linear in the range of 5 to 15 % soil water content.Under wetter conditions, i.e.SWC > 15 %, NEE values plateaued with increasing SWC and even showed a tendency to decline at high water content (>25 %).The scatter in NEE increased considerably with increasing SWC.In the dry season, both NEE and SWC were mostly invariant, thereby making it impossible to determine how much of the variability in NEE can be explained by SWC.

Carbon balance
In our analysis, we focused on the carbon exchange dynamics of these understudied sites, compared CO 2 exchange between the LG and EG sites in the first two years of measurement to demonstrate the effect of the past overgrazing and in the last two years to distinguish the impact of current livestock grazing (EG site under heavy grazing regime) (Fig. 6).After years I and II (sum of the years I and II), the difference in NEE between sites was statistically significant (84 g C m − 2 with an annual mean uncertainty of 42 g C m − 2 for the LG site and − 4 g C m − with an annual mean uncertainty of 41 g C m − 2 for the EG site), whereas after the last two years (sum of the years III and IV) the difference in carbon exchange was no longer statistically significant (− 7 g C m − 2 with an annual mean uncertainty of 39 g C m − 2 for the LG site compared to − 32 g C m − 2 with an annual mean uncertainty of 48 g C m − 2 for the EG site) (Fig. 6).Further, the cumulative GPP values in the years I and II were higher for the EG site, whereas in the years III and IV they were higher at the LG site (Table 1).
Carbon budgets were estimated on monthly, seasonal, and annual scales to demonstrate the carbon source/sink strength of the studied sites (Figs. 6, C3, Tables C3, 1).During the dry seasons, the sites were statistically neutral with the LG site releasing on average 6 g C m − 2 month − 1 , while the EG site emitted a net amount of 5 g C m − 2 month − 1 .The highest carbon release during the dry seasons for both sites was measured in December 2016 (year II).On a monthly basis, dry season carbon sequestration was observed only around the occurrence of high precipitation events (July-August 2016 with cumulative precipitation 35 mm), but this was not reflected in seasonal budgets (Table C3).Although only slight differences were found between the years I-IV in the dry seasons, significant fluctuations were observed during the growing seasons with mean seasonal carbon uptake (i.e., negative NEE) of 24 ± 15 g C m − 2 and 49 ± 19 g C m − 2 for the LG and EG sites, respectively.The carbon sequestration periods varied throughout the measurement period in length (2 to 6 months) and in the strength of carbon sequestration.During the growing seasons of years I and IV, carbon release was observed on a seasonal basis.In comparison, years II and III showed enhanced carbon sequestration with the highest carbon uptake in the year III for both sites (Table 1).
Annually integrated cumulative NEE varied from carbon sink to source over the four years measurement period (Fig. 6).Only year III of the four years had negative NEE (i.e. were net carbon sinks at the annual timescale) for the LG site and years II-III for the EG site.In the year I, the LG site had higher carbon release compared to the EG site (Table 1).In the year II, however, slight carbon sequestration was measured for the EG site compared to a carbon release at the LG site.In the year IV, impacted by drought, carbon uptake was limited by water availability even during the peak of the growing season in summer (Table 1).After the four years measurement period, the total cumulative NEE indicated a carbon release of 82 g C m − 2 with an annual mean uncertainty of 40 g C m − 2 for the LG site while the EG site was nearly carbon neutral with a sequestration rate of 36 g C m − 2 and an annual mean uncertainty of 44 g C m − 2 (Fig. 6e).Fig. 6f shows the negative trends during the growing seasons of the years I and II, while after the reintroduction of livestock grazing to the EG site we can observe positive peaks at the beginning of the growing season in the years III and IV indicating temporarily higher uptake at LG.

Table 1
Seasonal (growing (January-May) and dry (June-December) seasons) and annual net ecosystem exchange (NEE), gross primary production (GPP), ecosystem respiration (R eco ), and precipitation (P).Years I-IV are defined as a hydrological year (November-October).

NEE (g C m
LG   C4.

Primary and secondary drivers of daytime NEE during the growing season
The ANNs trained with all available 19 climatic variables produced a R 2 value of 75 % at the LG site and 72 % at the EG site.This benchmark means that over 70 % of the total variability of daytime NEE in the halfhourly dataset can be explained with this set of drivers (Fig. 7).
Most of the variability (about 50 %) can be explained by the radiation measurements (Fig. 7a, b).Thus R g can be considered the primary climatic driver of NEE at both sites.The significance of the other climatic controls was examined by training the ANNs using R g plus one additional driver at a time.It was found that the SWC was the most relevant secondary control for the daytime NEE response, with slightly greater importance at the LG site as shown in Fig. 7c and d.However, the basic patterns of the primary and secondary climatic controls are very similar between the two sites.This finding highlights the substantial similarity between the ecosystem responses of the examined sites, emphasizing their comparable reactions to the investigated climatic controls.

Discussion
This study provides new insights into the poorly understood NEE of CO 2 in semi-arid grazed South African ecosystems.Our original predictions, based on the literature, were only partially supported and in some cases were strongly refuted.Overall, the sensitivity of the EC systems to capture appreciable differences in NEE suggests that these could be used to test under even better controlled conditions the many idiosyncratic and anecdotal models linking grazing regimes to NEE.We discuss our three main predictions below.

Impacts of overgrazing via altered vegetation cover
After two years of measurements, our results indicate that the EG ecosystem was carbon neutral (− 4 g C m-2 with an annual mean uncertainty of 41 g C m − 2 ), while the LG ecosystem was a net carbon source (84 g C m − 2 with an annual mean uncertainty of 42 g C m − 2 ), thereby refuting the prediction that the LG site would be the stronger carbon sink.The significantly higher carbon sequestration rates at the EG site compared to the LG site in the years I and II indicate that a long resting period (2007-2017) from long-term disturbance, combined with changes in species composition caused by historical overgrazing, resulted in a somewhat higher carbon sink potential, even after heavy grazing was resumed.Although the differences in NDVI between EG and LG sites were not statistically significant, the EG site showed slightly higher NDVI peaks during the resting period, which suggests regrowth of biomass.This finding is consistent with reports from other studies about increases in canopy cover after a long resting period (Chen et al., 2014;Mofidi et al., 2013), but does not accord clearly with common understanding that greater production is associated with increasing biodiversity.There are some instances where a decrease in biodiversity may result in higher levels of biomass production.According to a study by Cardinale et al. (2011), the association between biodiversity and biomass production was intricate and relied on the unique composition of species in the particular ecosystem, and this study supports this contention.
The near-identical responses of micrometeorological drivers for both sites, as shown in Fig. 7, leads us to the assumption that the small differences observed between the sites may be attributed to different levels of biodiversity.It is likely that differences in plant growth form dominance were important, i.e. the relative abundance of grasses and shrubs may explain some of the differences in carbon sink strength between the sites, particularly via differences in response to pulse rain events.The natural state of the EG site was disturbed by heavy livestock grazing followed by a change in species composition to unpalatable grass species (nearly all palatable shrubs have been extirpated), which may have increased the resilience of NEE during the drought conditions (Year IV).Zhou et al. (2012) showed that the fine root biomass of perennial grasses was more abundant in the upper soil layer compared to shrubs.It implies that grasses are stronger competitors for water as they are able to absorb water faster due to the abundant fine-root biomass in the upper soil layer compared to shrubs with deep-root systems, especially in water-limited ecosystems with pulsed precipitation, and during short events when only little rainfall was recorded.This corresponds with our finding that the EG site turned into a carbon sink faster (years I and II) in response to high precipitation than the LG site (Fig. 6).The difference was likely due to the vegetation cover at the EG site, being dominated by grasses, and thus able to respond more quickly to rain events due to shallow-root systems that use discontinuous and erratic water sources in the upper soil layers compared with the shrubs deep-root systems that use water in deeper soil layers (Canadell et al., 1996;Zhou et al., 2012).Aristida diffusa is a drought-tolerant dominant grass at the EG site that managed to survive the intense summer drought, while net mortality of other grass species was observed (du Toit and O'Connor, 2020).Furthermore, Zhou et al. (2012) found that total soil organic carbon storage was higher for perennial grasses (higher soil organic carbon inputs through primary production and slower return of carbon through decomposition) than for shrubs.
Our findings do not accord with previous studies indicating that species richness positively impacts ecosystem biomass production (Balvanera et al., 2006;Cardinale et al., 2007).It has been suggested that the positive impacts may be due to higher number of productive species or by reciprocal effects between species, which allow the use of available resources from the environment more efficiently.In our study, the less diverse EG site was found to be more productive than the LG site, with a possible role for plant growth form composition differences overriding assumed biodiversity effects on NEE.It appears that the abundance of species of growth forms that are more resilient and recover more quickly after drought at the EG site is the determining factor in such a water-limited ecosystem.Adler et al. (2011) found no clear relationship between productivity and site diversity in 48 herbaceousdominated plant communities on five continents.While we measured higher carbon sequestration at the previously overgrazed EG site, it is important to note that this site is unfavorable from the perspective of the current livestock grazing system.Livestock grazing was defined as the main driver of plant species composition in the semi-arid Karoo ecosystems (Kraaij and Milton, 2006;Tidmarsh, 1951).Seymour et al. (2010) reported that 20 years of recovery period in the Karoo degraded ecosystems restored grazing potential, while not returning all palatable species.Leu et al. (2014) showed that heavily grazed shrubland (northern Negev, Israel) can be recovered within <16 years by implementing strict conservation management.

Impacts of lenient temporary vs. continuous heavy livestock grazing
The initially observed substantial differences in NEE between the studied sites established for years I and II were intriguing.However, when heavy grazing was resumed in the years III and IV, it became evident that NEE differences between the sites had decreased, yet the EG site continued to demonstrate high carbon sequestration.Again, these findings did not support our expectations.It appears that when grazing was resumed, the difference in NEE between sites was diminished, and differences in grazing intensity between the two sites did not reflect as differences in carbon fluxes.In addition, cumulative GPP was higher in the LG site in the years III and IV, whereas in the years I and II it was higher in the EG site (Table 1).Based on previous studies, we expected that the impact of heavier grazing would affect the carbon sequestration capacity via reduced canopy cover and altered species composition (Tagesson et al., 2015;Yan et al., 2017).Continuous grazing may also affect carbon fluxes via a reduction in aboveground biomass (AGB) (Na et al., 2018) or decreasing the soil C and N storage (Yan et al., 2014).Ma et al. (2019) found a strong correlation between grazing intensity and vegetation indexes (AGB and NDVI).In this study, no significant difference in NDVI between grazing intensities was found.The unexpected results may be due to the relatively high rainfall in the year III; the effects of grazing, which is in any case low-intensity compared to the systems where most previous studies are conducted, may be stronger when coupled with drought stress in the highly seasonal systems.Furthermore, it needs to be noted that the ground state of LG and EG at the beginning of the second comparison period, i.e., at the start of year III, was not the same.Legacy impacts of past heavy grazing were already visible through altered plant species composition, possibly supporting continued resilience of NEE in the EG site with resumed heavy grazing in years III and IV.A longer term approach would be needed to determine if the two sites will develop diverging CO 2 exchange if current grazing schemes are maintained.These findings for years III and IV are among the first to suggest such resilience in an apparently degraded ecosystem even after reintroduction of continuous heavy grazing at many times the recommended rate.

Water availability as the main driver of inter-annual variability of carbon fluxes
While we could attribute differences in CO 2 exchange between sites to grazing intensity, sums and their variability on an annual basis were still highly dominated by the amount of available water, as reported in previous studies conducted in similar systems (Merbold et al., 2009 and references therein).During the growing seasons, carbon fluxes were negative in the daytime (correlated with light intensity throughout the day) and positive at nighttime (switched to net carbon release after the sunset).At the beginning of the growing season, during the onset of summer precipitation events, an immediate response of R eco to even small wetting events and plant germination was observed (i.e., R eco began to increase rapidly with precipitation of 9 mm during the transition period) (Fig. 5).Additionally, these responses seem to be influenced not only by water availability, but also by the physiological recovery of the ecosystem from previous dry periods.The first precipitation event also caused an increase in GPP to some extent, but the NEE indicated a comparatively small carbon source during the plant growth initiation phase due to high respiration.However, during the peak of the growing season, GPP exceeded R eco and the ecosystems turned into a carbon sink.This is consistent with other findings in the semi-arid ecosystems (Kutsch et al., 2008;Merbold et al., 2009).During the dry seasons, the ecosystems were not physiologically active due to the limited water availability.These findings are in line with other studies conducted in African ecosystems (Brümmer et al., 2008;Räsänen et al., 2017).Mean carbon fluxes showed the highest carbon uptake and release rates in the year III, followed by the year II.This is in line with the temporal distribution of precipitation and NDVI indexes.Based on NDVI data, year II stands for the highest NDVI peak, whereas year III (the wettest year) represents the longest NDVI peaks, and years I and IV represent the lowest NDVI peaks due to the distribution (year I) and the deficit of precipitation (year IV).
At the annual scale, the carbon sink strength in the studied sites varied depending on the climatic conditions.Years I and II had similar precipitation rates (373 mm and 370 mm) (Fig. 2d, Table B1).The LG site acted as a carbon source in both periods but in the year II with lower strength compared to the year I, whereas the EG site was a carbon source in the year I and carbon sink in the year II (Fig. 6).The differences in the carbon sequestration rates between years I and II were determined by different precipitation distribution which resulted in fever days of continuous NEE < 0 in the year I.During the transition period (November-December) and the peak of the growing season (January-March) both sites received ~64 % of annual precipitation in the year I compared to ~88 % in the year II (Table B1).Thus, in the growing season of the year II, both sites were net carbon sink ecosystems, while in the year I, the LG site was a small source of carbon and the EG site was carbon neutral.At the end of the growing season (April-May), both sites received higher amounts of precipitation in the year I than in the year II (Fig. 4, Table B1).Fig. 4 shows that the same amount fell within two rainy periods in the year II causing a longer and delayed soil moisture decay leading to more days of continuous NEE < 0 than in the year I, where more equally distributed rain led to a faster dry up.In the year III, both sites acted as carbon sinks due to the increased precipitation and, as a result, the longest number of continuous NEE < 0 days.In the year IV, the studied ecosystems acted as carbon sources with the highest CO 2 release due to the precipitation deficit compared to the long-term mean annual precipitation (approximately 25 % lower).Our results agree with other, previously mentioned studies, which observed similar seasonal dynamics of NEE in their respective regions (Brümmer et al., 2008;Tagesson et al., 2016).In summary, antecedent vegetation state, soil condition, and precipitation distribution are evidently significant factors influencing ecosystem responses.
In our studied sites, mean annual NEE were 21 g C m − 2 year − 1 for the LG site and − 9 g C m − 2 year − 1 for the EG site.Similar annual NEE of − 98 g C m − 2 year − 1 to 21 g C m − 2 year − 1 were observed by Scott et al. (2010) in the Kendall grassland, USA with mean annual precipitation of 345 mm (63 % in the summer months).Räsänen et al. (2017) observed annual carbon budgets of − 85 g C m − 2 year − 1 , 67 g C m − 2 year − 1 and 139 g C m − 2 year − 1 (2011-2013) in the Welgegund atmospheric measurement station grassland ecosystem, South Africa (540 mm annual P).The annual NEE ranges at the South African Kruger National Park were from − 138 g C m − 2 year − 1 to 155 g C m − 2 year − 1 with a mean annual precipitation of 550 mm (Archibald et al., 2009).
In a broader context, Valentini et al. (2014) reported a small carbon sink of 0.61 ± 0.58 Pg C year − 1 (− 20 g C m − 2 year − 1 ) on the whole African continent on an annual basis by averaging out all the estimates.However, complete and accurate estimation of the carbon budget for the African continent is not available due to large gaps in carbon fluxes data from in situ measurements for many African ecosystems.Thus, understanding the dynamics of carbon fluxes in different types of ecosystems is crucial to underpin these larger-scale estimations.

Conclusions
The net exchange of CO 2 in semi-arid Karoo ecosystems in South Africa was studied at two adjacent sites under different management, but practically identical climate regimes.Our approach allowed for a comparison of past and current grazing effects.We found a significant difference between a site under conventional, i.e. temporary lenient grazing, and a site that had been overgrazed in the past but was rested for the previous 10 years.Intriguingly, the conventional site was a larger CO 2 source (84 g C m − 2 with an annual mean uncertainty of 42 g C m − 2 ) than the overgrazed site (− 4 g C m − 2 with an annual mean uncertainty of 41 g C m − 2 ) over the first two years of observation.When sheep were reintroduced to the previously overgrazed site at high density, differences in NEE between these two systems was reduced, but NEE in the EG site showed remarkable resilience, remaining comparable to that in the LG site with far lower grazing pressure.The role of plant species composition and their relative coverage is an important status indicator for the site, but plant growth form composition differences emerge as critical.It appears that optimally managed sites with high abundance of palatable grasses are favorable for sustainable grazing, but less efficient for net CO 2 sequestration.These findings may have important implications on the design of future land management strategies in the region, but need some increase in replication covering more ecosystem types and management systems.It is certainly confirmed that the amount and distribution of rainfall play a crucial role in how the Karoo ecosystems switch from net carbon sources in dry years to net carbon sinks in years with above-average rainfall.

CRediT authorship contribution statement
CB, GF, and OR conceived the study.JP, KM, JKJ, JdT installed and maintained the eddy covariance flux towers and collected the raw data.

Table C1
Summary of mean diurnal carbon fluxes in mg C m − 2 h − 1 of the dry (June-December) and growing (January-May) seasons for years I-IV.Years I-IV defined as hydrological year (November-October).

Table C2
Daily cumulative net ecosystem exchange (NEE), ecosystem respiration (R eco ), and gross primary production (GPP) in g C m − 2 d − 1 .Years I-IV defined as hydrological year (November-October).Fig. C3.Absolute monthly ratios of gross primary production (GPP) versus ecosystem respiration (R eco ) in the (black) lenient grazing (LG) and (grey) experimental grazing (EG) sites.GPP/R eco > 1 ratios revealing a net carbon uptake.

Table C3
Monthly cumulative net ecosystem exchange (NEE) in g C m − 2 .Years I-IV defined as hydrological year (November-October).
semi-arid South African ecosystems has not been extensively studied in relation to the Net Ecosystem Exchange (NEE) of carbon dioxide (CO 2

Fig. 1 .
Fig. 1.South African biomes and location of the studied sites marked as a red circle with wind roses on the right side and pictures in growing and dry seasons for the lenient grazing (LG) (top) and the experimental grazing (EG) sites (bottom) (OpenAfrica: Department of Agriculture Fisheries and Forestry, 2013).

Fig. 2 .
Fig. 2. Annual and seasonal variations of meteorological variables for the entire measurement period: (a) hourly (light red) and daily (dark red) means of air temperature (Tair), (b) hourly (light blue) and daily (dark blue) means of relative humidity (RH), (c) hourly (light yellow) and daily (dark yellow) means of global radiation (R g ), (d) daily means of soil water content (SWC, left) and cumulative precipitation (P, right), (e) normalized difference vegetation index (NDVI) produced on 16-day intervals with red and blue shadows indicating the grazing periods for the LG (blue) and EG (red) sites.Grey shadow indicates the growing season.

Fig. 3 .
Fig. 3. Mean diurnal carbon fluxes in the dry (a-d) and growing (e-h) seasons for years I-IV (top to bottom) in the (blue) lenient grazing (LG) and (red) experimental grazing (EG) sites.Shaded area indicates ± 1 standard deviation.

Fig. 4 .
Fig. 4. Daily net ecosystem exchange (NEE) and partitioned component fluxes (i.e., gross primary production (GPP), ecosystem respiration (R eco )) on the left axis, cumulative precipitation (P cum), monthly mean soil water content (SWC) and normalized difference vegetation index (NDVI) on the right axis across different grazing intensities for (a) lenient grazing (LG) and (b) experimental grazing (EG) sites.The grey background areas represent growing seasons.

Fig. 5 .
Fig. 5. Gross primary productivity (GPP) and ecosystem respiration (R eco ) response to precipitation (P) at the beginning of the growing season in the year I for the lenient grazing (LG) and experimental grazing (EG) sites.Vertical lines indicate changes in the R eco and GPP response to rainfall.

Fig. 6 .
Fig. 6.Annual cumulative net ecosystem exchange (NEE) (a-d) for the years I-IV with vertical lines indicating the highest net respiration point throughout the year when the ecosystem turns into a carbon sink and shaded areas that represent ± uncertainty (see Section 2.3.4),(e) four years cumulative NEE in the (blue) lenient grazing (LG) and (red) experimental grazing (EG) sites and (f) cumulative difference between sites (EG NEE − LG NEE) with shaded areas that represent livestock grazing periods at the LG (blue) and EG (red) sites.

Fig. 7 .
Fig. 7. R 2 performance of the ANN models trained with a single climatic driver at a time (a, b) and R g plus secondary climatic driver climatic driver (c, d) for the LG site (a, c) and EG site (b, d).The dotted line is the R 2 performance with using all 19 drivers as a benchmark.The (partly invisible) error bars show the standard deviation of 10 ANN training scenarios.For variable description see TableC4.

Fig. C1 .
Fig. C1.Mean diurnal carbon fluxes fingerprint in the growing seasons (January-May) in the lenient grazing (LG) and experimental grazing (EG) sites.Blue color represents net carbon uptake, while red indicates net carbon release.

Fig. C2 .
Fig. C2.Temporal dynamics of hourly carbon fluxes for the entire measurement period in the (a) lenient grazing (LG) and (b) experimental grazing (EG) sites.

Fig. C4 .
Fig. C4.Correlation between soil water content (SWC) and daytime net ecosystem exchange (NEE) in the growing seasons (January-April) of years I-IV for the (a) lenient grazing (LG) and (b) experimental grazing (EG) sites.

Fig. C5 .
Fig. C5.Correlation between (a) means of normalized difference vegetation index (NDVI) and cumulative net ecosystem exchange (NEE), (b) precipitation (P) and cumulative net ecosystem exchange (NEE) in the transition period and peak of the growing season (November-February) of years I-IV for the (blue) lenient grazing (LG) and (red) experimental grazing (EG) sites.

Table C4
Description of variables used in ANN models.