Economic evaluation of variable rate nitrogen management of canola for zones based on historical yield maps and soil test recommendations

Canola (Brassica napus L.) is a highly valuable crop for Canada’s economy, making the efficient management of canola a priority. A field-scale study was conducted at ten sites between 2014 and 2016 to evaluate the viability of site specific nitrogen (N) management zones (MZ) based on analysis of historical yield maps and soil test recommendations to improve canola productivity and profitability in western Canada. Treatments included factorial combinations of three canola yield zones (low, average, high) by four N rates, replicated four times at each site. The canola yield function had a quadratic form in each field but the effects of MZ varied between fields with positive effects in only a few fields. When ten site-years data were combined, MZ had positive effects on canola performance. On average, MZ of N fertilizer over ten fields generated between $28 to $65 ha−1 more net revenue (NR) relative to average yield management. Site-years, which reflect farm management and other farm characteristics had significant effects on yield and NR ranging from − $91 to $352 ha−1 compared to a baseline. Nitrogen application under MZs was only reduced by 8% compared to uniform rates. The potential for MZ does exist; however, its effectiveness is highly variable.


Results and discussion
Canola cost of production. Total average annual production costs over the ten sites were $898, $956, and $992 year −1 ha −1 for low, average, and high zones, respectively ( Table 1). The cost differences among zones were due mainly to differences in fertilizer application rates and machinery or transportation costs related to the yield dependent operations. For example, the high zone had higher crop yield and, therefore, higher yield dependent machinery and transportation costs. Fertilizer costs were higher in high zone due to higher application of N rates compared to the average or low zones (Table 1). Operating costs were 66.1, 68, and 69.35% of the total costs for low, average, and high ones, respectively. Fertilizer costs were 16.3, 20.8, and 23.3% of total costs for low, average, and high zones, respectively. Overall, mean fertilizer costs, which made up 20.2% of total cost, were in line with fertilizer costs reported for crop production in western Canada. N fertilizers can comprise more than 31% of the operating costs or 21% of total costs of traditional production systems in western Canada 17 . Canola yield, 27 sites 2014 to 2017. Canola yield was affected by an interaction between fertilizer treatments, management zones and farm site, in analyses of 27 farm fields studied from 2014 to 2017 (Table 2). Yield response to management of N fertilizer was affected by management zones, and field as influenced by covariates such as soil properties, topography, agronomic management, local precipitation and temperature at each siteyear. Canola yield within zones increased from the control to 50 and 100% of recommended fertilizer rates in combined analysis across all 27 sites (Fig. 1B). The yield response to N fertilizer was nonlinear, in the combined analysis for all fields, canola yield did not increase from 100 to 150% of the recommended fertilizer rates within the average and high zones (Fig. 1B). Canola yield was higher in the high zone compared to the low yield zones, with intermediate values for average zones.

Canola yield and NR by sites in 2014.
Canola yield and NR within zones in 2014 were affected by recommended N in Manitoba but not in Saskatchewan fields (Fig. 2). In all zones of 2014-MB1, yield and NR of canola were greater for the 100% recommended N treatment compared to 0% N (Fig. 2). In high and average yielding zones of 2014-MB2 and 2014-SK1, NR was maximized at 0% N and more Ns did not generate additional NR. In the low zone, N for maximum canola yield (and N for maximum NR) were 118 (118), 166 (66), and 36 (26) kg N ha −1 for 2014-MB1, 2014-MB2, and 2014-SK1, respectively. In the average yield zone, N for maximum canola yield (and N for maximum NR) were 208 (172), 177 (58), and 164 (23) kg N ha −1 for 2014-MB1, 2014-MB2, and 2014-SK1, respectively. In the high yield zone, N for maximum canola yield (and N for maximum NR) were 141 (141), 208 (61), and 193 (23) kg N ha −1 for 2014-MB1, 2014-MB2, and 2014-SK1, respectively. These Ns included both soil and N application. The Ns to achieve maximum NR were less than Ns for maximum yield.
Within 50-100% recommended N rates, the significant effect of management zone on canola yield and NR was found only for 2014-SK1 (Fig. 4). Yields in high and average zones (2068 and 1817 kg ha −1 , respectively) were higher than yield in low zone (1258 kg ha −1 ) for the 2014-SK1. Similarly, NRs in high and average zones ($58 and − $22 ha −1 , respectively) were higher than NR in low zone (− $219 ha −1 ).
The linear and quadratic response of predicted yield and NR to N rate (soil plus N applied) for 2014-MB1, 2014-MB2 and 2014-SK1 was regressed with average zone as a baseline (Table 3). Linear and quadratic terms for N were both significant for 2014-MB1. However, coefficients for high zone (difference between high zone and average zone) and low zone (difference between low zone and average zone) were not significant. At 2014-MB2, the effect of N fertilizer on yield was linearly significant along with effect of high zone on both yield (253 kg ha −1 ) and NR ($128 ha −1 ). At 2014-SK1, the regression coefficient for low zone was significant on both yield (− 408 kg ha −1 ) and NR (− $185 ha −1 ).

