Potato Yield, Net Revenue and Speciﬁc Gravity Responses to Nitrogen Fertilizer under Different Canadian Agroecozones

: Applying higher nitrogen (N) rates than required for optimum potato ( Solanum tuberosum L.) growth leads to economic and environmental losses. The extent to which the N rate associated with maximum potato yields differs from that maximizing net revenue (NR) or potato speciﬁc gravity is not fully understood. The objectives of this three-year study (2013–2015) conducted at ﬁve sites in three Canadian provinces (MB-1; MB-2; QC-1; QC-2; PEI) (15 site-years) were to: (i) assess potato marketable yield, NR, and speciﬁc gravity responses to increasing N application; (ii) calculate the N rate maximizing marketable (N max ) yield and NR using different statistical models. The year, N fertilizer, and their interaction were signiﬁcant on marketable yield and NR except at the MB-1 site where no signiﬁcant effect of N was observed. No signiﬁcant yield increases were observed at a N rate above 60 kg N ha − 1 at four site-years and above 120 kg N ha − 1 at ﬁve site-years, implying that the current recommended N rate could be reduced. All models ﬁtted the marketable and NR data equally based on R 2 , mean bias error or root mean square error and resulted in comparable predicted yield and NR values. However, N max values were different depending on the model with higher values being predicted by the quadratic- (161.4 to 191.9 kg N ha − 1 ) and the quadratic plateau models (60 to 191.9 kg N ha − 1 ), while lower N max values were obtained with linear plateau- (60.6 to 129.8 kg N ha − 1 ) and Mitscherlich–Baule plateau models (60.9 to 130. 9 kg N ha − 1 ). Nitrogen rate maximizing NR was on average 4% lower than the N rate maximizing marketable yields, except at one site where it was higher by 26 kg N ha − 1 when the quadratic plus plateau model was used. Speciﬁc gravity tended to decrease with the N rate. Our study conﬁrms trade-offs between the N rate maximizing yields or NR with that maximizing speciﬁc gravity. Nitrogen rate maximizing marketable yield and NR varies depending on the selected model. model. Speciﬁc gravity tended to decrease with increasing N fertilizer. Our results show trade-offs between fertilizing to maximize yield, NR, and speciﬁc gravity. Crop growth models that integrate pedo-climatic data could enhance the accuracy to estimate potato N requirements to identify the lower limit that allows maximizing NR and potato quality while limiting the environmental footprint.


Introduction
Potato (Solanum tuberosum L.) is the largest vegetable crop in Canada, with production reaching 4,770,521 tons in 2016 [1]. Worldwide, potatoes ranked as the third most consumed food, after rice (Oryza Sativa L.) and wheat (Triticum aestivum L.) [2]. Potatoes require a high amount of N due to their relatively shallow root system and low N use efficiency (NUE), which can be below 40% in humid conditions [3]. Furthermore, under cold and humid climatic conditions with a short growing season, high residual soil nitrates are encountered from the sensors were used to maintain soil moisture levels at approximately 65% to 75% available water capacity. The two soil types used in Manitoba were very different, which explains differences in the amount of water supplemented with irrigation. The MB-1 site was a fine sandy loam/sandy loam with rapid permeability and limited moisture-storage capacity in the upper (rooting) portion of the soil profile, while the MB-2 site was a clay loam texture with moderate permeability. Preceding crops prior to growing potatoes varied at each site and included canola (Brassica napus L.), corn (Zea Mays L.), barley (Hordeum vulgare L.), oat (Avena sativa L.), soybean (Glycine max L.), and wheat (Triticum aestivum L.) ( Table 1). Sensors were read two to three times weekly to make irrigation decisions. Average readings from the sensors were used to maintain soil moisture levels at approximately 65% to 75% available water capacity. The two soil types used in Manitoba were very different, which explains differences in the amount of water supplemented with irrigation. The MB-1 site was a fine sandy loam/sandy loam with rapid permeability and limited moisturestorage capacity in the upper (rooting) portion of the soil profile, while the MB-2 site was a clay loam texture with moderate permeability. Preceding crops prior to growing potatoes varied at each site and included canola (Brassica napus L.), corn (Zea Mays L.), barley (Hordeum vulgare L.), oat (Avena sativa L.), soybean (Glycine max L.), and wheat (Triticum aestivum L.) ( Table 1). Sensors were read two to three times weekly to make irrigation decisions. Average readings from the sensors were used to maintain soil moisture levels at approximately 65% to 75% available water capacity. The two soil types used in Manitoba were very different, which explains differences in the amount of water supplemented with irrigation. The MB-1 site was a fine sandy loam/sandy loam with rapid permeability and limited moisturestorage capacity in the upper (rooting) portion of the soil profile, while the MB-2 site was a clay loam texture with moderate permeability. Preceding crops prior to growing potatoes varied at each site and included canola (Brassica napus L.), corn (Zea Mays L.), barley (Hordeum vulgare L.), oat (Avena sativa L.), soybean (Glycine max L.), and wheat (Triticum aestivum L.) ( Table 1).

