Greenhouse gas emissions from New Zealand sheep and beef farms

This project was conducted to develop a dataset containing animal policies, production efficiencies and greenhouse gas (GHG) emissions from a large number of sheep and beef farms, and to examine the relationships between farm management and physical constraints, and GHG emissions. The farm-scale model Farmax was used to estimate feed inventories, livestock policies and GHG emissions from 170 New Zealand sheep and beef farms. Emissions were calculated from Farmax outputs using Agricultural Inventory methodology. A quantitative approach was used to group farms, based on physical constraints and management attributes. Mean annual biological GHG emissions from the modelled farms were 3,662 kg CO2 equivalents (CO2-e) per effective hectare (ha), and ranged from 157 to 7,096 kg CO2-e/ha. As stocking rate and animal product (wool + net carcass weight) per effective ha increased, GHG emissions increased. However, there was considerable variability in the data, whereby farms with GHG emissions of approximately 4,000 kg CO2-e/ha had an almost three-fold difference in animal production (range 129 to 360 kg/ha). This work provided a holistic assessment of the farmscale drivers of GHG emissions and a comprehensive database or baseline from which future trends in farmscale GHG emissions can be established.


Introduction
Agriculture was the single largest contributor (48.1%) to New Zealand's total greenhouse gas (GHG) emissions in 2019 (Ministry for the Environment 2021). Within the agricultural sector, sheep and beef livestock contributed almost 50% of emissions from enteric fermentation, and 18.3% of agricultural soils emissions, i.e., contribution of direct nitrous oxide emissions from urine and dung deposited by grazing animals (Ministry for the Environment 2021). However, the contribution of the sheep and beef sector to New Zealand's agricultural emissions has declined by 30% over the last 30 years, primarily associated with a reduction in sheep numbers (Beef + Lamb New Zealand Economic Service 2019). Despite this reduction, the sector has maintained similar levels of production and has doubled the value of its exports since 1990 (Beef + Lamb New Zealand 2019).
Although New Zealand is a small contributor to GHG emissions on a global scale, it exports the vast majority of its meat production and has an opportunity to demonstrate leadership internationally through innovation in the agricultural sector. Sheep and beef farms in New Zealand are diverse and located over a range of landscapes, where each farm has different natural and capital assets, due to climate, location, aspect, altitude, slope, soil type (natural assets), along with previous investment in fencing, stock water, capital fertiliser, pasture improvement and animal genetics (capital assets). Despite their significance, little information exists describing the impact of the range in on-farm natural and capital assets and management decisions, which are key drivers of GHG emissions across these diverse businesses (Mackay 2008;Harrison et al., 2016). The inherent differences between these critical factors across commercial sheep and beef farms means that a large dataset is required to represent the diversity of these systems.
He Waka Eke Noa is a recently established primary sector partnership (Primary Sector Climate Action Partnership) involving a five-year programme of work towards the implementation of a framework that aims to reduce GHG emissions while building farm resilience to climate change (https://hewakaekenoa.nz/ about/). This programme includes several milestones, e.g., by the end of 2022, all New Zealand farmers and growers need to know their annual on-farm biogenic GHG (methane + nitrous oxide) emissions number, and by 2025, all farms have a written plan to estimate and manage GHG emissions (https://environment.govt. nz/what-government-is-doing/areas-of-work/climatechange/he-waka-eke-noa-primary-sector-climateaction-partnership/#programme-milestones).
The objectives of this study were firstly to develop a dataset containing animal numbers, productive and reproductive efficiencies and GHG emissions from a large number of sheep and beef farms located throughout New Zealand (n = 170 farms). The second objective was to examine the relationships between variables that describe farm, physical constraints and GHG emissions. In more general terms, the aim was to provide an assessment of the farm-scale drivers of GHG emissions and a comprehensive baseline from which future trends in farm-scale GHG emissions can be established.

Modelling approach
Key farm characteristics affecting GHG emissions were assessed from a total of 170 sheep and beef farms using the farm systems model Farmax (Science Edition v. 7.2.2.46). Data were modelled for a single year for each farm, either 2015/2016 or 2016/2017. Animal numbers and animal performance (e.g., liveweight (LW) gains and reproductive performance), along with feed on offer, were used to parameterise Farmax, in search of attaining farm 'feasibility' (i.e., making sure that feed on offer matched that required for animal performance).
In the assessment of farm characteristics driving GHG emissions, a large number of explanatory variables were considered. Although closely related, these variables were broadly categorised as those associated with natural and capital assets (e.g., farm cultivatable areas as a proportion of total area) and management-driven (e.g., N fertiliser applied, animal performance).
Annual biological GHG emissions per effective ha (i.e., per unit of area in grazing or growing fodder; hereafter 'GHG emissions' (methane + nitrous oxide), in kg CO 2 -e/ha, unless stated otherwise) were calculated from Farmax outputs using the Agricultural Inventory Model (AIM) equations (Ministry for Primary Industries 2019). Sources of methane (CH 4 from enteric fermentation and manure) and nitrous oxide (N 2 O from direct and indirect sources) were calculated based on published emission factors used by New Zealand's National Inventory calculations. It is important to note that these calculations were made prior to the April 2020 changes to urine emission factors for hill country N 2 O emissions in AIM (Ministry for Primary Industries 2021). Global warming potentials of these GHG (25 and 298 for CH 4 and N 2 O, respectively) were used to convert emissions to CO 2 equivalents (CO 2 -e). Carbon dioxide emissions were not considered in this modelling exercise.

Farm data sources
Beef + Lamb New Zealand (B+LNZ) provided anonymised production and financial data from farms in its Sheep and Beef Farm Survey (B+LNZ; https:// beeflambnz.com/data-tools/sheep-beef-farm-survey), which was used to model 105 farms throughout New Zealand. A further 65 anonymised farm models were provided without assigned B+LNZ farm class nor financial benchmarking data. The 170-farm data set was used in the analysis of the physical relationships of the farm system, whereas the 105-farm subset of the data was used in the analysis of financial relationships. The 105 farms with reliable financial data were selected to ensure that they were representative of the range of New Zealand farm types by B+LNZ farm class (B+LNZ; https://beeflambnz.com/data-tools/ farm-classes). The focus of this paper was on certain natural and capital assets/constraints and managementdriven variables which may drive biological GHG emissions. Financial variables were not considered in this modelling exercise.

Farm clusters B+LNZ farm business classes
As mentioned above, B+LNZ has developed a set of industry standard farm business classes based on island location and business type. A total of 148 farms had been identified by their B+LNZ farm class and were grouped accordingly. A further 22 farms were not identified by their farm class and were not easily assignable due to the anonymity surrounding farm information.

Farm feed groups
New Zealand sheep and beef farms vary from extensive, high country breeding systems to irrigated finishing systems. Comparing such diverse farms reduced the ability to identify physical and efficiency indicators that may influence GHG emissions from the farm system. Therefore, clustering based on inputs and physical properties was necessary to aid in understanding the impact of farm natural and capital assets and efficiency indicators on GHG emissions. A quantitative approach was used to group farms based on either physical constraints or management attributes. A regression tree analysis was performed to rank the variables, based on their separation and magnitude of contribution to GHG emissions. Histograms were used to enable quantitative assessment of the highest-ranking variables and identify appropriate biological ranges for each grouping. Briefly, variables that represented physical constraints of the farm, climate, weather and livestock performance were evaluated to develop the 'farm type' groups. However, none of the groups based on combinations of physical constraints provided adequate separation in annual GHG emissions. This changed the focus onto clustering based on annual stocking rate and feed intake (herein by 'feed'). Criteria for clustering included annual stocking rate as stock units (SU) per effective ha (SU/ha; <5, 5-10, >10; 1 SU = 550 kg DM intake (Woodford and Nicol 2004)) and annual livestock feed intake (kg DM/ ha; <4500, 4500-6800, >6800). By combining these two variables, five feed groups were formed, ranging from low stocking rate/low feed intake (Feed Group 1) to high stocking rate/high feed intake (Feed Group 5).

Statistical analysis
A systematic assessment of the effect of farm variables on GHG emissions was conducted graphically using the statistical packages Genstat (VSN International 2015) and R v. 4.02 (R Core Team, 2013). Analysis of >150 explanatory variables was undertaken in an earlier iteration of the modelling (n = 125 farms), including post-weaning mortality for sheep and cattle, product per kg of LW wintered, LW gain for all stock classes per effective ha, and ratio of LW for breeding and trading livestock classes. As a result of this iteration, the search was narrowed to <20 explanatory farm variables to allow for both one-to-one relationships (one y, one x) and multivariate analyses (one y, several x). In a pathway to separate the impact of selected individual variables on GHG emissions, a correlation matrix heatmap was used to highlight the degree of association between variables on a one-to-one basis. Subsequently, principal component analysis (PCA) was performed using R v. 4.02 (R Core Team, 2013). This data reduction technique allowed for multivariate analysis using the same selected explanatory variables.

Results
Farms were grouped based on stocking rate and feed intake (Figure 1b), rather than farm class (Figure 1a). This allowed for the separation of farms into groups with similar GHG emissions and aided in the interpretation of the complex relationships between emissions, farm biophysical variables and management practices. Of the potential nine combinations (three stocking rates x three feed intakes) considered in the grouping process, all farms were in one of five groups: low-low (Group 1), medium-low (Group 2), medium-medium (Group 3), high-medium (Group 4) or high-high (Group 5) (stocking rate-feed intake, respectively) ( Table 1).
Mean annual agricultural GHG emissions (CH 4 + N 2 O) from the modelled farms were 3,662 kg CO 2 -e/ ha, and ranged from 157 to 7,096 kg CO 2 -e/ha ( Table  2) Annual GHG emissions and a) B+LNZ farm business class, and b) feed group (grouped by annual stocking rate and feed intake). B+LNZ farm business class: 1 South Island high country; 2 South Island hill country; 3 North Island hard hill country; 4 North Island hill country; 5 North Island finishing; 6 South Island finishing breeding; 7 South Island finishing; 8 South Island mixed finishing) for New Zealand sheep and beef farms (n=170 farms). Feed groups: low-low (Group 1), medium-low (Group 2), medium-medium (Group 3), high-medium (Group 4) or high-high (Group 5) (stocking rate-feed intake, respectively). beef production systems. As stocking rate and animal production (wool + net carcass weight; per effective ha) increased, GHG emissions rose. However, there was considerable variability in the data, whereby farms with GHG emissions of approximately 4,000 kg CO 2 -e/ha had an almost three-fold difference in animal product (range 129 to 360 kg/ha). Similar variation was seen in other efficiency metrics, e.g., live weight gain, lamb weaning percentage and ewe efficiency (Table 2).
To separate the impact of selected individual variables on GHG emissions, a correlation matrix heatmap was used to highlight the degree of association between variables on a one-to-one basis ( Figure 2). As expected, stocking rate (SU/effective ha) and total and pasture feed intake (kg DM/effective ha) were the strongest drivers of GHG emissions per ha (darkest shade of red within the green box associated with Pearson correlation coefficients >0.75) (Figure 2).
Two components (PC1 and PC2) captured most of the existing variance in the PCA. The resulting loading biplot covered two dimensions that explained most of the variation (i.e., PC1 explained 50% and PC2 11% of the variation (Figure 3). Feed groups were added in the background and are represented by the coloured bubbles (orange to purple, matching the feed group legend in Figure 3 and shown in Table 1). A number of key drivers of GHG emissions followed the same direction and extent as the big arrow in blue, denoting GHG emissions per effective ha. These key (and highly interrelated) variables included total feed intake and stocking rate (red dashed box in Figure 3). As mentioned above, PC1 explained 50% of the variation, and, within the 'purple' bubble (feed group 5), total DM eaten and stocking rate influenced that dimension.

Discussion
This work provided a holistic assessment of the farmscale drivers of GHG emissions and a comprehensive status or baseline from which future trends in farmscale GHG emissions can be established. No attempt was made to evaluate alternative mitigating scenarios Table 2 Mean (± standard error of the mean) annual GHG emissions and selected farm variables by feed group based on stocking rate and total feed intake (n = 170 farms). within each farm or group. The short-termed nature of the dataset collected (one farm, one year) implies that caution is required when extrapolating these results to longer periods of time, in circumstances other than those captured in the current dataset (weather, market or otherwise) and to all commercial sheep and beef farms in New Zealand (n = 11,300 farms; Beef + Lamb New Zealand Economic Service 2019). Since feed groups represented different stocking rates and feed intake, and these are known to directly drive methane emissions, average GHG increased linearly from feed groups 1 to 5. For many sheep and beef farmers, the natural and capital assets limited the capacity for these businesses to shift from one feed group to another. However, to what extent land managers have the opportunity to change or modify current systems or management to improve farm performance and/or profitability within their feed group (i.e., more production for the same emissions or maintain production while reducing emissions)?

Farm biophysical or efficiency indicator
On average, methane emissions from enteric fermentation and manure accounted for 80% of total GHG emissions from these farms. Methane is generated through ruminant digestion, and emissions are related to feed intake (Hammond et al., 2009;Niu et al., 2018).
Most nitrous oxide is emitted from soils and is affected by urinary nitrogen (N) excretion, other sources of N (fertiliser, biological fixation), vegetation type and cover and soil type (Di and Cameron 2006;de Klein et al., 2019). Given the low amounts of N fertiliser used, a strong relationship between methane and total emissions can be expected for these sheep and beef systems.
The strong relationship between feed intake and GHG emissions has been established in New Zealand and elsewhere, to the extent of providing the basis for GHG emissions calculations (Ministry for Primary Industries 2021). Feed is a consequence of animal demand and what is on offer. Feed demand is driven by physiological process taking place in the grazing animal and is estimated in Farmax as the sum of energy required for maintenance, pregnancy, lactation, growth and activity (CSIRO 2007). The quantity of feed on offer is largely dictated by the pasture and forage available to the grazing animal. This is because livestock are not typically fed other sources of nutrients, nor is N fertiliser used to increase pasture production. In addition to the expected diversity in herd composition and management, some variation in feeding at a given GHG emission level was observed, which indicated GHG efficiency gains on some farms. For the same GHG emissions, lower feed intake tended to be linked to farms with a higher proportion of nitrous oxide emissions, particularly from land areas in fodder crops. This was often offset by higher animal growth rates and better efficiency due to higher quality feed. Making farms more efficient would require a focus on increasing individual animal performance, in conjunction with a decrease in the number of animals to maintain feed intake at the same level or lower.
Ewe efficiency (kg of lamb weaned per kg of ewe mated, expressed as a percentage) ranged from 20 to 85% across the data set. Unlike the clear trend between feed groups and GHG emissions, the spread in ewe efficiency was similar across all the feed groups. The extent to which the spread in data represents management practices, breed differences or differences in genetic merit was not explored.
PCA has gained popularity to highlight strong patterns within complex biological datasets, by capturing essential relationships. These components convey most variation in the dataset, reducing the overwhelming number of dimensions to those that explain most variation. The group of variables rated as having an 'intermediate correlation' (i.e., a much weaker but positive correlation with GHG emissions compared with feed intake and stocking rate) are placed in the same right-hand side of Figure 3 (blue dashed box), but not in the same direction as the blue arrow. These variables included cultivatable area as a proportion of total land, N applied as fertiliser, LW sold per LW wintered, ewe efficiency and lamb weaning 8 3501 Two components (PC1 and PC2) captured most of the existing variance in the PCA. T resulting loading biplot covered two dimensions that explained most of the variation (i

Figure 2
Heatmap to visualize hierarchical clustering of selected farm variables, as they relate to each other and to GHG emissions per ha (in the green box). The dendrograms along the sides show how the variables on rows and columns are independently clustered. The patterns on the heatmap show the degree of association (Pearson correlation) between variables; these range from -1 (dark blue; negative correlation) to 1 (dark red; positive correlation). Values closer to 0 indicate that there is no linear trend between the two variables (n = 170 farms).

Figure 3
Principal component analysis (PCA) of selected farm variables including GHG emissions per effective ha (larger blue arrow) (n = 170 farms). Each individual number represents a farm. Feed groups: Feed groups: low-low (Group 1), medium-low (Group 2), medium-medium (Group 3), high-medium (Group 4) or high-high (Group 5) (stocking rate and feed intake, respectively) Journal of New Zealand Grasslands 83: 225-232 (2021) percentage. These variables provided a measure of farm constraints, N inputs, trading-to-breeding component, and reproductive efficiency (with underlying genetic and nutrition efficiencies). Changing some of these variables within a given farm system could reduce GHG emissions (Cruickshank et al., 2009), but, when considered across a large number of sheep and beef farms, these variables did not reduce emissions on their own; although some have been considered important in terms of GHG abatement by shortening the half-life of the growing livestock. However, for a snapshot (one year) across a large number of farms, this particular variable can become hidden. The current data instead provided a comprehensive baseline from which trends could be identified if additional years were to be modelled for these farms.
The shape and extent of the feed groups in the background vary, and a broader array of farms within feed group 1 was most likely a reflection of the broader range of values within each variable selected. The feed grouping criteria were based on the main drivers of the vectors moving horizontally along PC1 (X axis), with variables such as stocking rate, feed consumption and animal produce explaining most of the variability in GHG emissions. Conversely, effective area and, to a lesser degree, cultivatable area as a proportion of total land and feed conversion efficiency (FCE), explained most of the vertical variability along PC2 (Y axis; Figure  3). For these reasons, the shapes of the bubbles varied, whereby feed groups 2, 3 and 4 were narrower along PC1 but broad along PC2, which showed the variability of factors beyond those contemplated in PC1.
A number of variables were not selected for the correlation matrix heatmap and PCA analysis. Response-type variables, such as individual gases (methane and nitrous oxide) and emission intensity (emissions per unit of animal product) were not used in the final selection, and annual agricultural GHG emissions per effective ha (CH 4 + N 2 O) were the main response variables of interest. Following a similar correlation matrix and PCA approach, Becoña et al., (2014) showed that there is potential to reduce cow-calf GHG emissions through improved grazing management in Uruguay (n = 20 farms and 26 explanatory variables).
Variables relating to cattle were used not in the final selection due to less farms with significant numbers of breeding cattle (99 of 170 farms). The alternative of providing a measure of combined reproductive efficiency for both ewes and cows was considered, but a meaningful unifying metric could not be identified. Thus, different sample sizes (i.e., number of farms with breeding ewes versus breeding cows) and the lack of a meaningful unifying metric resulted in the decision against considering this variable further.
Difficulties in developing a meaningful breeding-to-finishing metric led to the inclusion of LW sold per LW wintered. Originally, ratios, such as traded-to-breeding LW and net gained-to-breeding LW, were analysed, but discarded and replaced with LW sold per LW wintered. The addition of the two trading-to-breeding ratios did not add clarity, are was most likely a reflection of the short-termed nature of the analysis. This and the lack of a (theoretical) steady-state assumption prevented the inclusion of other variables, such as replacement and mortality rates. Although specific mitigation strategies were beyond the scope of this work, the multivariate analysis was consistent with earlier studies showing that feed intake modification is the main option for GHG abatement (Reisinger et al., 2017). However, earlier analysis highlighted the variability in production per ha, which suggested the potential for increased efficiency to deliver to both animal production and gross margins. A continuation of improvements in production efficiency is necessary via a focus on increasing animal product and profit for the same emissions.

Conclusions
The 170 farms in the dataset were from across all B+LNZ farm classes and were selected to be representative of the New Zealand sheep and beef cattle livestock production/farming sector. Analysis of actual farm businesses in this study has provided insight into the complexity within and variability between sheep and beef farms. Total feed production, quality feed intake on sheep and beef farms drives stocking rate and animal production per ha, and these were highly correlated with GHG emissions. This work highlighted the complexity and trade-offs faced by farmers in considering a future low carbon economy. Opportunity exists for a re-analysis of farm emissions once the changes to urine emission factors for nitrous oxide on hill country slopes have been implemented within an appropriate calculator. This work provides a holistic assessment of the farm-scale drivers of GHG emissions and a comprehensive current state of affairs or baseline from which future trends in farm-scale GHG emissions can be established.
Maclean completed the multivariate analysis and Peter Green completed the initial multivariate analysis and the 'heat map' analysis.