Canola yield and NR by sites in 2016.
Statistical analysis showed that in the low zone of Manitoba sites in 2016, increasing N to 100% N increased canola yield compared to 0% N; however, increasing N rates did not increase NRs except at 2016-MB3 where increasing N to 100% recommended N generated more NR than 0% N. In the average zone, increasing N from 0 to 100% N increased canola yield at 2016-MB1 and at 2016-MB3 but not at 2016-MB2. Regardless of the locations, increasing N did not increase NR in the average zones of 2016 sites except at 2016-MB3. In all four sites in 2016, application of more N in high zones did not increase NRs, perhaps  Low  0  155  229  87  78  530  833   Low  50  155  230  128  78  574  877   Low  100  155  232  165  78  615  919   Low  150  155  232  204  78  657  961   Mean of low zone  155  231  146  78  594  The significant effect of management zone on canola yield and NR was found for 50-100% N rates ( Fig. 4). At 2016-MB1, canola yield was higher at the high zone (2796 kg ha −1 ) compared to the low zone (2381 kg ha −1 ). At 2016-MB2, canola yield was higher at high zone (1837 kg ha −1 ) compared to low and average zones (1132 and 1460 kg ha −1 , respectively), and yield at the average zone was higher than yield at the low zone ( Fig. 4). At 2016-MB2, NR was higher at the high zone (− $177 ha −1 ) than the low zone (− $405 ha −1 ). At 2016-MB3, yields at the high and average zones (2053 and 2014 kg ha −1 , respectively) were both higher than yield at the low zones (1717 kg ha −1 ) (Fig. 4).
At 2016-MB2, linear and quadratic coefficients for N, high zone (287 kg ha −1 ), and low zone (− 302 kg ha −1 ) were significant on yield (Table 4). At this location, only quadratic coefficients of N, high zone ($165 ha −1 ), and low zone (− $160 ha −1 ) were significant for NR. At 2016-MB3, linear and quadratic coefficients for N were also significant. were not significant at 2014-MB1, but the high zone in 2014-MB2 produced higher canola yield (253 kg ha −1 ) and higher NR ($128 ha −1 ) than the average zone. The low zone in 2014-SK1 produced less canola yield (408 kg ha −1 ) and less NR ($185 ha −1 ) compared to the average zone (Table 3). Results for sites in 2015 showed that zones were not significant except for the effect of low zone on NR in 2015-MB1 (Table 3). Similarly, zones were not significant in 2016 except in 2016-MB2 where high zone produced more yield (287 kg ha −1 ) and NR ($165 ha −1 ) and low zone less yield (302 kg ha −1 ) and less NR ($160 ha −1 ) than the average zone (Table 3). Results also showed that Ns to achieve maximum NR were less than Ns for maximum yield. Other studies have found similar mixed results. Beckie et al. 18 found net returns from precision farming to be about $10 ha −1 greater compared to uniform rate N fertilization of wheat (Triticum aestivum L.), canola, and flax (Linum usitatissimum L.) grown in black soil climate near Prince Albert, Saskatchewan, Canada. Similarly, in Inverleigh, Australia, Torpy 19 found that precision agriculture had a beneficial impact on investment as canola yield differences were higher in paddocks with high spatial variability. A study on the effect of variable rate application on cereal grains in Oklahoma www.nature.com/scientificreports/ and Virginia found that the increase in revenue from precision agriculture was more than sufficient to cover the technology costs 20 . Conversely, Long et al. 21,22 found net returns from precision agriculture of hard red spring wheat in northern Montana to be reduced by up to $27.97 ha −1 compared to net returns under uniform N man- www.nature.com/scientificreports/ agement. Boyer et al. 23 found no statistical differences in wheat yield and net return between variable rate N and uniform rate application in Oklahoma. A study in eastern Washington found no significant differences in winter wheat yield between uniform and variable rate N application but found that less N fertilizer was applied under precision application in two out of three management zones 24 . A three-year study in Saskatchewan, Canada found also little economic rationale for precision agriculture due to the unpredictable nature of wheat yield response to varied N rates 25 . Several studies have found that fertilizer cost savings have a greater impact on the economic benefits of precision farming than do increases in grain yield [26][27][28] .