Soil Analyses Prior to the Trial Establishment
Prior to the trial establishments, soils (0-20 cm depth) were sampled in each of four blocks to characterize the sites. Soil samples were taken with a Dutch Auger at five to six randomly selected sampling points and were mixed to constitute one composite sample per each block. The soil series in each province were representative of those normally used to grow potatoes. Soil textural groups were loam, loam sand, silt loam or sandy loam, with the MB-2 and QC-1 sites having the highest sand content (Table 2). Soil total organic carbon (TOC) and total N were determined using the dry combustion method on an Elementar analyzer (Vario Max, Elementar Analyzer, Hanau, Germany). The soil pH was measured in water (1:1 ratio) [33]. The cation exchange capacity (CEC) was determined using the ammonium acetate method [34]. The soil bulk density was measured using the core method [35] in three replicates at each of the four blocks. Macro-nutrients (P, K, Ca, Mg) were extracted using Mehlich-3 solution [36].

Initial Soil Properties
The sites were different in terms of TOC, with values ranging from 11.2 to 27.1 g kg −1 in 2013, from 11.3 to 25.9 g kg −1 in 2014, and from 15.4 to 29.3 g kg −1 in 2015 (Table 2)

Nitrogen Fertilizer Treatments
At each site, five increasing N rates were applied (0, 60, 120, 180, 240 kg N ha −1 , referred as N0, N60, N120, N180, and N240, respectively). Sources of N fertilizer and placements were different based on local common practice, and fertilizer was broadcast in MB-1 and MB-2 and banded in QC-1, QC-2 and PEI (Table 1). Split N application was performed at sites located in MB-1 and MB-2 (as urea) and QC-1 and QC-2 (as a combination of ammonium sulfate and calcium ammonium nitrate), whereas all N fertilizer was applied once (as ammonium nitrate) at planting in PEI (Table 1). Other nutrients and agricultural practices to control weeds and pests were applied to ensure optimum crop growth [37,38].

Temperature and Moisture Trend in Comparison to 30-Year Average
At all sites, the temperature was higher than the 30-yr average (Table 3). Total precipitation was higher than the 30-yr average at sites located in MB, whereas it was lower than the 30-yr average at sites located in QC and PEI. At irrigated sites in MB, the irrigation represented between 26% and 54% of total water input, depending on the year. Climate was drier on sites located in Manitoba than on sites located in Quebec and PEI based on the 30-yr average precipitation (Table 3). Table 3. Total amount of water input from planting to harvest and mean temperature during the growing seasons (May-September) measured at different sites in comparison with the 30-yr normal.

Potato Yield Calculation
The potato variety was Russet Burbank for all sites. Total yield was calculated by harvesting two middle rows 5 m long. Potato marketable yields were calculated based on local guidelines. Tuber size categories included Culls: <38 mm; Can1-small: 38 to 51 mm; Can1: 51 to 89 mm; Can1-large: 89 to 114 mm; Jumbo: >114 mm. In Manitoba and Quebec, marketable yield included Can1 and Can1-large, while it included Can1-small, Can1, and Can1-large for PEI. Finally, potato specific gravity was determined on a subsample of 30 tubers from marketable potatoes per plot as "[weight in air/(weight in air-weight in water)]" [39].

Net Revenue Calculation and Statistical Analyses
Agronomic data from the field trials, combined with price and cost data for machinery and other inputs, were used to develop potato budgets, expressed in CAD ha −1 . The economic analysis included all inputs used in field operations from pre-seeding to harvest activities. Total costs included all monetary costs (i.e., seed, fertilizer, chemicals, fuel, and oil, repairs, irrigation if applied, transportation, miscellaneous expenses, land taxes, and interest cost on variable inputs), ownership costs on machinery and buildings (depreciation, interest on investment, insurance, and housing), and labor. No allowance was made for interest costs related to land equity. Labor cost was assumed to be CAD 25 h −1 . Costs for field irrigation schedules and equipment used were determined for irrigated potato production in Manitoba based on 2018 guidelines for estimating irrigated processing potato costs in Manitoba [40]. Enterprise budgeting technique was used assuming a similar base cost for potato production in all sites, then cost differences based on specific fertilizer application, field activities, and marketable potato yield at each site were accounted for. Table 4 provides a summary of different fertilizer and potato prices in Manitoba, Quebec, and PEI. Fertilizer prices from various sources were used to estimate the costs of potato production. Three total costs were calculated based on three scenarios of fertilizer prices as explained in Table 5. Similarly, potato marketable tuber prices from various sources were used to represent the value of potato [41]. Three gross revenues were calculated based on three scenarios of potato prices (Table 5). NR was defined as the gross revenue remaining after paying for total costs as described by Zentner et al. [42] and Khakbazan et al. [43]. Nine combinations of NRs based on three potato and three fertilizer price scenarios were calculated (Table 5). Since the results for all scenarios were the same, only results from scenario net12 were used (net12, Table 5).  The MIXED procedure of SAS [44] was used to analyze the effect of N rates on marketable yield, NR, and specific gravity. Block was treated as a random factor, and treatments were treated as fixed effects. Because the nutrient recommendations are based on multiple year studies, the year was treated as a random factor to assess potato response to N fertilizer independently to yearly variations and, thus, allow generalization of potato N response at a specific site. On the other hand, annual variations are also important, and the year should be treated as a fixed factor. To capture both sets of information, an analysis of variance was performed per site, in which year was first treated as a fixed factor and then as a random factor. This allowed identification of sites where the interaction between the N rate and the year was significant and whether the N rate was still significant when data was combined over 3 years per site. The protected least significant difference (LSD) multiple comparison test was carried out to compare means after significant effects (5% probability level) were found in the ANOVA analysis.
To identify N rate maximizing marketable yield and NR, and corresponding maximum yield and NR, marketable yield and NR were fitted with five different models (quadratic, quadratic plateau, linear plateau, Mitscherlich-Baule, Mitscherlich-Baule plateau) using raw data. Coefficients for the selected models were calculated using MATLAB R2019a statistics and machine learning toolbox and are described by Equations (1)-(7) below.
The linear plateau model is described by Equations (1) and (2) where Y is the potato yield (Mg ha −1 ) or NR (CAD ha −1 ), a is the intercept, b is the linear coefficient and x is the rate of N application (kg N ha −1 ), C is the critical rate of fertilization where Y is the potato yield (Mg ha −1 ) or NR (CAD ha −1 ) and x is the rate of N application (kg N ha −1 ), a, b, and c are the coefficients for quadratic function. To obtain optimum N rate maximizing yield and NR, the first derivative of the response Equations was set to zero and then was solved for x. The quadratic plateau model is described by Equations (4) and (5) where C is the critical rate of fertilization that occurs at the intersection of the quadratic response and the plateau line. Note that the right side of Equation (5) is a constant value and represents plateau yield.
The Mitscherlich-Baule model where Y is the marketable yield in Mg ha −1 or NR (CAD ha −1 ); Y max,MB . is the maximum yield or maximum NR under the Mitscherlich-Baule function; a N is a coefficient directly affecting the slope of MB function. N in is a coefficient derived from curve fitting and represents the initial amount of nitrogen in the soil (kg N ha −1 ). Mitscherlich-Baule plateau: where Y max,MB−P is a coefficient representing the maximum achievable yield or NR under the Mitscherlich-Baule plateau function if there was no constraint on maximum yield imposed by the plateau function, and b N is a coefficient directly affecting the slope of the Mitscherlich-Baule plateau function. In this case, the maximum yield or NR is calculated by substituting the C into the Mitscherlish-Baule function. Two errors were calculated for all models. Root mean square error (RMSE) representing the noise in data, and the mean bias error (MBE) representing the average bias in prediction. The units for both errors are the same as the unit for dependent variable. A high RMSE value indicates the data are highly scattered around the predicted curve. A high positive or negative MBE indicates predictions are systematically higher or lower than measurements.
Nitrogen rate maximizing NR can be comparable to the economic optimum N rate (EONR, kg N ha −1 ), which is normally calculated as the application of N where CAD 1.00 of extra N fertilizer returns CAD 1.00 of potatoes [45]. The latter assumes that fertilizer N cost is the only variable cost and that all other costs are fixed. The EONR is calculated by setting the first derivative of the Y response curve equal to the ratio between the cost of N fertilizer (CAD kg N −1 ) over the price of potatoes (CAD Mg −1 tuber). Calculating EONR using the latter approach or calculating N rate maximizing NR yielded comparable results because, in both cases, the main difference across treatments at a specific site remains the fertilizer Agronomy 2021, 11, 1392 9 of 21 price. Net revenue not only considers N fertilizer cost but also integrates all costs at the farm level, including additional yield dependent transportation costs. Given that we had an application of two N sources in Quebec, it is preferable to calculate N rate maximizing NR instead of EONR in order to assess how NR varies across different sites where potato management slightly differs in terms of N fertilizer type and irrigation regimes. It can be justifiable to omit zero N rate and the lowest N rate, as it is expected that the optimum economic N rate is not reached at low rate [46], but it was reported that omitting the zero N rate resulted in increased confidence intervals and the parabola vertices lying either far left for the convex parabola, or far right for the concave parabola [46]. Therefore, our analysis included all N rates, including the check (N0) treatment.