Canola yield and NR, annual regression analyses for 2014, 2015, and 2016. The quadratic
response of predicted yield to N application for each individual year, assuming fixed effects of zones and location within each year, were regressed (Table 3). In 2014, although the canola yield function was quadratic in form similar to results for yield in combined analyses for 2014-2017 ( Fig. 1B), only the linear coefficient and the coefficients for locations were significant. 2014-MB1 produced 1250 kg ha −1 higher canola yield than 2014-MB2, while 2014-SK1 produced 58 kg ha −1 less canola than 2014-MB2. Canola yield for both 2014-MB1 and 2014-MB2 sites were higher than yield at 2014-SK1 due, in part, to higher growing season precipitations in Manitoba sites. In 2014, growing season precipitation in Manitoba was 17% higher than the growing season precipitation of the same year in Saskatchewan (data not shown). NR for 2014 showed similar patterns with only linear term for N and locations being significant. 2014-MB1 generated $510 ha −1 higher canola NR than 2014-MB2, while 2014-SK1 generated $34 ha −1 less canola NR than 2014-MB2. In 2015, coefficients for linear and quadratic effects of nitrogen and high zone on canola yield and NR and effects of locations were significant but not the effects of low zone. In 2016, both yield and NR functions were quadratic in form and the effects of all coefficients of independent variable on canola yield and NR were significant except the linear effect of N on NR. The high zone produced 120 kg ha −1 higher canola yield or $58 ha −1 higher canola NR than the average zone. Conversely, the low zone produced 142 kg ha −1 less canola yield or $70 ha −1 less canola NR than the average zone. Growing season precipitation in 2016 in Manitoba sites were 68% higher than the growing season precipitation of the same year in Alberta; however, this did not result in more canola yield in Manitoba (data not shown). In all three years (2014, 2015, 2016), coefficients for sites which captured combination effects of site management and site specifics were significant, reflecting different farm management and specifics and including soil and environment at each  www.nature.com/scientificreports/ www.nature.com/scientificreports/ site. Overall, regression results showed that the coefficients for zones though all had expected signs, were only significant for high zone in 2015 and for both high and low zones in 2016. This could be explained by variability of fertilizer rates for MZs between fields, by uniform yield for some sites, and by factors such as soil properties. The ANOVA results for the effects of nitrogen fertilizer, canola management zones, and locations on canola yield and NR were somewhat different than the regression analysis (Table 2). In regression, zones and locations were fixed effects with indicator variables; however, in the ANOVA analysis, no baselines were selected and only the effects of these variables on canola yield or NR were tested. ANOVA results for each individual year showed that the p values for the effects of nitrogen, zones, and locations on canola yield and NR were significant except the effects of nitrogen on NR in 2014 and zones on NR in 2015 were not significant.
Canola yield and NR, combined analyses for 2014 to 2016. The quadratic response of predicted yield to N application for all site-years was regressed (Table 3). Both the linear and quadratic effects of nitrogen on canola yield were found to be significant. The high zone produced 135 kg ha −1 higher, while yield in the low zone was 91 kg ha −1 lower than the average zone. 2016-MB1 was assumed as a baseline and other site-years were compared to this baseline. Coefficients for site-years which capture a combination of effects of year and site management, and site specifics were significant. This reflects different farm management, site specifics and soil and environmental aspects at each site. For example, 2016-MB2 produced 95 kg ha −1 less canola yield than 2016-MB1, while 2015-SK1 produced 177 kg ha −1 more canola than 2016-MB1. Growing season precipitation in 2016 in Manitoba was 14% higher than the 2014 growing season precipitation in Manitoba and 20% higher than the 2015 growing season precipitation in Saskatchewan (data not shown); however, this higher precipitation did not result in more canola yield compared to sites in 2014 and 2015 in Manitoba and Saskatchewan.
For all 3-year (10 sites) combined, the linear and quadratic coefficients for N on canola NR were significant ( Table 3). The high zone generated $65 ha −1 higher canola NR than the average zone. The effect of the low zone on canola yield and NR as compared to average zone was negative but was significant only at 10% p value. The low zone reduced canola NR by $37 ha −1 compared to NR of canola in the average zone. Coefficients for farm sites as compared to 2016-MB1 were significant for all farms, reflecting contribution of different farm management and soil and environmental specifics at each site. For example, 2016-MB2 generated $47 ha −1 less canola NR than 2016-MB1, while 2015-SK1 generated $106 ha −1 more canola NR than 2016-MB1.  (Table 4). For example, in 2014-MB2, 123 kg ha −1 N was applied at 100% N treatment over the three zones, but 134 kg ha −1 N was applied under a uniform application. Or, in 2016-MB2, 128 kg ha −1 N was applied at 100% N treatment over the three zones, but 136 kg ha −1 N was applied under a uniform application. As results showed, the N application was reduced on average by only 7 kg N ha −1 under MZ. This was based on the actual amount of N applied. If, however, N applied was based on a yield maximizing scenario (Figs. 2 and 3), reduction in N and therefore NUE will be more significant. For a yield maximizing scenario, N was on average 18 kg N ha −1 less under management zones compared to uniform application. Studies on the environmental benefits of precision agriculture as compared to conventional fertilizer application found significant benefits with precision agriculture 20,27,29-32 . Holzapfel et al. 30 reported that sensor-based N management of canola grown in Saskatchewan, Canada increased agronomic NUE 33% of the time. Similarly, Raun et al. 20 found that precision agriculture increased NUE by over 15% on average for cereal grains grown in Oklahoma and Virginia. Diacono et al. 31 reported up to 368% increased in NUE for wheat using real-time sensor-based N management systems as compared to conventional farming practices. www.nature.com/scientificreports/