Effect of Year and N Rate on Marketable Potato Yield
At the MB-1 site, only year effect was significant when year was treated as a fixed effect (Table 6). When year was treated as a random effect, N rate effect was not significant (Table 6). Significantly lower marketable yields were observed in 2014, but yields were comparable between 2013 and 2015 (Table 7). Mean marketable yield averaged across 3 years ranged from 35.6 to 40.2 Mg ha −1 (Figure 1a). At the MB-2 site, when year was treated as a fixed effect, the year, N rate, and their interaction were significant (Table 6). By slicing the interaction, no significant effect of N rates was observed in 2013, but N rate effect was significant in 2014 and 2015, with no significant yield increases observed above N180 (Table 7). When year was treated as a random effect, N rate was significant (Table 6). Overall, marketable yield means ranged from 30.3 to 42.1 Mg ha −1 (Figure 1a). At the QC-1 site, when year was treated as a fixed effect, the year, N rate, and their interaction were significant on marketable yield (Table 6). No significant yield increases were observed at N rate above N180 in 2013, whereas no significant yield increases were observed at N rate above N120 in 2014 and 2015 (Table 7). Nitrogen rate was significant when year was treated as a random factor. Overall mean values ranged from 10.8 to 31.0 Mg ha −1 (Figure 1a). At QC-2, when year was treated as a fixed factor, the year, N rate, and their interaction were significant (Table 6). By slicing the interaction, N rate was significant in all years, with no significant yield increases observed at N rate above N60 in 2013 and 2014 and no significant yield increases at N rate above N120 in 2015 (Table 7). Nitrogen rate was also significant when year was treated as a random factor (Table 6). Overall, mean values ranged from 17.4 to 30.8 Mg ha −1 (Figure 1a). At the PEI site, the year, N rate, and their interaction were significant when year was treated as a fixed effect (Table 6). No significant yield increases were observed at N rates beyond N60 in 2013 and 2014 (Table 7). Nitrogen rate was significant when year was treated as a random effect (Table 6). Overall, means ranged from 33.8 to 42.3 Mg ha −1 (Figure 1a). Table 6. Analysis of variance of the effects of N rate on marketable yield, specific gravity and NR, treating year as a fixed or random factor.

Year As a Fixed Factor
Year As a Random Factor

Effect of Year and N Rate on Net Revenue
Similarly to what was observed for marketable yields, when year was treated as a fixed effect, the effects of year, N rate, and their interaction were significant on NR at all sites except at the MB-1 site where only year effect was significant (Table 6). When year was treated as a random factor, N rate was significant at all sites except at the MB-1 site (Table 6). Based on the overall means across years, the highest NR was associated with N60 (2200 CAD ha −1 ) at MB-1, with N180 (2457 CAD ha −1 ) at MB-2 site, then N120 (2226 CAD ha −1 ) at QC-1, N180 (2199 CAD ha −1 ) at QC-2, and N180 (2936 CAD ha −1 ) at PEI (Figure 1b). Yearly variations were observed at all sites. At the MB-1 site, negative NR values were observed in 2014 with all N rates. In 2013, the highest NR was observed at N60 (3984 CAD ha −1 ) and at N120 (3524 CAD ha −1 ) in 2015 (Figure 2a). At the MB-2 site, the highest NR (2060 CAD ha −1 ) was observed at N180 in 2013, with N120 in 2014 (2305 CAD ha −1 ) and N180 in 2015 (3098 CAD ha −1 ) (Figure 2b). At the QC-1 site, negative NR was observed in 2013 with all N rates except at N240 (491 CAD ha −1 ) and the highest NR was observed at N120 in 2014 (5411 CAD ha −1 ) and N240 in 2015 (2431 CAD ha −1 ) (Figure 2c). At the QC-2 site, negative NR was observed in 2014 independent of the N rate and the highest NR was associated with N180 in 2013 (2985 CAD ha −1 ) and with N240 in 2015 (7740 CAD ha −1 ) (Figure 2d). At the PEI site, the highest NR was associated with N60 in 2013 (3494 CAD ha −1 ), with N120 in 2014 (1326 CAD ha −1 ), and with N180 in 2015 (6231 CAD ha −1 ) (Figure 2e).

Effect of Year and N Rate on Specific Gravity
At the MB-1 site, year and N rate were significant, while their interaction was not significant ( Table 6). Close values were observed in 2013 and 2015 and higher numerical values were observed in 2014, with a trend towards high value with the check treatment (Figure 3a). When year was treated as a random factor, N rate effect was significant (Table 6). Specific gravity values ranged from 1.083 to 1.094 (Figure 3a), with grand mean values ranging from 1.086 to 1.092 (Figure 4a). At the MB-2 site, year and N rate effects were significant, and their interaction was not significant (Table 6). Close values were observed in 2013 and 2015, and a high numerical value was observed in 2014 (Figure 3b). When overall means were compared, lowest value was observed at N240 and highest value at N60, but this was comparable with N120 ( Figure 4b). Nitrogen rate was significant when year was treated as a random effect (Table 6). Specific gravity values ranged from 1.088 to 1.099 (Figure 3b), with grand mean values ranging from1.090 to 1.094 (Figure 4b). At the QC-1 site, year, N rate, and their interaction were not significant. Specific gravity values ranged from 1.075 to 1.089 (Figure 3c), while grand mean values ranged from 1.080 to 1.083 (Figure 4c). At QC-2, N rate effect was significant when year was treated as a fixed effect (Table 6), with values decreasing at N rate above N60 (Figure 3d). Nitrogen rate was not significant when year was treated as random effect (Table 6). Values ranged from 1.073 to 1.097 (Figure 3d), with grand mean values ranging from 1.085 to 1.093 (Figure 4d). At the PEI site, year, and the interaction between year and N rate were significant (Table 6). When the interaction was sliced, N fertilizer effect was only significant in 2013, with values significantly declining at N rate above N60 (Figure 3e). Nitrogen rate was not significant when year was treated as a random effect (Table 6)

Effect of Year and N Rate on Specific Gravity
At the MB-1 site, year and N rate were significant, while their interaction was not significant (Table 6). Close values were observed in 2013 and 2015 and higher numerical values were observed in 2014, with a trend towards high value with the check treatment (Figure 3a). When year was treated as a random factor, N rate effect was significant (Table  6). Specific gravity values ranged from 1.083 to 1.094 (Figure 3a), with grand mean values ranging from 1.086 to 1.092 (Figure 4a). At the MB-2 site, year and N rate effects were significant, and their interaction was not significant (Table 6). Close values were observed in 2013 and 2015, and a high numerical value was observed in 2014 (Figure 3b). When overall means were compared, lowest value was observed at N240 and highest value at cant when year was treated as a random effect (Table 6). Values ranged from (Figure 3e), with grand mean values ranging from 1.083 to 1.095 (Figure 4e)