Conclusions
Variable application of N fertilizer based on management zones was not effective in all fields, due to variability between sites. Overall, when ten fields were combined, management zones of N fertilizer increased canola NR by up to $65 ha −1 compared to average yield management. Site-year or locations of the fields over 2014-2016 which reflect farm management and other specifics like soil and environment at each site had significant effects on yield and NR when compared to a baseline ranging from − $91 to $352 ha −1 . Rates of N fertilizer applied under MZ were about 8% less than N applied under uniform rates as nutrient supply matched more effectively with crop demand. We attribute this modest N application between uniform and management zones to similar yield goals for different zones, and spatial distribution of MZs, where fertilizer rates in high zones likely offset low rates in low zones. The potential for spatial management of N fertilizer does exist; though the application of three management zones based on historical yield maps combined with soil test recommendations does not account for variability in all fields and is not consistently effective.

Materials and methods
Experimental fields. Research   . These ten sites were selected, for analysis of yield and net revenue, from a larger agronomic data set of 27 sites in the Canadian Prairies from a research project conducted from 2014 to 2017. This economic subset was selected based on the availability of suitable economic data for the fields, as sufficient cost of suitable production and agronomic management data were not available for the remainder of the sites. Producers set yield goals for each field that were adjusted for each zone based on historical data. Fields were selected from those volunteered by producers based on the availability of 3-5 years of yield maps, variable rate equipment and GPS equipped combines, and logistical constraints such as travel distance for sample collection. Terrain attributes were calculated for each field, and soil properties including pH, conductivity, soil organic carbon were sampled. Meteorological data were obtained from nearby Environment Canada meteorological stations, and weed and plant disease were measured in each field. Canola varieties were identified by most producers in the project. However, detailed spatial analysis of yield, terrain attributes, canola varieties, soil properties such as pH, conductivity and organic carbon, weeds and plant disease and meteorological data will be addressed in subsequent publications, in analyses which are outside the scope of the objectives for economic analysis. The low, average, and high yield zones were determined for all fields based on analysis of 3-5 years yield maps with K-means cluster spatial analysis of normalized yields in SMS advanced 34 . The cluster analysis separated the fields into zones with low, average and high normalized yield, based on historical data. The number of zones (3) was predetermined prior to cluster analysis, in order to provide sufficient space for plots, replicates and logistics for field-scale seeders and combines. Four nitrogen fertilizer treatments nested within each zone in each of four replicates (A, B, C and D) in each field (Fig. 1A). Nitrogen fertilizer treatments as urea or blends  Uniform  103  134  129 134  123 118 101 116  136  178  126   MZ  103  123  119 127  118 107 98  110  128  164  119   Difference  0  11  10  7  5  11  3  6  8  14  7 N (soil + applied) at maximum yield www.nature.com/scientificreports/ were based on soil tests for N (2 M KCl extract) and producer yield goals for each zone. Fertilizer rates were determined for soil samples collected in the fall or spring prior to seeding canola, from each plot and subplot in each field, and were based on recommendations in Manitoba 33 . These rates were applied in the spring based on the producer's seeding and production system for all canola varieties. Target yields were identified by producers for the field average based on their experience, and were closely related to those determined from the historical yield maps. Target yields were further adjusted for low and high zones, relative to the average. Phosphorus, potassium (K), and sulphate (S) fertilizer rates were based on soil test P, K and S, were uniformly applied in all treatments, and judged sufficient for the study. Mono-ammonium phosphorus (MAP) was applied based on soil test, to all treatments. The experimental design addresses the yield response and economics of canola for farm fields in response to two fixed factors, fertilizer rates (ordinal) nested within management zones (ordinal) based on soil test analysis, in a strip split plot design in linear models which specify replicates, treatments, fields and interactions. Experiments were conducted by producers at the field scale, canola varieties, equipment, agronomic management and fertilizer used for production varied between farms. All canola production data were collected with GPS equipped yield monitors and were calibrated on an annual basis with commercial scales or weigh wagons. Yield monitor variables, including site accuracy and outliers, and alignment of GPS points along harvest paths were reviewed for quality control. Yield monitor data, with values greater than the 97.5% quantile, and less than the 0.1% quantile, were removed. The remaining points were averaged by plot for statistical analysis of economic return, but retained for regression analysis of yield response to N fertilizer.
Economic analysis. Net revenue (NR), as calculated by Zentner et al. 36 , was defined as the revenue remaining after paying for all costs. This included variable costs (i.e. seed, nutrients, weed and disease control, transportation, fuel and oil, repairs, miscellaneous expenses, land taxes, interest cost on variable inputs), machinery and buildings land investment and ownership costs (depreciation, interest on investment, insurance and housing), and labour costs. Costs and NR were estimated using 2019 input prices (Table 1). For each site, the NR was averaged and expressed in CND $ ha −1 .
Initial crop budgets were determined using field trial agronomic data, and inputs and machinery price and cost data. All field operations from pre-seeding to harvest were included in the economic analysis. Input and output prices are summarized in Table 1. The net farm-gate price of $485 tonne −1 (t) for canola (i.e. the net cost of rail transportation and elevator handling) was used for canola gross income calculation 37 . The mean values for N and P (P 2 O 5 ) were $1240 t −1 and $1200 t −1 , respectively (Table 1).
We assumed that all farm machinery was operated by farmers and no custom work costs were included. Crop insurance, utilities and property taxes were not included in the total cost. For the effects of management zones and fertilizer treatments to be seen, participation in the Canada-provincial Crop Insurance Program was not included. No allowance was made for interest costs related to land equity. The skilled labour cost was assumed to be $25 h −1 for crop production. Farm practices followed by producers in western Canada, including equipment usage and field operation schedules, were established for each site. An average farm size of 907 ha was used to calculate work rates and farm machinery sizes 38 . For example, for fertilizer application, purchased price of machinery or equipment minus residual value was spread over optimal life of equipment and annual hours (h) of use (cost h −1 ) and investment and housing cost were added and then divided by work rate of machinery (hectares (ha) h −1 ) to calculate fixed cost for this operation per ha 38 . Fixed cost includes the cost of depreciation of the equipment due to use and years in service. It also includes an investment cost (e.g. the cost to borrow money to purchase the equipment and/or the lost interest revenue if that money had been invested), and housing/insurance costs 38 . Repair and maintenance, fuel and lubrication, and labour costs were also divided by work rate of the equipment to estimate each of these costs. Similar methods were used for each of the other field operations and added up to calculate field activity costs. Fertilizer costs, chemical costs, and seed costs were calculated by multiplying application rates and prices of each fertilizer, chemical, or seed. Combining and transportation costs from field to farm storage and from farm to first point of delivery were calculated in similar fashion assuming different truck sizes 38 . Grain truck capacity was assumed 9.5 t for field to yard and $9.5 t −1 for the first 5 km plus $0.92 t −1 for each additional kilometer (Table 1). Grain truck capacity was assumed 13 t for farm to first point of delivery with 100-km distance and 60 km h −1 speed. The grain truck was used 150 h per year. Transportation costs including fixed, repair and maintenance, fuel and lubrication, and labour costs based on crop yields were calculated similar to the other field operations. Costs and therefore NRs not only are affected by fertilizer treatments but also by the yield dependent costs as yield differences among zones or N treatments will affect transportation costs. Costs for auger and associated tractor use were also incorporated and it was assumed 20 min time was required to load and 5 min to unload a truck (Table 1). Interest rate on variable costs was assumed at 7%. This rate was used to reflect the range of public (5%) and private (10%) opportunity costs of capital. The analysis included the costs and returns for each site and year.
Yield data for 27 sites, from research conducted during 2014-2017 in Alberta, Saskatchewan and Manitoba, were related to replicate, yield zone, fertilizer treatment, and site-year in a linear model. Canola yield was analyzed for the period from 2014 to 2017 with a mixed model, where years, farms 39 www.nature.com/scientificreports/ treatments were defined as fixed effects and nested based on the experimental design, while experimental and observational error were expressed as random 40 . All interactions were specified in the model, with the exception of replicate. The model specified yield as the dependant continuous variable, with fertilizer treatment as ordinal, yield zone and site-year as nominal independent variables (fixed effects) and replicate as a nominal variable (random effect). Each field was considered as a site year in the linear model, which includes variability due to seasonal precipitation and temperature at that site. Canola yield was calculated from yield monitor data, with 4-7 subplots for each combination of replicate, yield zone, and fertilizer treatment in each field. Plot dimensions varied from field to field as influenced by seeder and swather width, and distribution of yield zones. Residuals were assessed for normality and homogeneity of variance, yield was log (base 10) transformed to meet these assumptions. Parameters for the linear model were calculated with maximum likelihood in a mixed model in statistical software JMP 35 . Denominator degrees of freedom for the linear model were calculated from subplots and independent variables with the Kenward-Roger first order approximation in JMP. Tukey's HSD was used to compare means following regression analysis. The economic analyses of site specific nitrogen management zones were based on both regression and statistical contrasts or analysis of variance (ANOVA). Estimated canola yield and NR were captured through change in nutrient, zones, and other variables such as farm management (site-year). Linear models of 2014-2016 data were specified and estimated as defined by Woolridge 41 . The regression model was specified as follows: where y ft consists of crop yield or NR for field f over time t, N is nitrogen available (soil plus applied N), a is the constant coefficient, b and c are coefficients of N, zone represents spatial variability, and site-year represents site and temporal trend effects of each field over 2014-2016. λ and θ are coefficients of zones and site-year variables. Regression analysis was done by introducing categorical variables for zones or site-year and quantitative values for N rates. Modeled yield, crop prices and crop costs, and potential NR changes were then calculated for canola 42 . The Least Squares regression between yield or NR and zone as well as N rate and site-year was conducted using Eviews Econometric Software 43 for each site-year or for all sites-years combined. The coefficients of the regression were tested and reported in Table 3.
Canola yield and NR were also analyzed with Proc Mixed 44,45 for each site-year or combination of sites-years, and the error term was used in subsequent comparisons of means with the Tukey method in Proc PLM. The management zone and N rate, expressed as a percentage of recommended fertilizer rates for the target yield in the field, were fixed effects; replicate was random. Interactions with field, zone and N rate were not significant and not included in the final analyses in the results and discussion. Residuals from the model were reviewed to confirm normality and homogeneity of variance. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.