Nitrogen Rate Maximizing Marketable Yield and Net Revenue Using Different Models
Different models predicted comparable maximum yields, and all five models seeme to fit data equally based on R 2 , MBE, or RMSE (Table 8). In general, the highest Nmax wer associated with the quadratic model (from 161.4 to 191.9 kg N ha −1 ) followed by the quad ratic plateau (60.7 to 191.9 kg N ha −1 ), and the lowest values were associated with the linea plus plateau model (from 60.6 to 129.8 kg N ha −1 ) and the Mitscherlich-Baule platea model (60.9 to 130.9 kg N ha −1 ) ( Table 8).

Nitrogen Rate Maximizing Marketable Yield and Net Revenue Using Different Models
Different models predicted comparable maximum yields, and all five models seemed to fit data equally based on R 2 , MBE, or RMSE (Table 8). In general, the highest N max were associated with the quadratic model (from 161.4 to 191.9 kg N ha −1 ) followed by the quadratic plateau (60.7 to 191.9 kg N ha −1 ), and the lowest values were associated with the linear plus plateau model (from 60.6 to 129.8 kg N ha −1 ) and the Mitscherlich-Baule plateau model (60.9 to 130.9 kg N ha −1 ) ( Table 8). Nitrogen rates maximizing NR were very close to those maximizing marketable yields (Table 9) and were on average 4% lower than those predicted using marketable maximizing marketable yield when data was fitted with the quadratic plus plateau model (Table 8 vs.  Table 9).

Marketable Yield Response to Nitrogen Rate
Omitting lower yields observed at the QC-2 site in 2014, which were most likely due to drought, marketable yield ranges observed at different sites are comparable to those reported in previous studies in Eastern Canada. Yield ranges are comparable to those reported in New Brunswick, Canada by Nyiraneza et al. [47] with yields up to 35 Mg ha −1 , by Zebarth et al. [48] with yields ranging from 13 to 49 Mg ha −1 , and by Snowdon et al. [49] who reported average values of 34.4 Mg ha −1 . Our values are also close to those reported by Cambouris et al. [20] in a study conducted in Quebec with yields ranging from 13 to 38 Mg ha −1 .
All sites responded to N application, except the MB-1 site. This site was characterized by high TOC, ranging between 25 and 29 g kg −1 , and the highest CEC, which was above 30 meq 100 g −1 . The lack of responsiveness to N fertilizer at the MB-1 site may be explained by high fertility status in combination with supplemental irrigation, which may have enhanced the soil N supplying capacity. Total organic carbon at QC sites was comparable to that at the MB-1 site but yield response to increasing N rates was observed at the former site. Lower yields at sites located in Quebec than in Manitoba may be attributed to lower soil pH and a lack of irrigation. The soil pH at sites in Quebec was below 5.5, except the QC-2 site, which had a soil pH of 5.9 in 2013. A soil pH below 5.5 can restrict P availability in podzolic soils where aluminum tends to bind with P, decreasing its availability [50].
Yield associated with the check treatment represented 89%, 72%, 35%, 56%, and 76% of the highest yield obtained at MB-1, MB-2, QC-1, QC-2, and PEI sites, respectively. Therefore, the sites located in Quebec responded the most to N applications but were also associated with lower yields than other sites. Overall, no effect of N rate was observed at the MB-1 site and MB-2 site in 2013. In four site-years (QC-2-2013, 2014; PEI-2013, 2014), applying N60 resulted in comparable yields as when applying N240 (Table 7). Comparable yield between N120 and N240 was also observed in five site-years (MB-2-2014; QC-1-2014, 2015; QC-2-2015; PEI-2015). This corroborates previous findings that the fertilizer N rate could be reduced below the recommended N rate without affecting yield. A 3-yr study in Manitoba reported an absence of N fertilizer response when it was applied at 80% of the recommended rate and concluded that current N recommendations over-estimate sitespecific potato N requirements [17]. In Manitoba, N rate recommendations under irrigated land ranged from 0 to 291 kg N ha −1 [51] and changed under dry land and irrigated land and depended on yield goals. In PEI, the recommended potato N fertilizer for Russet Burbank was 185 kg N ha −1 [52], whereas in Quebec it varied between 120 kg N ha and 170 kg N ha −1 [53]. A 5-yr study that compared two potato varieties' (Russet Burbank and Kennebec) response to increasing N rates in PEI reported no yield advantage at N rates above 134 kg N ha −1 , and decreased specific gravity was observed with increasing N rates [16]. Nitrogen rates ranging from 90 to 130 kg N ha −1 were reported to be optimal for processing potato varieties in Atlantic Canada [13][14][15]. Another study that assessed Russet Burbank potato response to increasing N fertilizer in Michigan, reported that N rates used by most growers could be reduced to 112 kg N ha −1 without significant effect on tuber yields [12]. Conversely, in another study conducted in Michigan, tuber yield increases were reported with N application ranging from 225 to 270 kg N ha −1 [54,55]. A study in Maine tested increasing N rates up to 270 kg N ha −1 and reported an optimum N rate of 126 kg N ha −1 and 136 kg N ha −1 for Russet Burbank and Shepody cultivars, respectively [9].

Nitrogen Rate Maximizing Marketable Yield Versus That Maximizing Net Revenue
Considering the main goal of a producer is to maximize profit, crop response to increasing N rates should not solely be interpreted based on statistical significance or highest yield values, and it is thus important to identify the N rate that maximizes NR. Nitrogen rate maximizing marketable yield and NR varied depending on the model used, and N rates maximizing NR were slightly lower than those maximizing marketable yield. This corroborates previous findings [28,29] that the N rate maximizing yield may be different than that maximizing NR.
Nitrogen rate maximizing yield and NR varied depending on the model. The predicted maximum yield provided little basis for selecting one model over another. The reason for selecting one model over others influences the recommended N rate. Other previous studies have reported N maximizing yield or economic return to vary depending on the models used [20,[56][57][58][59]. Corn response to N fertilizer was assessed, and linear plateau, quadratic plateau, exponential, and square root models were used to fit the data to calculate the economic optimum N rate and it was reported that predicted maximum yields were comparable, but different economic optimum N rates were obtained [56]. The latter authors reported that R 2 was not a satisfactory goodness of fit, similarly to this study, where R 2 were comparable independent of the model used. Low R 2 values were obtained when all data were pooled together (n = 60) for each site, reflecting the annual variations (Tables 8 and 9), but R 2 values increased when overall means were used (i.e., quadratic model Figure 2). Our results are aligned with those of Cambouris et al. [20], who reported higher N rate maximizing yield with the quadratic model and quadratic plus model than with the linear plateau model.
By fitting NR raw data with the same models used to fit marketable yield data, we again observed that predicted maximum NR values were closely independent of the model used but N rate maximizing NR values were different, with the quadratic and quadratic plus models being associated with higher N max . Predicted maximum NR ranged from 1875 to 3066 CAD ha −1 . Our results confirm that N rate maximizing marketable yield or NR fall within a very wide range. It was reported that the economic profit would not change significantly for an N application ranging from 150 kg N ha −1 to 250 kg N ha −1 [46], but excess N fertilizer is expected to be associated with an environmental cost.
Models have their own weaknesses and may fail to describe the reality accurately. For instance, the quadratic model, which is frequently used to fit yield response to nutrient applications, assumes that yield will increase up to a certain maximum and then decline with additional N application and assumes a symmetry between yield increases below the optimal and yield decreases above the optimal [60]. However, it is very rare to find symmetry in nature. A strong decline in potato yield with increasing N rate is not always observed [61]. In addition, the constant slope of the linear plateau model is not natural [46]. According to Yang et al. [60], uncertainty in precisely identifying optimum N rates is due to the fact that statistically significant trends are expected only when N rates are below the optimal N rate and that yield responses are indirect estimations of N sufficiency. This study corroborates previous studies concluding that tools are still needed to precisely identify optimum N rates.
According to Cruz-Medina et al. [58], accurate discrimination between the models would require experiments with more levels of fertilizer (increasing N rates with small increments) and including N rates higher than normally used. Another alternative for model selection is using simulated yield response to small N increments obtained from processbased crop models with detailed processes on the N cycle within the soil-crop-atmosphere system. A previous study used this approach and simulated corn in the predominant Canadian corn production ecozone for 30 plus years [62]. The N responses obtained from simulations were then used to evaluate the performance of different models (Mitscherlich-Baule and Mitscherlich-Baule-plateau) in predicting the slope of yield functions. The simulations also indicated significant differences in corn yield response to N rates from year to year. Potato N requirements also change from year to year, and they are expected to be lower or higher than what is recommended depending on the growing season conditions. We observed yield fluctuations from one year to another even at the irrigated site. Nitrogen can only be used efficiently when soil moisture is not limiting. Applying most of the N fertilizer at planting can reduce N uptake and recovery [29,30]. In Manitoba, N recommendations are based on fall soil nitrate contents [51]. In QC and PEI Canada, due to wet falls and springs characterized by a high risk of nitrate leaching, estimating N supplied by the soil at the beginning of the season is still challenging. Bélanger et al. [45] reported that, under wet conditions, pre-plant soil nitrate was not a good predictor of potato N needs in Atlantic Canada and reported that mid-season soil nitrate was a better index to estimate potato N requirement. Remediating N deficiency in the middle of growing may be challenging due to the waiting time needed to receive results from the laboratory. Models need to take into account crop growth stage, environmental conditions, and N sources and sinks; smart decision-support systems [61], precision agriculture and site-specific management to fine-tune N applications that meet a crop's N requirements.

Potato Specific Gravity Response to N Fertilization
At four sites out of five, specific gravity tended to decrease with N rate. It was previously reported that excess N fertilization was associated with lower specific gravity [9,32]. In cold climates, indeterminate potato cultivars such as Russet Burbank may require lower N, as the short growing season does not allow the late varieties to reach their potential. Reduction in tuber maturity is often associated with lower specific gravity. In addition to decreasing specific gravity with increasing N rate, higher residual soil N was reported at higher N rates, and these were not translated into increased potato yield [6]. Our study demonstrated trade-offs between fertilizing to achieve highest potato yield, NR, and specific gravity. The challenge is to identify the lower limit of optimum N rate to minimize the environmental footprint while maximizing NR.

Conclusions
This study investigated potato marketable yield, NR, and specific gravity responses to increasing N rates at 15 site-years under different Canadian agroecozones. The study also determined N rate maximizing marketable potato yield and NR using different models. At one site, marketable yield did not respond to N applications. At the rest of the sites, no yield advantages were obtained at N rates lower than currently recommended, implying the possibility of reducing N rate without impacting yields. Nitrogen rates maximizing marketable, and NR were close and varied depending on the models used, with the quadratic model predicting higher N rates than the linear plus plateau model and the Mitscherlich-Baule plateau model. Current recommended N rates at each site fall within the range of estimated N maximizing marketable yield and NR. Our results confirm that N maximizing yield and NR falls within a large range as it varies depending on the selected model. Specific gravity tended to decrease with increasing N fertilizer. Our results show trade-offs between fertilizing to maximize yield, NR, and specific gravity. Crop growth models that integrate pedo-climatic data could enhance the accuracy to estimate potato N requirements to identify the lower limit that allows maximizing NR and potato quality while limiting the environmental footprint.
Author Contributions: J.N. carried out manuscript writing, data analysis, and figure and table preparation; J.N., A.N.C., A.N. and J.L. were site leads of the trials conducted at different sites and coordinated the field activities and data collection; M.K. calculated net revenue; M.M. contributed to statistical analyses by calculating N max for marketable yield and NR using different models; I.P. compiled initial data from all sites; N.Z. coordinated soil analyses to characterize each site. All authors have read and edited the first version of the manuscript and agreed to the submitted version of the manuscript.
Funding: This research was funded by Agriculture and Agri-Food Canada through an A-Base program (Project J-000106, Development of decision-support components for optimal N fertilizer applications at the field scale).