Interplay of urbanization and agricultural modernization shapes nitrogen use in global croplands

Domain: neuroscience
Published: 2026-03-27
Editorial outcome: minor_revision

Abstract
Urbanization reshapes agricultural systems through labor and land-use changes, interacting with modernization processes including farm size expansion, mechanization, and irrigation to drive nonlinear trends in cropland nitrogen use. Using a 61-year dataset from 139 countries, here we show that the association between urbanization and nitrogen outcomes is profoundly nonlinear and contingent on development stages. In low-income countries, urbanization initially increases fertilizer use while suppressing nitrogen yield and efficiency, though larger farm sizes mitigate these early losses. As countries reach upper-middle-income levels, modernization enhances nitrogen efficiency but introduces trade-offs between environmental gains and yield growth. In high-income countries, advanced modernization mitigates adverse urban impacts, reversing nitrogen use efficiency from a 4% decline to a 12% gain at high urbanization levels. These findings indicate that there is no universal sustainability pathway. Instead, integrating land consolidation, mechanization, and precision irrigation can transform urbanization into a catalyst for sustainable management and resilient food systems.

Introduction
Urbanization, one of the most significant global trends of the 21st century, profoundly reshapes population distribution and land-use patterns, with implications for agricultural systems and global food security1-3. A critical nexus in this transformation is the agricultural nitrogen cycle. While nitrogen-driven agricultural intensification has been instrumental in global food security 4, 5, its inefficient management has triggered severe, cascading environmental challenges6-8. Consequently, understanding the drivers that shape nitrogen management practices is paramount for navigating the trade-offs between productivity and sustainability in an increasingly urbanized world.

Existing literature has extensively explored the multifaceted associations of urbanization on agriculture with a dualistic view. One highlights the challenges, such as the irreversible loss of cropland and disruptions to rural labor markets2, 3, 9-11. Another suggests that urbanization can be a catalyst for agricultural modernization by fostering technology adoption and farm consolidation12-14. While insightful, these studies often focus on estimating a net average effect, potentially overlooking that the association of urbanization with agricultural outcomes is neither linear nor uniform. A critical knowledge gap exists in understanding this conditional heterogeneity: how does the association evolve with socioeconomic development in a country, and how is this evolution linked to the pace of its own agricultural modernization? Specifically, as urbanization coincides with rural labor migration and land consolidation, how do the resulting changes in farm size, machinery intensity, and irrigation intensity act as moderating variables that can either amplify or buffer the outcomes associated with urbanization at different stages? Answering these questions requires moving beyond a simple net-effect estimation to a more dynamic, interactive framework that can systematically map these nonlinear and context-dependent pathways.

To fill this gap, the central conceptual innovation of this study lies in integrating urbanization, agricultural modernization, and nitrogen trends into a unified, nonlinear moderating framework. Prior research has largely examined these domains in isolation. While linear global models are insightful, they often focus on estimating a net average effect, thereby masking critical threshold effects and structural phase transitions. Conversely, regional case studies offer rich granular details but often lack the macro-level perspective necessary to generalize how these relationships evolve across distinct stages of economic development. By capturing these joint nonlinearities, our approach yields macro-systemic insights that neither linear models nor regional studies can provide alone.

We hypothesize that the association of urbanization with nitrogen fertilizer application, nitrogen yield, and nitrogen use efficiency (NUE) is not uniform, but is profoundly nonlinear and critically contingent upon the stage of economic development and the prevailing levels of farm size, machinery intensity, and irrigation intensity. In low-income countries, the primary challenge is whether farm size can grow fast enough to compensate for labor loss and create a viable scale for machinery intensity15, 16. In contrast, in high-income countries, the challenges shift to managing the trade-offs of a highly optimized system , questioning whether further advances in machinery intensity and irrigation intensity can provide resilience against diminishing returns17. To test these hypotheses, we compiled a comprehensive panel dataset spanning 139 countries from 1961 to 2021. Recognizing the heterogeneity described above, our analytical approach moves beyond a single global model. We employ a dynamic, nonlinear analytical framework based on two-way fixed-effects models, a strategy specifically designed to robustly control for temporal dynamics and unobserved country-specific heterogeneity. This study provides a systematic, data-driven assessment of how the interplay between urbanization and agricultural modernization, contingent on development stage, jointly shapes the sustainability of global food production.

Results
Between 1961 and 2021, global patterns of urbanization and nitrogen use underwent profound transformations. The proportion of the global population living in urban areas rose from 31% to 56% (Fig. 1a), reflecting a worldwide transition from rural to urban societies. Although this urban shift occurred across all income groups, it was particularly pronounced in developing regions that started with low urbanization levels in 1961. East and Southeast Asia, Sub-Saharan Africa, West Asia, and Latin America experienced the most rapid increases (Fig. 1b–d), driven largely by demographic expansion in these regions18. In parallel, key indicators of agricultural modernization, farm size, machinery intensity, and irrigation intensity, also underwent significant and spatially heterogeneous changes (Fig. S1).

Over the past six decades, nitrogen fertilizer application on global cropland has risen dramatically. The global mean application rate increased from 9 kg N ha⁻¹ in 1961 to 67 kg N ha⁻¹ in 2021 (Fig. 1e). Most regions followed a steady upward trajectory, reflecting the broad diffusion of synthetic fertilizers after World War II. A notable deviation occurred in Eastern Europe, where nitrogen fertilizer application fell sharply in the 1990s after the dissolution of the Soviet Union (Fig. 1e, h). Despite favorable agricultural conditions provided by extensive black soils and high natural fertility19, the collapse of state-subsidized systems, rising input costs, and disrupted production and supply chains during the economic transition led to a substantial reduction in nitrogen fertilizer application20, 21. In contrast, China and India showed substantial increases in nitrogen fertilizer application, paralleling rapid agricultural intensification and sustained gains in crop production (Fig. 1f–h).

Global crop nitrogen yield, defined as nitrogen harvested per hectare based on model estimates, rose from 15 kg N ha⁻¹ in 1961 to 46 kg N ha⁻¹ in 2021 (Fig. 1i). This increase reflects advances in crop breeding22, improved management practices23, and the expanded use of nitrogen fertilizers5. The most substantial gains occurred in China, India, Europe, and North America, where intensive farming and technological innovations supported sustained productivity growth24. In contrast, Central Asia, particularly countries such as Kazakhstan and Uzbekistan, experienced stagnation and even sharp declines in nitrogen yield during the 1990s (Fig. 1j–l). The breakdown of Soviet-era agricultural systems, difficulties in transitioning to market economies, and restricted access to essential inputs severely constrained productivity21. Environmental stressors such as arid climates, soil degradation, and long-standing irrigation challenges further amplified these declines25, 26.

Despite the global surge in nitrogen inputs and yields, nitrogen use efficiency (NUE) declined markedly from roughly 45% in the 1960s to about 35% in the 1990s (Fig. 1m). This decline indicates that nitrogen fertilizer application growth far exceeded nitrogen yield improvements, producing diminishing returns27. Since the 1990s, however, NUE has gradually recovered to around 45% in the 2020s, driven largely by improved fertilizer management and more efficient production systems, especially in developed countries28. Some nations achieved higher nitrogen yields while reducing nitrogen fertilizer application. For example, Denmark29 and the Netherlands30 have substantially improved NUE through precision farming and optimized nutrient management28. Conversely, many rapidly industrializing countries such as China31 and India32 experienced declining NUE during periods of input-driven intensification, as nitrogen fertilizer application expanded faster than nitrogen yield gains. These contrasting trajectories underscore the intertwined evolution of urbanization, nitrogen use, and nitrogen losses across regions (Fig. 1n–p).

To elucidate the mechanisms underlying the macro-level trends, we employ a path analysis framework grounded in dynamic panel models (see Methods for detailed specifications). In this context, “average transmission pathways” refer to the baseline, unmoderated direct and indirect effects of urbanization on the nitrogen cycle. This approach enables the estimation of both direct and indirect pathways through which urbanization and agricultural modernization jointly influence the agricultural nitrogen cycle, providing a baseline for understanding the more complex nonlinear dynamics.

The direct association between urbanization and nitrogen yield and use efficiency evolves with economic development. In low- and lower-middle-income countries, a one percentage point increase in the urbanization rate corresponds to an average 0.67% increase in nitrogen fertilizer application, while its average direct associations with nitrogen yield and NUE remain statistically insignificant (Fig. 2a). This pattern suggests that, in early development, urbanization primarily drives input intensification rather than efficiency gains. As countries advance to upper-middle-income status, urbanization's average direct effects on nitrogen fertilizer, yield, and NUE become non-significant. While the former two lack a clear association, the neutrality for yield is due to a significant U-shaped curve (quadratic term: 0.97, p<0.05) where negative and positive effects offset each other (Fig. 2b). In high-income countries, the direct effects of urbanization on nitrogen yield and NUE also remain insignificant overall, but both exhibit significant U-shaped relationships, supported by positive quadratic coefficients for nitrogen yield (1.45, p < 0.05) and NUE (1.68, p < 0.05) (Fig. 2c). The emergence of this robust nonlinear pattern marks a key transition in advanced stages of development.

The mediating role of farm size also transforms fundamentally across income levels. Its association shifts from lower nitrogen yield in low-income countries to a central role in reducing machinery intensity, irrigation intensity, and nitrogen fertilizer application in high-income economies. The link between urbanization and farm size peaks during the upper-middle-income stage, when a one percentage point increase in urbanization corresponds to a 0.11% decrease in farm size, likely reflecting intensified land competition and peri-urban fragmentation. The association of farm size with lower input intensity per hectare becomes increasingly pronounced across development stages. In low-income countries, a 1% increase in farm size is associated with a significant 0.06% decrease in machinery intensity, while its associations with other inputs are not statistically significant (Fig. 2a). In the upper-middle-income stage, the same increase in farm size corresponds to significant declines of 0.05% in machinery intensity, 0.07% in irrigation intensity, and a marginally significant 0.09% reduction in nitrogen fertilizer application (Fig. 2b). By the high-income stage, this pattern becomes comprehensive: a 1% increase in farm size corresponds to decreases in machinery intensity, irrigation intensity, and nitrogen fertilizer application by 0.06%, 0.03%, and 0.10%, respectively (Fig. 2c). While farm size is weakly linked to lower nitrogen yield in low-income countries (−0.09%), it becomes positively associated with NUE in upper-middle-income economies (+0.05%).

The pathways governing nitrogen yield and NUE also reconfigure with development. In low- and lower-middle-income countries, nitrogen yield responds to multiple inputs: a 1% increase in machinery intensity, irrigation intensity, and nitrogen fertilizer application is associated with 0.02%, 0.03%, and 0.02% increases in nitrogen yield, respectively. In contrast, NUE is primarily constrained by nitrogen fertilizer application, with a 1% increase in nitrogen fertilizer application corresponding to a 0.08% decrease in NUE (Fig. 2a). As economies advance to upper-middle-income levels, the contribution of machinery intensity becomes negligible, while irrigation intensity emerges as a dominant factor, with a 1% increase in irrigation intensity yielding a 0.11% increase in nitrogen yield and a marginally significant 0.05% increase in NUE, indicating a potential synergy between productivity and efficiency (Fig. 2b). In high-income countries, the associations of both machinery intensity and irrigation intensity become insignificant, leaving nitrogen yield almost entirely dependent on f nitrogen fertilizer application, where a 1% increase in nitrogen fertilizer application corresponds to a 0.14% increase in nitrogen yield. Yet, this comes at a cost: the same increase in nitrogen fertilizer application is linked to a 0.28% decline in NUE, underscoring the growing trade-off between productivity and efficiency that characterizes advanced agricultural systems (Fig. 2c).

While the path analysis captures the average transmission channels between urbanization and nitrogen outcomes, it simplifies the inherently nonlinear and conditional nature of these relationships. To systematically examine these dynamics, we constructed conditional nonlinear marginal association curves (Fig. 3). These curves reveal two interrelated processes: the evolution of urbanization’s marginal association across its full spectrum and the systematic reshaping of this trajectory as a function of agricultural modernization. To illustrate this moderation, we compared "Low Modernization" and "High Modernization" scenarios, defined by holding farm size, machinery, and irrigation intensity at specific lower and upper bounds within each income group (see Methods for detailed visualization construction and robustness checks).

In low- and lower-middle-income countries, higher agricultural modernization corresponds to a more favorable trajectory in the association between urbanization and nitrogen outcomes. For nitrogen fertilizer application, higher modernization tempers the initial surge in nitrogen fertilizer application linked to urbanization and shifts the vertex of the inverted U-shaped curve from an urbanization level of 37% to 53% (Fig. 3a, d). For nitrogen yield, urbanization exerts a prolonged suppressive effect under low modernization. This aligns with recent evidence of widespread yield stagnation among African smallholders33, 34. This stagnation is driven by rural outmigration and capital constraints that hinder mechanization35. Our framework reveals that larger farm size acts as a critical structural buffer against these early losses. Scaling-up of farm size helps overcoming capital barriers, as productivity rises significantly only after farms reach a viable size for mechanization36. Consistent with this technological threshold, the impact of urbanization on nitrogen yield reverses dramatically under high modernization: the turning point of the U-shaped curve moves from 52% to 34% urbanization, indicating that the shift to a positive association occurs much earlier (Fig. 3b, e). Similarly, for NUE, higher modernization mitigates the initial efficiency losses associated with early urbanization. Although the turning point shifts from 30% to 50%, the entire curve is elevated upward, implying that the initial decline is less severe (Fig. 3c, f).

In upper-middle-income countries, higher modernization brings clear environmental co-benefits but also introduces nuanced trade-offs with efficiency gains. Under high modernization, the vertex of the U-shaped curve for nitrogen fertilizer application shifts markedly rightward, from 44% to 72% urbanization, prolonging the phase during which urbanization is associated with declining nitrogen fertilizer application (Fig. 3a, d). For nitrogen yield, the curvature of the relationship remains similar, but the entire curve is higher in the early stages and lower in the later stages of urbanization, crossing at the global mean (Fig. 3b, e). This pattern suggests a shift from maximizing output toward optimizing resource use. Under high modernization, the positive association between urbanization and NUE weakens notably: the vertex of the U-shaped curve shifts from beyond the observed range to an early 27% urbanization level, and the later-stage improvement in NUE becomes visibly smaller compared with the low modernization scenario (Fig. 3c, f). These results indicate that once agricultural modernization reaches a high level, further efficiency gains rely less on structural urban transitions and more on internal technological and managerial improvements.

In high-income countries, the relationship reveals a structural bottleneck, where modernization primarily mitigates the adverse effects of advanced urbanization. For nitrogen fertilizer application, the association follows a U-shaped curve, but higher modernization shifts the vertex rightward (from 68% to 77%), delaying the onset of recovery (Fig. 3a, d). For both nitrogen yield and NUE, urbanization is initially associated with positive outcomes that turn negative at very high urbanization levels. High modernization exerts a dual moderating effect: it dampens early positive associations while slowing the rate of decline at later stages. Consequently, the curve lies lower in early stages but higher in later stages, reflecting a substantial mitigation of eventual negative outcomes (Fig. 3b, c, e, f). Around the vertices, the lowest points of the curves are consistently higher under high modernization, signaling reduced vulnerability to adverse effects. This resilience is especially pronounced for NUE in highly urbanized contexts. At 75% urbanization, high modernization elevates the marginal effect on NUE from a negative value to near zero, while at 90% urbanization, it shifts the effect from –4% to +12%, representing a swing of more than 16 percentage points. This striking reversal underscores the pivotal capacity of agricultural modernization to buffer, and even counteract, the negative impacts of advanced urbanization. Identifying the specific modernization components responsible for this resilience is therefore a central focus of the subsequent analysis.

To identify which components of agricultural modernization most strongly reshape urbanization’s association with nitrogen outcomes, we decomposed the interaction terms into a moderation heatmap (Fig. 4). By evaluating urbanization’s marginal effects at early (25%), middle (50%), and late (75%) transition stages across low and high levels of each moderator, this framework reveals which specific technological or structural factor exerts the most dominant influence in a given developmental context (see Methods for analytical decomposition and robustness checks).

In low- and lower-middle-income countries, farm size emerges as the dominant moderator in the early phase of urbanization, mitigating adverse effects on nitrogen yield and NUE. In later stages, agricultural modernization contributes to buffering the negative marginal associations for nitrogen fertilizer application. Early in the urban transition, larger farm size plays a crucial stabilizing role, reversing the initially negative association of urbanization with nitrogen yield and NUE, thereby acting as a structural safeguard against early disruptions (Fig. 4). However, this influence is largely stabilizing rather than growth-inducing, as the effects hover near zero (Fig. S2). Machinery intensity and irrigation intensity exert limited moderating effects during this early phase but gain importance later, when all three factors collectively buffer the decline in nitrogen fertilizer application (Fig. 4). This coordinated response highlights a key developmental insight: while early urbanization challenges can be mitigated by scaling up farm size, sustained improvement requires multifaceted modernization. Such a combined buffering effect is particularly relevant where nitrogen fertilizer application remains below agronomic recommendations, making nitrogen yield stability more critical than further input reduction.

In upper-middle-income countries, the early stages of urbanization exhibit a more complex interplay of moderating effects. Machinery intensity and irrigation intensity intensities enhance the nitrogen fertilizer application-reducing effect of urbanization, producing clear environmental benefits. This reflects the input substitution phase where urban labor pull drives severe fertilizer overuse to compensate for lacking fixed inputs37, 38. However, our results demonstrate that this trade-off can be optimized. Consistent with recent empirical studies, land consolidation must be coupled with mechanization to achieve genuine efficiency gains and positive regional spillovers39, 40. This synergy provides the necessary scale to significantly reduce fertilizer intensity41. However, they also intensify the negative association of urbanization with nitrogen yield. Irrigation intensity plays a particularly dual role: it accelerates nitrogen fertilizer application reduction, offering strong environmental co-benefits, while simultaneously amplifying nitrogen yield losses (Fig. 4). In contrast, machinery intensity consistently acts as a stabilizer, mitigating early nitrogen yield declines and reinforcing later gains. Farm size expansion, however, exerts a persistently adverse moderating influence, intensifying downward pressures on both nitrogen yield and NUE (Fig. 4; Fig. S3). This divergence underscores a key policy challenge for middle-income economies: distinct modernization strategies carry distinct trade-offs. Policies prioritizing irrigation for environmental benefits may sacrifice nitrogen yield, while emphasizing consolidation may undermine efficiency. Effective navigation of this transitional stage requires balanced modernization portfolios, such as coupling irrigation investments with mechanization to sustain nitrogen productivity while maintaining environmental gains.

In high-income countries, modernization dynamics fundamentally shift. Farm size expansion, once a driver of efficiency, now acts as a negative moderator, exacerbating the decline in both nitrogen yield and NUE (Fig. 4). This reversal indicates diminishing or even negative returns to scale. It aligns with recent evidence that highly optimized and large-scale operations face yield stagnation and limited environmental benefits42-44. In contrast to further scale expansion, technological modernization through machinery intensity and irrigation intensity emerges as the principal source of resilience. These components consistently buffer the adverse impacts of high urbanization levels, albeit with a characteristic trade-off: slightly lower peak performance in early stages in exchange for greater stability in later phases (Fig. S4). In the advanced stages, machinery intensity becomes particularly decisive, being associated with both strengthened reductions in nitrogen fertilizer application and, crucially, a shift in the marginal association between urbanization and NUE from neutral to positive (Fig. 4). This nuanced finding refines the average positive role of farm size by delineating its limits. For high-income economies, sustainable nitrogen management in the face of continued urbanization will depend less on further farm consolidation and more on technology-driven intensification, emphasizing mechanization and irrigation efficiency to balance productivity, resource use, and environmental resilience.

Methods
We compiled national-level data from 139 countries spanning 1961 to 2021, sourced from FAOSTAT72, the USDA database73, FAO's Aqua-stat database74, World Bank75, and the CHANS model31, widely recognized datasets for agricultural analysis, to investigate how urbanization is associated with nitrogen use and nitrogen losses in global croplands. Using nitrogen as an indicator of agricultural intensification and environmental impact, we incorporated farm size, machinery intensity, and irrigation intensity as mediating variables and developed empirical regression models stratified by mediator levels to estimate how nitrogen fertilizer application, nitrogen yield and NUE respond to urbanization. Subsequently, we constructed an impact pathway, modeled using SEM, to trace how urbanization relates to nitrogen use indicators through its association with the mediators. Finally, we conducted a counterfactual analysis based on historical data and aggregated the results by income groups to assess how the estimated effects of urbanization on nitrogen dynamics vary across different economic contexts.

Data sources and processing. In this study, urbanization level, population, and cropland area data for the period 1961–2021 were obtained from FAOSTAT. PGDP (GDP per capita) data is obtained from World Development Indicators Database of World Bank75. Farm size was estimated as the cropland area per agricultural laborer. This metric provides a more direct measure of agricultural operational scale by linking land resources to the active agricultural workforce, thereby avoiding potential endogeneity that could arise from using total rural population as a denominator. Data for the agricultural labor force were obtained from the FAOSTAT database.

We conducted a regression analysis on the samples where ‘Cropland area per agricultural laborer’ and ‘Actual farm size’ could be matched by country and year, using the same fixed-effects settings as the main model (Table S1). Data for the latter, obtainable from the World Census of Agriculture76, is scarce, unevenly distributed, and lacks consistent statistical standards across countries, making it unsuitable for the main model's formal analysis. The results show that the coefficient for "Cropland area per agricultural laborer" is positive and highly significant. This indicates that after controlling for all country-level fixed characteristics and common time trends, the changes in "Cropland area per agricultural laborer" and "Actual farm size" are strongly and positively correlated. Furthermore, the goodness-of-fit for the model was 96.7%, suggesting that our proxy, in combination with country and year fixed effects, can account for 96.7% of the fluctuations in "Actual farm size". This demonstrates that the dynamic changes in "Cropland area per agricultural laborer" can indeed effectively capture the important trends in the changes in agricultural operation scale. For clarity and impact in discussing our results and policy implications, we refer to this proxy as 'farm size' throughout the text, with its precise definition and validation detailed herein.

Information regarding agricultural machinery was extracted from the September 2023 release of the USDA’s International Agricultural Productivity database73. This dataset supplies national-level estimates of total machinery stock, expressed in metric horsepower (CV), spanning the years 1961 through 2021. To capture the concentration of mechanized capital investment independent of population shifts, we calculated machinery intensity by dividing the total machinery power by the corresponding cropland area (CV ha⁻¹).

Irrigation statistics were obtained from the FAO AQUASTAT database. We explicitly utilized the metric 'actually irrigated area under full control irrigation'74. This specific category purposefully omits spate irrigation, equipped lowlands, and areas possessing inactive irrigation infrastructure. We omitted these elements because their reliance on floodwater means they remain highly vulnerable to climate fluctuations, unlike full-control systems. Moreover, their inclusion would artificially lower the calculated cropping intensity, distorting the actual agronomic advantages of active irrigation77,78. We then determined irrigation intensity by calculating the ratio of this actively controlled irrigated area to the total cropland area (%), serving as a proxy for the prevalence of advanced water management practices.

For the nitrogen budget, we relied on the CHANS model, which employs a comprehensive mass balance approach31. This framework combines top-down regional and global reactive nitrogen data with bottom-up flux estimations across 14 distinct subsystems (such as croplands, livestock, forests, wastewater, and atmospheric deposition). By synthesizing these diverse inputs and outputs, the CHANS model offers a holistic, multi-scalar evaluation of nitrogen flows and cycling.

In this study, the term 'nitrogen fertilizer application' specifically refers to the amount of nitrogen (N) applied to cropland from synthetic fertilizers, measured in kg N ha⁻¹. This is denoted as 𝑁 in our nitrogen budget calculations and is a core input flux derived from the CHANS model, which synthesizes data from sources including the International Fertilizer Association (IFA) and FAOSTAT. The total nitrogen output from cropland (𝑁 ) includes harvested nitrogen (𝑁 ), which is calculated as the sum of nitrogen content in the harvested grain and straw components. This calculation is based on annual crop production data from FAOSTAT and crop-specific nitrogen content coefficients established in the scientific literature. Our 'nitrogen yield' indicator, termed 'nitrogen harvested per hectare,' is calculated as the total modeled nitrogen harvested in crop grain and straw (𝑁 ) divided by the cropland area. This output flux is calculated within the CHANS model based on annual crop production data (from FAOSTAT) and crop-specific nitrogen content coefficients.

The nitrogen budget for cropland used in this study encompasses all crops (including staple crops and cash crops) from 1961 to 2021 and is calculated using the following equations (1)-(4):

𝑁 = 𝑁 +𝑁 +𝑁 +𝑁 +𝑁
(1)

𝑁 = 𝑁 +𝑁 +𝑁 +𝑁
(2)

𝑁 = 𝑁 +𝑁
(3)

𝑁𝑈𝐸 = 𝑁 / 𝑁
(4)

where 𝑁 represents the total nitrogen input to cropland, encompassing contributions from all crops. This includes 𝑁 , which accounts for synthetic nitrogen fertilizer applications; 𝑁 ,denoting organic nitrogen fertilizer use; 𝑁 , referring to biological nitrogen fixation; 𝑁 , representing nitrogen supplied through irrigation; and 𝑁 , signifying atmospheric nitrogen deposition.

𝑁 denotes the total nitrogen output from cropland, where 𝑁 refers to the total harvested nitrogen, encompassing both grain nitrogen (𝑁 ) and straw nitrogen (𝑁 ). Additionally, 𝑁 indicates nitrogen gas emissions from cropland, which include NH₃, NOx, N₂O, and N₂. 𝑁 reflects nitrogen runoff emissions from cropland, while 𝑁 pertains to the leaching of nitrogen into groundwater.

Statistical analysis

The urbanization level and PGDP data were standardized to a mean of 0 prior to statistical analysis to ensure comparability of regression coefficients. Additionally, due to the inherent right-skewed distribution of the raw PGDP data, a natural logarithm transformation was applied before standardization to better approximate a normal distribution. For the visualization of response curves, while the model inputs remained standardized values, the horizontal axis (urbanization level) was plotted using the original data scale to facilitate the practical interpretation of curve inflection points. This dual approach preserved computational accuracy in model estimation while retaining intuitive interpretability of results within the original measurement context. All results of main model are shown in Tables S2-S5.

Dynamic panel interaction model

To investigate the complex and potentially heterogeneous relationships between urbanization, agricultural modernization, and nitrogen outcomes, we developed and applied a dynamic panel interaction model. This analytical framework is designed to capture two critical features of the underlying data generating process. First, the 'dynamic' component acknowledges that nitrogen outcomes (e.g., nitrogen yield, NUE) exhibit strong temporal persistence due to systemic inertia, such as soil nitrogen legacy and the path-dependency of farming behaviors. We therefore incorporate a one-year lagged value of the dependent variable to account for this state dependence. Crucially, the inclusion of this autoregressive term directly dictates the interpretational role of our estimates. Because the lagged dependent variable absorbs the historical accumulation of effects, the estimated coefficients for urbanization and modernization must be interpreted strictly as short-run (contemporaneous) marginal effects. They quantify the immediate shift in nitrogen outcomes conditional on the system's state in the previous year, rather than the cumulative long-run equilibrium effects, thereby yielding a more conservative and precise measure of contemporary dynamics. Second, the 'interaction' component allows us to model how the effects of urbanization on nitrogen use indicators vary continuously and nonlinearly with the levels of key agricultural modernization mediators: farm size, machinery intensity, and irrigation intensity.

The urbanization level and all continuous control variables were centered (i.e., demeaned) prior to analysis. This ensures that the estimated coefficients are comparable and mitigates potential multicollinearity when including quadratic and interaction terms. As mentioned, PGDP was log-transformed prior to centering to address its right-skewed distribution.

Our core empirical specification is a multi-dimensional fixed-effects model, estimated using the following general form:

𝑌 = 𝑏 +𝑐 𝑌 + 𝑐 𝑈 +𝑐 𝑈2 +∑(𝑑 𝑀 + 𝑒 𝑀2 )+∑(𝑓 𝑈 ∙𝑀 )+𝑔′𝑊 +𝜂 +𝜏 +𝜁 +𝜈
(5)

where the subscripts i and t denote country and year, respectively. 𝑌 represents the dependent variables: the natural logarithms of nitrogen fertilizer application, nitrogen yield, and NUE. 𝑌 is the one-year lag of the dependent variable, with its coefficient 𝐜 capturing the degree of temporal persistence or autoregression. 𝑈 is the centered urbanization rate. Its quadratic term, 𝑈2, is included to capture the nonlinear (U-shaped or inverted U-shaped) effects of urbanization. 𝑀 is the k-th centered mediator of agricultural modernization, including farm size (cropland area per agricultural laborer), machinery intensity (CV ha⁻¹), and irrigation intensity (proportion of irrigated cropland). Their quadratic terms, 𝑀2 , are also included to account for their own nonlinear impacts. 𝑈 ∙𝑀 represents the interaction term between urbanization and each modernization mediator. The coefficients of primary interest, 𝑓 , capture how the marginal effect of urbanization on 𝑌 is moderated by the level of agricultural modernization. W is a vector of time-varying control variables, including the natural logarithm of per capita GDP (lnPGDP), precipitation, and temperature, along with their quadratic terms. To further account for system dynamics, this vector also includes cross-dynamic controls (e.g., lagged nitrogen fertilizer application, N , when modeling nitrogen yield, N ) to capture inter-temporal feedback between different nitrogen cycle components. 𝜂 , 𝜏 , and 𝜁 are country, year, and crop-type fixed effects, respectively. These absorb unobserved time-invariant heterogeneity at the country level (e.g., geography, culture), common time shocks affecting all countries (e.g., global technological trends), and variations specific to dominant crop systems. 𝜈 is the idiosyncratic error term, clustered at the country level to produce robust standard errors.

To systematically examine heterogeneity across development stages, we partitioned the full sample into three sub-groups based on World Bank income classifications: (1) Low- and Lower-middle-income, (2) Upper-middle-income, and (3) High-income countries. We estimated the complete interaction-term model independently for each sub-group. This classification is based on the World Bank's income group designations for the most recent year of our dataset (2021). This approach ensures that each country remains in the same group throughout the analysis, allowing us to interpret the results as the average developmental trajectory for countries that have achieved a certain income status. We acknowledge that some countries transitioned between groups during the 60-year study period, and we explore the implications for these transitioning economies in our Policy Implications section. This strategy allows all model coefficients, including those for quadratic and interaction terms, to vary freely across income groups, thereby avoiding the imposition of homogeneity constraints. This approach enables a granular examination of how the mechanisms linking urbanization to agricultural nitrogen cycles differ fundamentally across countries at varying levels of economic development.

The results from these models form the basis for our main findings, which are presented through a series of visualizations designed to illustrate the conditional marginal effects of urbanization across different scenarios of agricultural modernization and income levels.

To operationalize this for Fig. 3, the moderating variables (farm size, machinery, and irrigation intensity) were held constant at either the 10th percentile (defining the "Low Modernization" scenario) or the 90th percentile (defining the "High Modernization" scenario) within each income group. While this discrete comparison simplifies the visualization, it effectively represents the main trend of what is, in reality, a continuous moderation process. We verified this by plotting the smooth shifts and deformations of the effect curves across the full 1st-to-99th percentile spectrum of modernization levels (Figs. S5–S7), confirming that the 10th and 90th percentiles serve as robust representative bounds.

Furthermore, we rigorously tested the stability of these nonlinear patterns against historical outliers. A secondary analysis excluding countries that experienced major structural breaks (e.g., severe armed conflicts or state dissolution) yielded highly consistent curve shapes and vertices (Supplementary Note S3; Figs. S8-S10; Tables S6-S10). This confirms that the estimated nonlinear trajectories are driven by fundamental economic and agricultural transition mechanisms rather than anomalous geopolitical shocks. To address potential omitted variable bias from time-varying environmental factors, we conducted a second robustness check to ensure our findings are not sensitive to interannual climate variability (Fig. S11-13, Table S11-14).

To construct the moderation heatmaps (Fig. 4), we systematically decomposed the multidimensional interaction terms (𝑈 ∙𝑀 ) from the main empirical models. This involved evaluating the analytical derivative of the outcome with respect to urbanization at three representative urban transition stages (25%, 50%, and 75%), conditioned on the 10th and 90th percentiles of each specific moderator. The relative moderating strength visualized by the color intensity of each cell, was calculated as the difference in the marginal effect of urbanization between the high and low modernization scenarios, scaled by the baseline effect. The robustness of these distinct moderating patterns is further corroborated by the continuous nature of the underlying interactions, as mapped across the full data distribution (Figs. S5–S7).

Dynamic path analysis

To analyze the causal pathways linking urbanization to nitrogen outcomes while rigorously addressing the temporal dynamics inherent in our panel data, we replaced the conventional structural equation modeling (SEM) framework with a path analysis constructed from a system of dynamic panel regression models. This approach involves estimating a separate regression equation for each endogenous variable in our conceptual model (e.g., farm size, machinery intensity, nitrogen fertilizer application, nitrogen yield, and NUE). The estimated coefficients from this system of equations are then used to construct the final path diagram (Fig. 2).

Prior to estimation, all continuous variables were centered by subtracting their respective sample means to aid in the interpretation of coefficients. To handle extreme outliers, key agricultural modernization indicators (farm size, machinery intensity, and irrigation intensity) were winsorized at the 99th percentile.

The general form for each equation in the system is a two-way fixed-effects model with a lagged dependent variable:

𝑌 = 𝜌𝑌 +𝛽 𝑈 +𝛽 𝑈2 +𝛿′𝑍 +𝜇 +𝜀
(6)

where the subscripts i, t, and j denote country, year, and dominant crop type, respectively. 𝑌 is the centered outcome variable for a specific equation in the path model. 𝑌 is its one-year lag, included to capture the dynamic adjustment process and mitigate bias from serial autocorrelation. Empirically, our estimations yielded consistently high self-lagged coefficients (typically ranging from 0.6 to 0.9; Tables S2–S3). These high values reflect strong systemic inertia, such as the persistence of local farming habits and soil nutri ent carry-over. By explicitly absorbing these historical accumulation effects into the lagged term, our model successfully isolates the immediate net impact of urbanization. 𝑈 is the centered urbanization rate. Its quadratic term, 𝑈2 , is included to capture the nonlinear effects of urbanization on each endogenous variable. 𝑍 is a vector of control variables specific to each equation. It includes other endogenous variables that are causally prior in the path diagram (e.g., farm size is included when modeling machinery intensity), as well as a set of exogenous controls: centered ln(PGDP), precipitation, temperature, and their quadratic terms. To further account for system dynamics, this vector also includes lags of other relevant endogenous variables (e.g., lagged nitrogen yield and nitrogen fertilizer application when modeling NUE). These cross-lagged terms explicitly account for complex feedback loops where past outcomes influence current agricultural management decisions. 𝜇 , 𝜏 , and 𝜅 are country, year, and crop-type fixed effects, respectively. These absorb unobserved heterogeneity at the country level (e.g., geography), common time shocks (e.g., global technological trends), and variations specific to crop systems. This two-way fixed-effects specification also substantially mitigates concerns about non-stationarity by de-trending the series. 𝜀 is the idiosyncratic error term. Standard errors for all models are clustered at the country level to ensure robust statistical inferences in the presence of arbitrary within-country serial correlation and heteroskedasticity.

This entire system of equations was estimated separately for three distinct income sub-groups (Low- & Lower-middle, Upper-middle, and High-income), allowing all path coefficients to vary across different stages of economic development.

All coefficients reported in the path diagram (Fig. 2) represent the average marginal effect, which is calculated for a typical case where predictor variables are at their mean values. Because the dynamic specification controls for historical states, these coefficients must be strictly interpreted as short-run marginal effects, capturing the contemporaneous changes within a single year. This value quantifies the expected change in the outcome variable for a one-unit change in the predictor.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "I found two problems related to this paper. Inconsistent terminology and insufficient justification for proxy variables. Methodological limitations arise as conventional linear SEM fails to capture nonlinear urbanization effects on nitrogen use and inadequately addresses temporal dependencies in panel data, potentially biasing estimates. The details comments can be described as follow.\nIntroduction\n1. There are some inconsistencies in key terminology throughout the content. For instance, ‘nitrogen application’ appears in the Introduction section, while ‘fertilizer application' is used in the remaining content, including the Methods section. As commonly recognized, chemical fertilizers contain various nutrients such as nitrogen (N), phosphorus (P₂O₅), and potassium (K₂O). Please ensure logical consistency in keyword usage.\n2. Regarding the statement “Global crop yields, measured as nitrogen harvested per hectare”, I question how the authors quantified nitrogen harvested per hectare, as it is methodologically challenging to separate nitrogen harvest from overall crop yield measurements.\nMethodology\n1. The study uses “per capita cropland area” as a proxy for farm size due to data limitations. However, this approach raises concerns: 1) Why continue calling it ‘farm size’ when using a different metric? This could mislead readers. 2) The proxy’s validity depends on unverified assumptions about land distribution. While validation is mentioned, key details (like correlation strength) are buried in supplements. For clarity, the paper should either: a) consistently use ‘per capita cropland area’ throughout, or b) better justify why this proxy equals actual farm size in the main text.\n2. This paper use fertilizer application, yield and NUE to characterize nitrogen use and loss in global croplands, however, the quantification of fertilizer application and yield are omitted in the data sources and processing parts. The relationship between nitrogen and fertilizer application and yield are unclear in this paper.\n3. The paper employs structural equation modeling (SEM) to analyze longitudinal data spanning 1961–2021, I question the methodological appropriateness of applying conventional SEM to panel data without addressing critical temporal dynamics. SEM typically assumes cross-sectional covariance structures, yet panel data inherently involve autocorrelation, time-varying heterogeneity (i.e., unobserved factors that vary across both individuals and time periods, such as technological shocks or policy regimes), and potential non-stationarity—issues that, if unaddressed, may invalidate parameter estimates that was listed in Fig.4.\n4. The study’s core focus on nonlinear urbanization effects (e.g., phase-dependent NUE changes) conflicts with its use of linear SEM. Traditional SEM cannot capture the hypothesized inflection points or divergent regional responses, risking biased estimates. For example, SEM cannot capture threshold effects, saturation points, or the U-shaped relationships implied by the findings (e.g., low-income vs. high-income country divergence). It is urged testing nonlinear SEM specifications (e.g., interaction terms) or switching to inherently nonlinear methods.", "reviewer_2": "Review of the manuscript \"Nonlinear effect of urbanization on nitrogen use and loss in global croplands\"\nJournal: Nature Communications\nA. Manuscript Summary\nThis study investigates the nonlinear and context-dependent impacts of urbanization on nitrogen dynamics in global croplands, drawing on a panel dataset covering 139 countries from 1961 to 2021. Employing segmented regression and structural equation modeling, the authors disentangle the direct and mediated effects of urbanization on nitrogen fertilizer application, crop yield, and nitrogen use efficiency (NUE), with farm size, mechanization, and irrigation intensity as key intermediaries. The analysis reveals distinct phases in the urbanization-nitrogen nexus: in early urbanization, rising nitrogen inputs and declining NUE dominate, particularly in low-income settings with constrained technological infrastructure; whereas in more advanced urbanized economies, farm consolidation and technological adoption contribute to improved nitrogen management and yield gains. The study underscores the heterogeneous nature of urbanization’s agricultural impacts and calls for income-stratified, region-specific nitrogen management strategies to reconcile food security and environmental sustainability goals under ongoing global urban transitions.\nB. Comments\n1. Threshold Robustness and Generalizability Across Income Groups\nWhile the authors employ a data-driven segmented regression approach using AIC for threshold selection (e.g., 1 ha for farm size, 0.2 CV/cap for mechanization, 0.04 ha/cap for irrigation, and 55% urbanization rate), the decision to apply identical thresholds across all income strata, despite evident economic heterogeneity, remains methodologically debatable. The dataset spans 139 countries with wide disparities in structural agricultural characteristics, demographic transitions, and technological baselines.\nIn fact, the study’s own results, such as those presented in Fig. 5 and Extended Data Fig. 3, demonstrate that the impact of urbanization on fertilizer use, yield, and NUE differs markedly across income groups. For example, fertilizer application intensity increased by 37% in low-income countries but decreased by 15% in high-income countries, while yield responses range from −1% to +7% depending on the group. These differential trends suggest that the assumed threshold values may not function uniformly across contexts and may reflect pooled global averages more than context-specific inflection points.\nTo enhance robustness and ensure better policy relevance, we recommend the authors:\n• Re-estimate the thresholds separately within each income group (low, middle, high) using the same AIC-based segmented regression framework. This would directly test whether the thresholds are truly “universal” or vary by economic context.\n• Report and compare the group-specific thresholds and discuss any significant deviations from the pooled estimates. Divergence would indicate that applying uniform thresholds could obscure within-group dynamics.\nFinally, to strengthen the internal validity of the threshold identification, the authors could consider using alternative threshold detection techniques as a robustness check, such as Bayesian change point analysis, classification and regression trees (CART), or generalized additive models (GAMs) with breakpoint detection. These methods offer more flexible, non-parametric alternatives to segmented regression and may uncover threshold structures that AIC-based models miss or oversimplify.\nWithout these additional steps, the current threshold choices risk overgeneralizing the relationships across structurally diverse systems. Given the global scope and policy implications of the study, addressing this heterogeneity is critical for both scientific rigor and real-world relevance.\n2. Data Heterogeneity\nThe authors work with a large, heterogeneous panel of 139 countries. While fixed-effects models account for time-invariant heterogeneity, there is little discussion on how missing data, reporting inconsistencies, or structural breaks (e.g., due to war, institutional reforms) might bias results. This should be discussed more transparently, possibly supported by a robustness check (e.g., excluding countries with major conflicts, or sensitivity to reporting gaps).\n3. Environmental Heterogeneity and Omitted Variable Bias\nWhile the authors include a rich set of controls, including PGDP, crop type, agricultural investment, and fixed effects at the country and year level, the model may still omit time-varying environmental variables such as interannual climate variability (e.g., rainfall, temperature anomalies) and dynamic soil conditions, which are known to affect nitrogen use efficiency and yield response (e.g., Ray et al., 2015; Leng et al., 2016; Zhang et al., 2015). The use of fixed effects absorbs average differences across countries, but it does not address short-term agro-environmental shocks. A robustness check incorporating agro-ecological zone dummies or historical climate indices could help assess whether urbanization effects are overstated due to omitted variable bias.\n4. Y-Axis Scale Consistency for Comparative Interpretation\nIn Figures 5 and Extended Data Fig. 3, the Y-axis scales differ across income groups for the same indicators, which makes it difficult to visually compare the magnitude of urbanization-related effects across low-middle, upper-middle, and high-income countries.\nTo enhance interpretability and facilitate fair cross-group comparison, we recommend adopting consistent Y-axis scales across the three income groups for each indicator. This standardization would allow readers to more clearly assess the relative impact of urbanization across contexts and support one of the central aims of the study: to understand differential urbanization effects globally.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "This study offers a global-scale analysis of the complex interactions among urbanization, farm size, technological advancement, and irrigation infrastructure in shaping nitrogen use and agricultural productivity. The findings highlight the importance of region-specific management strategies that account for the nonlinear effects of urbanization, and provide targeted recommendations for nitrogen management across different stages of urban development. However, some aspects—particularly regarding summary of key findings, indicator selection, and the interpretation of results—would benefit from further analysis, clarification, and refinement. My comments are as follows:\n1. Although the data and modeling work support the points made in the second paragraph of the Introduction, the findings primarily reinforce these existing arguments rather than delivering substantially new perspectives. This reduces the perceived depth and novelty of the results. I suggest the authors more explicitly analyze and articulate the novel aspects of the specific findings derived from the calculation.\n2. Regarding indicator selection, the calculations of farm size, mechanization ownership, and irrigated area per capita all use rural population as the denominator. Since urbanization level is also determined by the proportion of rural population, these indicators may have endogenous relationships with urbanization itself. Does this imply that, all else being equal, increases in urbanization level are necessarily associated with increases in farm size, mechanization, and irrigation? How do the authors account for this potential endogeneity and its influence on the results?\n3. The first two sentences of the abstract involve multiple concepts. Specifically, the first sentence states that urbanization impacts land use and population distribution, while the second emphasizes a particular focus on mechanization and irrigation. The logical relationships among these keywords are not clear. I recommend refining these sentences to make them more focused and concise, and explicitly defining the dimensions of urbanization addressed in this study. The same issue exists in the first three sentences of the main text. In other words, in lines 62-63, when stating that urbanization “accounts for... socioeconomic dimensions,” it would be helpful to specify which dimensions are being considered.\n4. The results section of the abstract would be strengthened by including key quantitative findings, especially regarding changes in NUE or significant thresholds identified in the study.\n5. Lines 53-56: What specific characteristics of urbanization in low-income countries (as defined by the conceptual framework above) lead to potentially different impacts on nitrogen use compared to high-income countries? A discussion of this point would help readers better understand the rationale for conducting a global analysis and the importance of paying special attention to the transition from low- to high-income contexts.\n6. Lines 89-91 state that the decline in Eastern Europe potentially stems from its naturally high soil fertility and extensive black soils. I am unsure about the specific timeframe during which this decline occurred. Could this decline also be related to socioeconomic developments, such as economic transformation, agricultural decline, or environmental policies?\n7. I suggest that the authors include the analysis of typical countries within the results section to support and illustrate the main conclusions. For example, in lines 101-103, when referring to \"some countries,\" it would be helpful to specify which countries are involved and describe their experiences with unfavorable climatic conditions, soil degradation, or limited access to critical agricultural inputs. This would help readers gain a clearer understanding of the issues.\n8. Please check the numerical data presented in line 197.\n9. During the study period, some countries experienced significant income level changes (e.g., transitioning from low- to middle-income). How are these countries categorized in the temporal analysis? Given that the authors mention the transition from low- to high-income contexts in the introduction, could the authors use typical examples of countries that have undergone such transitions to further elucidate their findings?\n10. Please discuss the limitations of this study and suggest directions for future research.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to Referees’ Comments\nReferee #1\nComment Response and Change\nI found two problems related to this paper. Inconsistent terminology and R1-1: Thank you for your insightful and clear summary of the key issues in\ninsufficient justification for proxy variables. Methodological limitations our manuscript. We sincerely appreciate your valuable feedback. We have\narise as conventional linear SEM fails to capture nonlinear urbanization addressed your specific comments in the following sections.\neffects on nitrogen use and inadequately addresses temporal\ndependencies in panel data, potentially biasing estimates. The details\ncomments can be described as follow.\nIntroduction R1-2: Thank you for flagging the inconsistency. To improve clarity and consistency, we standardize\n1. There are some inconsistencies terminology to fertilizer N application, which refers to the N component of fertilizer inputs and serves as our\nin key terminology throughout indicator of agricultural intensification and associated environmental impacts. In the Methods, we define\nthe content. For instance, fertilizer N application as the mass of N applied to cropland via mineral and organic fertilizers, aligned with our\n‘nitrogen application’ appears in N budget framework. Throughout the manuscript and figures, instances of fertilizer application or nitrogen\nthe Introduction section, while application have been replaced with fertilizer N application.\n‘fertilizer application' is used in\nthe remaining content, including 1. Standardization of terminology: To ensure the highest degree of precision, we have standardized the\nthe Methods section. As term to \"nitrogen fertilizer application\" throughout the main text, including the Introduction, Results,\ncommonly recognized, chemical and Discussion sections. This change makes it unequivocally clear that our analysis is focused specifically\nfertilizers contain various on the nitrogen (N) input from synthetic fertilizers.\nnutrients such as nitrogen (N),\nphosphorus (P₂O₅), and 2. Exception in the abstract: The only exception is in the abstract. Due to strict word count limitations, after\npotassium (K₂O). Please ensure first establishing the focus on nitrogen, we occasionally use the shorter term \"fertilizer \" for conciseness. We\nlogical consistency in keyword are confident that the context provided within the abstract prevents any ambiguity.\nusage.\n3. Explicit definition in methods: To fundamentally anchor this terminology, we have also added a precise\ndefinition in the Methods section, under the 'Data sources and processing' subsection (Page 10, Line 482).\nThe sentence reads:\n\"In this study, the term 'nitrogen fertilizer application' specifically refers to the amount of nitrogen (N) applied\nto cropland from synthetic fertilizers, measured in kg N ha⁻¹. This is denoted as N in our nitrogen budget\nfertilizer\ncalculations.\"\n2. Regarding the statement R1-3: Thank you for noting the ambiguity. In our analysis, nitrogen harvest is estimated with a standard mass\n“Global crop yields, measured balance that integrates FAOSTAT crop production data and crop-specific nitrogen content coefficients from the\nas nitrogen harvested per literature. The revision now states this distinction explicitly.\nhectare”, I question how the\nauthors quantified nitrogen 1. Clarification in the main text:\nharvested per hectare, as it is We have revised the sentence in the Results and discussion section to immediately specify the modeled nature\nmethodologically challenging to of this metric.\nseparate nitrogen harvest from Revision in the manuscript (Page 3, Line 105):\noverall crop yield \"Global crop nitrogen yield, defined as nitrogen harvested per hectare based on model estimates...\"\nmeasurements.\n2. Detailed explanation in the Methods section:\nWe have added a comprehensive explanation of the calculation process in the Methods section under \"Data\nsources and processing\".\nNew text in the Methods section (Page 10, Lines 486-493):\n\"The total nitrogen output from cropland (N ) includes harvested nitrogen (N ), which is calculated as\noutput harvest\nthe sum of nitrogen content in the harvested grain and straw components. This calculation is based on annual\ncrop production data from FAOSTAT and crop-specific nitrogen content coefficients established in the\nscientific literature. Our 'nitrogen yield' indicator, termed 'nitrogen harvested per hectare,' is calculated as the\ntotal modeled nitrogen harvested in crop grain and straw (N ) divided by the cropland area. This output\nharvest\nflux is calculated within the CHANS model...\"\nMethodology R1-4: Thank you for these thoughtful and important comments regarding our proxy for 'farm size'. We have\n1. The study uses “per capita substantially revised the MS with incorporating your suggestions to enhance the paper's clarity and rigor. We\ncropland area” as a proxy for have: 1) upgraded to a more precise proxy, 2) moved the robust validation (showing a strong, significant\nfarm size due to data limitations. correlation and high explanatory power) into the main text, and 3) implemented a clear and consistent\nHowever, this approach raises communication strategy regarding the terminology.\nconcerns: 1) Why continue\ncalling it ‘farm size’ when using 1. Regarding the terminology ('Farm Size' and 'Cropland area per agricultural laborer'):\na different metric? This could To retain the term 'farm size' in the main narrative (e.g., Abstract, Discussion) is a deliberate choice aimed at\nmislead readers. 2) The proxy’s ensuring our findings are accessible and impactful for a broad audience, including policymakers. 'Farm size' is\nvalidity depends on unverified a cornerstone concept in agricultural policy, and using a more technical proxy name throughout would hinder\nassumptions about land the communication of our core, policy-relevant conclusions. To ensure complete academic transparency, we\ndistribution. While validation is have implemented a robust framing strategy throughout the manuscript.\nmentioned, key details (like\ncorrelation strength) are buried 1) Clear definition: In the Methods section, we now explicitly define our proxy. We have also upgraded our\nin supplements. For clarity, the proxy to the more precise 'cropland area per agricultural laborer' for greater theoretical soundness.\npaper should either: a) Revision in the Methods section (Page 9, Lines 430-432 & 449-451):\nconsistently use ‘per capita \"Farm size was estimated as the cropland area per agricultural laborer. This metric provides a more direct\ncropland area’ throughout, or b) measure of agricultural operational scale...\"\nbetter justify why this proxy \"...For clarity and impact in discussing our results and policy implications, we refer to this proxy as 'farm size'\nequals actual farm size in the throughout the text, with its precise definition and validation detailed herein.\"\nmain text.\n2) Consistent clarification: To continuously remind the reader, we have added a clarifying note to the caption\nof every figure and table that uses this variable.\nExample of revision in figure captions (e.g., Fig. 2, Line 683):\n\"...'Farm size' is proxied by cropland area per agricultural laborer.\"\nThis approach allows us to maintain a clear and impactful narrative while upholding rigorous standards of\ntransparency.\n2. Regarding the validity and justification of the proxy:\nWe agree that the evidence for our proxy's validity should be a prominent part of the main text, not \"buried in\nsupplements.\" We have moved the full validation analysis into the main Methods section. We now demonstrate\nnot that the proxy is 'equal' to actual farm size, but that it robustly captures the dynamic trends in actual farm\nsize, which is the requirement for a valid proxy in a panel data context.\nNew paragraph added to the Methods section (Page 9, Lines 437-449):\n\"We conducted a regression analysis on the samples where ‘Cropland area per agricultural laborer’ and 'Actual\nfarm size' could be matched by country and year, using the same fixed-effects settings as the main model (Table\nS1). Data for the latter, obtainable from the World Census of Agriculture, is scarce, unevenly distributed, and\nlacks consistent statistical standards across countries, making it unsuitable for the main model's formal analysis.\nThe results show that the coefficient for \"Cropland area per agricultural laborer\" is positive and highly\nsignificant. This indicates that after controlling for all country-level fixed characteristics and common time\ntrends, the changes in \"Cropland area per agricultural laborer\" and \"Actual farm size\" are strongly and\npositively correlated. Furthermore, the goodness-of-fit for the model was 96.7%, suggesting that our proxy,\nin combination with country and year fixed effects, can account for 96.7% of the fluctuations in \"Actual\nfarm size\". This demonstrates that the dynamic changes in \"Cropland area per agricultural laborer\" can indeed\neffectively capture the important trends in the changes in agricultural operation scale.\"\nThe key statistical evidence from this analysis is now clearly presented for the reviewer's convenience (and is\ndetailed in Supplementary Table S1):\nSummary of Validation Results (from Table S1):\nTable S1 Correlation between cropland area per agricultural laborer and actual farm size.\nEach column represents a separate regression model. All equations include country and year fixed effects. The\nasterisks indicate the statistical significance level based on P values: *p ≤ 0.1, **p ≤ 0.05, ***p ≤ 0.01,\nwhich assessed using a two-sided t-test.\nLn Actual farm size\nLn Cropland area per agricultural laborer 0.417***\ns.e. (0.051)\nP values 0.000\nCountry Yes\nYear Yes\nN 418\nAdjust-R2 0.967\nReferences:\nFAOSTAT. World Census of Agriculture Database; https://www.fao.org/faostat/en/#data/WCAD. (2025).\n2. This paper use fertilizer R1-5: Thank you for highlighting this lack of clarity. To rectify this, we have substantially revised the Data\napplication, yield and NUE to sources and processing subsection within the Methods to provide precise, nitrogen-based definitions for all\ncharacterize nitrogen use and key variables.\nloss in global croplands,\nhowever, the quantification of 1. Quantification of \"fertilizer application\" and \"yield\":\nfertilizer application and yield We have added new paragraphs to the Methods section to explicitly detail the quantification of these terms.\nare omitted in the data sources 1) For \"Fertilizer Application\": We now clarify that this term exclusively refers to synthetic nitrogen\nand processing parts. The fertilizer.\nrelationship between nitrogen Revision in the Methods section (Page 10, Lines 482-486):\nand fertilizer application and \"In this study, the term 'nitrogen fertilizer application' specifically refers to the amount of nitrogen (N)\nyield are unclear in this paper. applied to cropland from synthetic fertilizers, measured in kg N ha⁻¹. This is denoted as N in our\nfertilizer\nnitrogen budget calculations and is a core input flux derived from the CHANS model, which synthesizes\ndata from sources including the International Fertilizer Association (IFA) and FAOSTAT.\"\n2) For \"Yield\": We now clarify that our yield metric is a measure of nitrogen harvested, not simply crop\nbiomass.\nRevision in the Methods section (Page 10, Lines 490-493):\n\"Our 'nitrogen yield' indicator, termed 'nitrogen harvested per hectare,' is calculated as the total modeled\nnitrogen harvested in crop grain and straw (Nharvest) divided by the cropland area. This output flux is\ncalculated within the CHANS model based on annual crop production data (from FAOSTAT) and crop-\nspecific nitrogen content coefficients.\"\n2. Clarifying the relationship between nitrogen, fertilizer, and yield:\nWith these explicit definitions now in place, the relationship between these key terms within our nitrogen\nbudget framework is unambiguous. This relationship is mathematically defined by the nitrogen budget\nequations already present in our Methods section, which we respectfully highlight here for clarity.\nFrom the Methods section (Pages 10-11, Lines 499, 501, 505):\nN = N + N + ... (1)\ninput fertilizer manure\nN = N + N + ... (2)\noutput harvest gas\nNUE = N / N (4)\nharvest input\nAs these equations and our new definitions illustrate:\nWith these explicit definitions now in place, the relationship between these key terms within our nitrogen budget\nframework is unambiguous:\n\"Nitrogen fertilizer application\" is the primary synthetic nitrogen input (N ).\nfertilizer\n\"Nitrogen yield\" is the nitrogen output captured in the harvested crop (N ).\nharvest\nNitrogen Use Efficiency (NUE) is therefore consistently calculated as the ratio of this nitrogen output (N )\nharvest\nto the total nitrogen input (N ), as shown in our nitrogen budget equations.\ninput\n3. The paper employs structural R1-6: Thank you for these insightful and specific technical questions regarding the methodological\nequation modeling (SEM) to appropriateness of applying conventional SEM to our longitudinal dataset. To address these issues, we\nanalyze longitudinal data have replaced the conventional SEM framework with a more appropriate and robust approach: a path\nspanning 1961–2021, I question analysis constructed from a system of dynamic panel regression models. This allows us to retain the intuitive\nthe methodological path-diagram structure while rigorously modeling the temporal dynamics inherent in our data. The results of\nappropriateness of applying this new approach are presented in the revised Fig. 2, and the full methodology is detailed in the new \"Dynamic\nconventional SEM to panel data path analysis\" subsection of the Methods. Each path in our model is now estimated using a two-way fixed-\nwithout addressing critical effects (TWFE) model that includes a lagged dependent variable. We have detailed this specification in\ntemporal dynamics. SEM the Methods section (Lines 647-691).\ntypically assumes cross-\nsectional covariance structures, The general form for each equation is:\nyet panel data inherently involve\nautocorrelation, time-varying From the Methods section (Page 13, Line 611):\nheterogeneity (i.e., unobserved\nfactors that vary across both 𝑌 = 𝜌𝑌 +𝛽 𝑈 +𝛽 (𝑈 )2 +𝛿′𝑍 +𝜇 +𝛾 +𝜅 +𝜀\n𝑖𝑡 𝑖,𝑡−1 1 𝑖𝑡 2 𝑖𝑡 𝑖𝑡 𝑖 𝑡 𝑗 𝑖𝑡\nindividuals and time periods,\nsuch as technological shocks or This framework systematically resolves the econometric challenges you identified:\npolicy regimes), and potential\nnon-stationarity—issues that, if 1. Regarding Autocorrelation:\nunaddressed, may invalidate Our dynamic panel model addresses serial correlation in two robust ways:\nparameter estimates that was\nlisted in Fig.4. 1) Lagged dependent variable: We explicitly model persistence by including the lagged dependent variable\n(𝑌 ) in every equation.\n𝑖,𝑡−1\nAs stated in the Methods section (Page 13, Lines 614-616):\n\"𝑌 is its one-year lag, included to capture the dynamic adjustment process and mitigate bias from serial\n𝑖,𝑡−1\nautocorrelation.\"\n2) Clustered standard errors: To handle any remaining arbitrary serial correlation, we cluster standard errors\nat the country level.\nAs stated in the Methods section (Page 13, Lines 627-629):\n\"Standard errors for all models are clustered at the country level to ensure robust statistical inferences in the\npresence of arbitrary within-country serial correlation and heteroskedasticity.\"\n2. Regarding Time-Varying Heterogeneity:\nOur TWFE specification is the standard and most effective approach for controlling for unobserved time-varying\ncommon shocks (e.g., technological waves, global price shocks).\nAs explained in the Methods section (Page 13, Lines 624-625):\n\"These [𝜇 , 𝛾 , and 𝜅 ] absorb unobserved heterogeneity at the country level (e.g., geography), common time\n𝑖 𝑡 𝑗\nshocks (e.g., global technological trends), and variations specific to crop systems.\"\n3. Regarding Non-Stationarity:\nWe acknowledge the potential for non-stationarity in our 60-year time series. The TWFE estimator is well-suited\nto mitigate this issue. By including country fixed effects (𝜇 ), we effectively de-trend the data for country-\n𝑖\nspecific trends. By including year fixed effects (𝛾 ), we control for common stochastic trends affecting all\n𝑡\ncountries.\nAs noted in the Methods section (Page 13, Lines 626-627):\n\"This two-way fixed-effects specification also substantially mitigates concerns about non-stationarity by de-\ntrending the series.\"\n4. The study’s core focus on R1-7: Thank you for this critical methodological point. To address this, we have replaced the linear SEM with\nnonlinear urbanization effects a comprehensive, multi-stage approach using dynamic panel models that are inherently nonlinear. We first use\n(e.g., phase-dependent NUE a dynamic panel model with quadratic terms to establish the average nonlinear pathways (Fig. 2), and then\nchanges) conflicts with its use of employ a full dynamic panel interaction model to explore how these pathways are moderated, explaining the\nlinear SEM. Traditional SEM complex heterogeneity in our data (Figs. 3 & 4). Our new analytical framework addresses nonlinearity in two\ncannot capture the hypothesized distinct, sequential steps, using two different model specifications to answer different questions.\ninflection points or divergent\nregional responses, risking Step 1: Establishing the average nonlinear pathways (for Fig. 2)\nbiased estimates. For example, First, to establish the foundational transmission pathways, we developed a system of dynamic panel\nSEM cannot capture threshold models where each equation includes a quadratic term for urbanization. This allows us to capture\neffects, saturation points, or the the average U-shaped or inverted U-shaped relationships, correcting the misspecification of the previous linear\nU-shaped relationships implied model.\nby the findings (e.g., low-\nincome vs. high-income country As specified in the \"Dynamic path analysis\" subsection of the Methods (Lines 608-617):\ndivergence). It is urged testing \"The general form for each equation in the system is a two-way fixed-effects model...:\nnonlinear SEM specifications\n(e.g., interaction terms) or 𝑌 = 𝜌𝑌 +𝛽 𝑈 +𝛽 (𝑈 )2 +𝛿′𝑍 +𝜇 +𝛾 +𝜅 +𝜀\n𝑖𝑡 𝑖,𝑡−1 1 𝑖𝑡 2 𝑖𝑡 𝑖𝑡 𝑖 𝑡 𝑗 𝑖𝑡\nswitching to inherently\nnonlinear methods. ...Its quadratic term, (𝑈 )2, is included to capture the nonlinear effects of urbanization on each endogenous\n𝑖𝑡\nvariable.\"\nThe path diagram in Fig. 2 is constructed from this system of equations. For clarity, the arrows are annotated\nwith the linear term's coefficient (𝛽 ) to represent the average marginal effect across the data range, while the\n1\npresence of significant quadratic terms (𝛽 ), detailed in Supplementary Tables S2-S3, confirms the underlying\n2\nnonlinearity.\nStep 2: Deconstructing the conditional nonlinear dynamics (for Figs. 3 and 4)\nSecond, to answer the deeper questions about how these nonlinear relationships change and what drives the\n\"divergence\" you mentioned, we estimated a more complex dynamic panel interaction model. This model not\nonly includes the quadratic term for urbanization but also interacts urbanization with our key modernization\nvariables (farm size, machinery, irrigation).\nAs specified in the \"Dynamic panel interaction model\" subsection of the Methods (Lines 549-552):\n\"Our core empirical specification is a multi-dimensional fixed-effects model...:\n𝐾\n𝑌 = ⋯+ 𝑐 𝑈 +c 𝑈2 +⋯+∑(𝑓 𝑈 ∙𝑀 )+⋯\n𝑖𝑡 1 𝑖𝑡 2 𝑖𝑡 𝑘 𝑖𝑡 𝑘𝑖𝑡\n𝑘=1\nThis interaction model is the engine for our main nonlinear findings, which are visualized in Figures 3 and 4:\nFig. 3 (Conditional nonlinear effect curves): This figure is generated from the interaction model by calculating\nthe marginal effect of urbanization at different levels of modernization.\nAs described in the Results section (Lines 195-201):\n\"...we constructed conditional nonlinear marginal association curves (Fig. 3). These curves reveal two\ninterrelated processes: the evolution of urbanization’s marginal association across its full spectrum and the\nsystematic reshaping of this trajectory as a function of agricultural modernization. To illustrate this moderation,\nfarm size, machinery intensity, and irrigation intensity were held constant at either the 10th percentile (Low\nModernization scenario) or the 90th percentile (High Modernization scenario) ...\"\nFig. 4 (Heatmap of moderating effects): This figure further deconstructs the interaction terms (𝑓 ) from the\n𝑘\nsame model to quantify the moderating effects at specific stages of urbanization.\nAs explained in the Results section (Lines 253-255):\n\" To identify which components of agricultural modernization most strongly reshape urbanization’s association\nwith nitrogen outcomes, we decomposed the interaction terms from the main model. The resulting heatmap (Fig.\n4) ...\"\nReferee #2\nComment Response and Change\nA. Manuscript Summary R2-1: Thank you for\nThis study investigates the nonlinear and context-dependent impacts of urbanization on nitrogen dynamics in global providing such a clear and\ncroplands, drawing on a panel dataset covering 139 countries from 1961 to 2021. Employing segmented regression and accurate summary of our\nstructural equation modeling, the authors disentangle the direct and mediated effects of urbanization on nitrogen fertilizer manuscript's core\napplication, crop yield, and nitrogen use efficiency (NUE), with farm size, mechanization, and irrigation intensity as key objectives, methodology,\nintermediaries. The analysis reveals distinct phases in the urbanization-nitrogen nexus: in early urbanization, rising and findings. We have\nnitrogen inputs and declining NUE dominate, particularly in low-income settings with constrained technological addressed your specific\ninfrastructure; whereas in more advanced urbanized economies, farm consolidation and technological adoption contribute comments in the following\nto improved nitrogen management and yield gains. The study underscores the heterogeneous nature of urbanization’s sections.\nagricultural impacts and calls for income-stratified, region-specific nitrogen management strategies to reconcile food\nsecurity and environmental sustainability goals under ongoing global urban transitions.\nB. Comments R2-2: Thank you for this profound and highly constructive critique. In the revised MS, we have abandoned\n1. Threshold Robustness and the threshold-based segmented regression approach. Our new framework directly models the continuous\nGeneralizability Across Income nature of these relationships using a dynamic panel interaction model, estimated separately for each income\nGroups group. This new strategy, which forms the core of our revised manuscript, is detailed in the Methods\nWhile the authors employ a data-driven section and the Results section.\nsegmented regression approach using\nAIC for threshold selection (e.g., 1 ha 1. Methodological upgrade: From discrete thresholds to continuous interactions\nfor farm size, 0.2 CV/cap for Our new core analytical tool is a dynamic panel model that includes quadratic terms and,\nmechanization, 0.04 ha/cap for crucially, interaction terms between urbanization and each modernization variable. This allows us to ask a\nirrigation, and 55% urbanization rate), more nuanced question: \"How does the nonlinear effect of urbanization continuously change as the level of\nthe decision to apply identical agricultural modernization varies?\"\nthresholds across all income strata,\ndespite evident economic The core of our new specification is detailed in the Methods section (Lines 549-564):\nheterogeneity, remains 𝐾\nmethodologically debatable. The 𝑌 = ⋯+ 𝑐 𝑈 +c 𝑈2 +⋯+∑(𝑓 𝑈 ∙𝑀 )+⋯\n𝑖𝑡 1 𝑖𝑡 2 𝑖𝑡 𝑘 𝑖𝑡 𝑘𝑖𝑡\n𝑘=1\ndataset spans 139 countries with wide \"...The coefficients of primary interest, 𝑓 , capture how the marginal effect of urbanization on 𝑌 is\n𝑘 𝑖𝑡\ndisparities in structural agricultural moderated by the level of agricultural modernization...\"\ncharacteristics, demographic\ntransitions, and technological Most importantly, to fully address your concern about heterogeneity, we do not estimate one pooled model.\nbaselines. Instead, we estimate this full interaction model independently for each of the three income groups.\nIn fact, the study’s own results, such as As explained in the Methods section (Lines 575-586):\nthose presented in Fig. 5 and Extended \"To systematically examine heterogeneity across development stages, we partitioned the full sample into\nData Fig. 3, demonstrate that the three sub-groups based on World Bank income classifications... We estimated the complete interaction-term\nimpact of urbanization on fertilizer use, model independently for each sub-group... This strategy allows all model coefficients, including those for\nyield, and NUE differs markedly across quadratic and interaction terms, to vary freely across income groups, thereby avoiding the imposition of\nincome groups. For example, fertilizer homogeneity constraints.\"\napplication intensity increased by 37%\nin low-income countries but decreased This approach is superior because it allows the entire functional form of the relationship to differ across\nby 15% in high-income countries, income groups, a much more flexible approach than simply allowing a single threshold point to vary.\nwhile yield responses range from −1%\nto +7% depending on the group. These 2. New visualizations to display heterogeneous and nonlinear effects\ndifferential trends suggest that the We have designed a new set of visualizations to present the rich findings from this model, directly addressing\nassumed threshold values may not the \"differential trends\" you noted.\nfunction uniformly across contexts and\nmay reflect pooled global averages Fig. 3 (\"Agricultural modernization reshapes the nonlinear effects of urbanization... \"): This figure now\nmore than context-specific inflection serves as the centerpiece of our results. It plots the conditional nonlinear effect curves, showing how the\npoints. shape, vertex, and position of the urbanization effect curve differ dramatically across income groups and\nmodernization levels. This directly visualizes the heterogeneity that a uniform threshold would have\nTo enhance robustness and ensure obscured.\nbetter policy relevance, we recommend\nthe authors: Fig. 4 (\"Agricultural modernization reshapes the effects of urbanization by exerting divergent\nRe-estimate the thresholds separately moderating roles...\"): This heatmap deconstructs the interaction terms from our model, quantitatively\nwithin each income group (low, showing which modernization factors are most important at different stages of urbanization for each income\nmiddle, high) using the same AIC- group. It provides the mechanistic explanation for the divergent pathways.\nbased segmented regression\nframework. This would directly test Supplementary Figs. S1-S3 (\"Continuous moderation spectrum...\"): To show the full continuous nature\nwhether the thresholds are truly of the moderation, these figures visualize a spectrum of 100 curves, illustrating how the effect of urbanization\n“universal” or vary by economic\ncontext. smoothly \"twists\" and \"shifts\" as modernization levels change. This provides the full, continuous evidence\nthat underpins the discrete scenarios shown in Fig. 3.\nReport and compare the group-specific\nthresholds and discuss any significant 3. Regarding alternative Methods (e.g., GAMs)\ndeviations from the pooled estimates. Finally, we wish to address your thoughtful suggestion regarding the use of alternative threshold detection\nDivergence would indicate that techniques like Bayesian change point analysis, CART, or GAMs. We fully agree that these non-parametric\napplying uniform thresholds could methods offer valuable flexibility for exploratory data analysis. We adopt a parametric interaction model,\nobscure within-group dynamics. rather than pursuing these alternatives, was a deliberate methodological choice based on three key\nconsiderations:\nFinally, to strengthen the internal\nvalidity of the threshold identification, 1) Our goal shifted from \"finding a breakpoint\" to \"modeling a continuous process\":\nthe authors could consider using Your initial critique rightly pointed out the weakness of a single, discrete threshold. This led us to realize that\nalternative threshold detection the core research question is not \"At what point does the relationship change?\" but rather \"How does the\ntechniques as a robustness check, such relationship change continuously as modernization progresses?\" Methods like change point analysis or\nas Bayesian change point analysis, CART are excellent for the first question (finding discrete breaks), but our new interaction model is\nclassification and regression trees specifically designed to answer the second, more nuanced question. It models the marginal effect of\n(CART), or generalized additive urbanization as a continuous function of the modernization variables, which aligns perfectly with our revised\nmodels (GAMs) with breakpoint theoretical framing. In essence, we have moved beyond the need for \"breakpoint detection\" altogether.\ndetection. These methods offer more\nflexible, non-parametric alternatives to 2) Theoretical interpretability and hypothesis testing:\nsegmented regression and may uncover Our study is grounded in economic and agronomic theory, which strongly suggests the existence of U-shaped\nthreshold structures that AIC-based or inverted U-shaped relationships (e.g., due to initial resource strain followed by economies of scale, or\nmodels miss or oversimplify. initial gains followed by diminishing returns). Our quadratic specification (β₁U + β₂U²) allows for a direct\nand formal test of these theoretically-derived hypotheses. While GAMs can flexibly fit such shapes, they do\nWithout these additional steps, the not provide easily interpretable parameters (like the sign and significance of β₂) to confirm or reject these\ncurrent threshold choices risk specific theoretical propositions. Our parametric approach provides a clearer link between theory, model, and\novergeneralizing the relationships results.\nacross structurally diverse systems.\nGiven the global scope and policy\n3) Quantifying moderating effects for policy implications:\nimplications of the study, addressing A primary goal of our paper is to provide actionable policy insights by quantifying how much each component\nthis heterogeneity is critical for both of modernization (farm size, machinery, irrigation) alters the urbanization pathway. Our interaction model\nscientific rigor and real-world\n(𝑓 𝑈 ∙𝑀 ) is perfectly suited for this. The interaction coefficients (𝑓 ) and the analyses presented in our\n𝑘 𝑖𝑡 𝑘𝑖𝑡 𝑘\nrelevance.\nheatmap (Figure 4) provide a direct, quantitative measure of these moderating effects. Achieving this level\nof quantitative decomposition and clear attribution with more flexible, non-parametric models like GAMs\nwould be significantly more complex and less direct.\nIn summary, while we appreciate the power of the methods you suggested, our fundamentally upgraded\nstrategy no longer requires \"threshold detection.\" Instead, we have chosen a parametric (but highly flexible\nand interactive) framework because it better aligns with our revised research goal of modeling a continuous\nprocess, allows for direct hypothesis testing of U-shaped relationships, and yields more directly interpretable\nquantitative insights for policy. We believe this approach strikes the most effective balance between\nflexibility, robustness, and interpretability for our specific research questions.\n2. Data Heterogeneity R2-3: Thank you for this critical suggestion regarding the potential impact of data heterogeneity and\nThe authors work with a large, structural breaks. To address your concern, we have undertaken a two-pronged approach: first, by adding a\nheterogeneous panel of 139 countries. more transparent discussion of data limitations in the manuscript, and second, by conducting a comprehensive\nWhile fixed-effects models account for robustness check where we exclude countries that experienced severe structural breaks.\ntime-invariant heterogeneity, there is\nlittle discussion on how missing data, 1. Transparent discussion of data limitations:\nreporting inconsistencies, or structural We have added a new paragraph in the \"Limitations and future research\" section to explicitly\nbreaks (e.g., due to war, institutional acknowledge the limitations you mentioned, including reporting inconsistencies and the impact of major\nreforms) might bias results. This reforms.\nshould be discussed more\ntransparently, possibly supported by a New text in the \"Limitations and future research\" section (Lines 378-399):\nrobustness check (e.g., excluding \"While this study provides a comprehensive global analysis... several limitations should be acknowledged...\ncountries with major conflicts, or First, our analysis relies on country-level panel data. This macro-level approach... inevitably masks\nsensitivity to reporting gaps). significant sub-national heterogeneity... Second, due to data limitations across the long time series, we\nemployed proxy variables for key concepts... Third, while our dynamic panel models account for key\neconometric challenges, they cannot fully eliminate all potential endogeneity concerns...\"\n2. Comprehensive robustness check excluding countries with major structural breaks:\nTo formally test whether our results are driven by countries with extreme historical events, we conducted a\nrigorous robustness check. We identified 16 countries/entities that underwent severe, prolonged armed\nconflicts or state dissolution, representing major structural breaks. We then re-estimated our entire set of main\nanalyses on a restricted sample excluding these countries.\nWe have created a dedicated new section in the Supplementary Information (Section S3) to present this\nentire analysis. We introduce it as follows:\nNew text in Supplementary Information, Section S3.1 (Lines 269-277):\n\"To address the concern that our main findings might be disproportionately influenced by a small number of\ncountries that experienced exceptionally severe and prolonged historical disruptions, we conducted a\nrigorous robustness check. The inclusion of countries that underwent major structural breaks due to large-\nscale conflicts or state dissolution during the study period (1961–2021) could potentially introduce bias.\nSuch events can cause systemic shocks to agricultural systems and data reporting that are unrelated to the\ntypical processes of urbanization and agricultural modernization we aim to study.\"\nThe results of this extensive check, now presented in Supplementary Figures S4-S6 and Tables S6-S9,\noverwhelmingly confirm the stability of our main findings. We summarize the key outcomes below:\nRobustness of transmission pathways (corresponding to Fig. 2): The re-estimated path analysis, now\nshown in Supplementary Fig. S4, remains highly consistent with our main findings.\nFrom the discussion of the robustness check (Supplementary Information, Lines 292-311):\n\"First, the re-estimated dynamic path analysis reveals that the key transmission pathways, their signs, and\ntheir statistical significance levels remain largely consistent with the main analysis (Fig. S4; full results in\nTables S6-S7). For instance, the negative association between farm size and machinery intensity in low-\nincome countries and the positive association in high-income countries both remain statistically significant.\nThis confirms that the identified structural differences... are a general feature of the data, not an artifact of\noutlier countries.\"\nRobustness of nonlinear trajectories (corresponding to Fig. 3): The characteristic nonlinear curves, re-\nestimated and presented in Supplementary Fig. S5, are remarkably stable.\nFrom the discussion of the robustness check (Supplementary Information, Lines 313-319):\n\"Second, the re-estimation of our main dynamic panel interaction model shows that the identified nonlinear\npatterns and moderating roles are highly robust. The characteristic U-shaped and inverted U-shaped\nrelationships... persist across all income groups and modernization scenarios (Fig. S5). The locations of the\ncurve vertices (turning points) are remarkably stable; for example, the turning point for the U-shaped yield\ncurve in the high-modernization scenario for low-income countries shifted only slightly from 34% to 38%.\"\n3. Environmental Heterogeneity and R2-4: Thank you for this crucial point regarding potential omitted variable bias from time-varying\nOmitted Variable Bias environmental factors. To address this, we have implemented a two-layered enhancement for the climate\nWhile the authors include a rich set of variables, which we believe is more targeted to the issue of time-varying shocks than the use of agro-\ncontrols, including PGDP, crop type, ecological zone dummies. We have also acknowledged the limitation regarding dynamic soil conditions.\nagricultural investment, and fixed\neffects at the country and year level, the 1. Integration of climate controls and robustness check for climate anomalies:\nmodel may still omit time-varying 1) We have now updated all baseline specifications to include annual average precipitation and\nenvironmental variables such as temperature as core time-varying controls.\ninterannual climate variability (e.g.,\nrainfall, temperature anomalies) and This revision is now reflected in the Methods section (Lines 564-566):\ndynamic soil conditions, which are \"W is a vector of time-varying control variables, including... precipitation, and temperature, along with\nit\nknown to affect nitrogen use efficiency their quadratic terms...\"\nand yield response (e.g., Ray et al.,\n2015; Leng et al., 2016; Zhang et al., 2) To specifically address \"short-term agro-environmental shocks\", we conducted a stringent robustness\n2015). The use of fixed effects absorbs check using \"climate anomaly\" variables. This entire analysis is presented in Supplementary\naverage differences across countries, Information (Section S4), and the results (shown in Figs. S7-S9 and Tables S11-S14) overwhelmingly\nbut it does not address short-term agro- confirm that our core findings are robust to interannual climate shocks.\nenvironmental shocks. A robustness\ncheck incorporating agro-ecological The methodology is described in Section S4.1 (Lines 438-447):\nzone dummies or historical climate \"To address potential omitted variable bias... we conducted a second robustness check... To account for these\nindices could help assess whether shocks more explicitly, we constructed two 'climate anomaly' variables (annual deviations from the country's\nurbanization effects are overstated due long-term average precipitation and temperature). In this robustness check, we replaced the standard climate\nto omitted variable bias. control variables... with these two climate anomaly variables...\"\nThe results of this check are highly compelling and demonstrate that our core findings are remarkably robust,\neven when specifically controlling for short-term climate shocks.\nThe conclusion of this robustness check is summarized in Section S4.2 (Lines 450-455 & 487-489):\n\"The results of this robustness check provide strong evidence that our core findings are robust and not biased\nby the omission of interannual climate variability. The re-estimated dynamic path analysis shows that the\nfundamental transmission mechanisms are not confounded by short-term climate shocks, as the pathways\nremain highly consistent with the main analysis (Fig. S7...).\n\"In conclusion, the relationships we have identified are a robust feature of the data, independent of annual\nclimate fluctuations, which strengthens the credibility of our findings and the policy implications derived\nfrom them.\"\nWe chose this approach over using agro-ecological zone (AEZ) dummies for a critical reason. AEZ\nclassifications are typically time-invariant. Including them as dummies would mean their effect would be\nabsorbed by the country fixed effects already in our model, as most countries do not change their dominant\nAEZ over time. Therefore, AEZ dummies would not help control for the time-varying nature of climate\nshocks (e.g., a specific drought year), which was the core of your concern. Our use of annual climate anomaly\nvariables, however, directly captures these year-to-year environmental shocks within each country, making it\na more appropriate and powerful test for this specific type of omitted variable bias.\n2. Regarding dynamic soil conditions:\nWe also wish to address your insightful point about dynamic soil conditions. We fully acknowledge that\nfactors like changes in soil organic matter, nutrient depletion, or degradation over time can influence nitrogen\ndynamics. Unfortunately, compiling a consistent, internationally comparable panel dataset for dynamic soil\nproperties spanning 139 countries over 60 years is not currently feasible due to severe data limitations and\nlack of standardized measurement protocols across nations and time.\nWe have taken steps to mitigate this issue to the extent possible within our framework:\nCountry fixed effects: Our model's country fixed effects control for all time-invariant soil characteristics,\nsuch as inherent soil type, texture, and baseline fertility, which account for a significant portion of the cross-\ncountry variation in soil properties.\nAcknowledged as a limitation: Recognizing that we cannot fully capture the time-varying component of\nsoil health, we have explicitly included this as a limitation and an avenue for future research in our revised\nmanuscript.\nNew text in the \"Limitations and future research\" section (Lines 399-405):\n\"...our model does not explicitly control for time-varying soil conditions. While the country fixed effects in\nour model account for time-invariant, baseline soil characteristics (e.g., soil type and texture), dynamic\nchanges in soil health over the 60-year period, such as degradation or improvements in soil organic matter,\ncould also influence nitrogen dynamics. Constructing a consistent, multi-decade global panel dataset for\nsuch dynamic soil properties remains a major challenge, but future research could explore this important\ndimension, perhaps at smaller regional scales where such data may be available.\"\n4. Y-Axis Scale Consistency for R2-5: Thank you for this excellent and very practical suggestion. You are absolutely right that inconsistent\nComparative Interpretation Y-axis scales hindered the direct visual comparison of effect magnitudes across income groups.\nIn Figures 5 and Extended Data Fig. 3,\nthe Y-axis scales differ across income We have now revised all relevant figures to ensure that for any given outcome variable (e.g., nitrogen\ngroups for the same indicators, which fertilizer application, nitrogen yield, NUE), the Y-axis scale is standardized and consistent across all panels\nmakes it difficult to visually compare representing the different income groups. This change has been implemented in our main conditional effect\nthe magnitude of urbanization-related plots, which correspond to the figures you mentioned in the original manuscript:\neffects across low-middle, upper-\nmiddle, and high-income countries. For example, Fig. 3 (which presents the main nonlinear effects under low vs. high modernization scenarios).\nTo enhance interpretability and\nfacilitate fair cross-group comparison,\nwe recommend adopting consistent Y-\naxis scales across the three income\ngroups for each indicator. This\nstandardization would allow readers to\nmore clearly assess the relative impact\nof urbanization across contexts and\nsupport one of the central aims of the\nstudy: to understand differential\nurbanization effects globally.\nFig. 3 | Agricultural modernization reshapes the nonlinear effects of urbanization on nitrogen fertilizer\napplication, nitrogen yield and NUE.\nReferee #3\nComment Response and Change\nI co-reviewed this manuscript with one R3-1: Thank you for letting us know about the co-review process. We appreciate the thorough attention that\nof the reviewers who provided the your collaborative review has brought to our manuscript.\nlisted reports. This is part of the Nature\nCommunications initiative to facilitate\ntraining in peer review and to provide\nappropriate recognition for Early\nCareer Researchers who co-review\nmanuscripts.\nReferee #4\nComment Response and Change\nThis study offers a global-scale analysis of the complex interactions among R4-1: Thank you for your insightful and clear summary of the key issues\nurbanization, farm size, technological advancement, and irrigation in our manuscript. We sincerely appreciate your valuable feedback. We\ninfrastructure in shaping nitrogen use and agricultural productivity. The will now address your specific comments in the following sections.\nfindings highlight the importance of region-specific management strategies\nthat account for the nonlinear effects of urbanization, and provide targeted\nrecommendations for nitrogen management across different stages of urban\ndevelopment. However, some aspects—particularly regarding summary of\nkey findings, indicator selection, and the interpretation of results—would\nbenefit from further analysis, clarification, and refinement. My comments\nare as follows:\n1. Although the data and modeling R4-2: Thank you for this critical feedback. We agree completely that our manuscript must move beyond\nwork support the points made in the simply reinforcing existing arguments and clearly articulate its novel contributions. To achieve this, we have\nsecond paragraph of the substantially revised the Introduction. Instead of merely reviewing the dualistic (positive vs. negative) views\nIntroduction, the findings primarily in the literature, we now argue that the true knowledge gap lies in understanding the conditional and\nreinforce these existing arguments heterogeneous nature of urbanization's effects.\nrather than delivering substantially\nnew perspectives. This reduces the The revised Introduction now frames our study's novelty by asking (Lines 53-61):\nperceived depth and novelty of the \"A critical knowledge gap exists in understanding this conditional heterogeneity: how does the association\nresults. I suggest the authors more evolve with socioeconomic development in a country, and how is this evolution linked to the pace of its own\nexplicitly analyze and articulate the agricultural modernization? Specifically... how do the resulting changes in farm size, machinery intensity,\nnovel aspects of the specific and irrigation intensity act as moderating variables...? Answering these questions requires moving beyond a\nfindings derived from the simple net-effect estimation to a more dynamic, interactive framework...\"\ncalculation.\nBy reframing the study in this way, our findings deliver several substantially new perspectives, which we\nhave now sharpened and explicitly articulated throughout the manuscript.\n1. We move from 'What' to 'How' and 'When': Our study provides the first systematic global\nquantification of how and when the roles of modernization variables change. Our heatmap analysis (Figure\n4) pinpoints that for low-income countries, farm size is the key moderator in early urbanization, whereas for\nupper-middle-income countries, irrigation's role is dominant early on, but machinery becomes the sole driver\nof yield gains in late urbanization.\n2. We reveal divergent and opposing pathways, not a single path: Our central novel finding is that there\nis no universal pathway. We have highlighted this core contribution directly in the abstract.\nAs stated in the revised Abstract (Lines 28-32):\n\"These results indicate there is no 'one-size-fits-all' pathway for aligning urbanization with nitrogen\nsustainability. Navigating divergent pathways by strategically integrating land consolidation,\nmechanization, and precision irrigation can transform urbanization from a challenge into a catalyst for\nsustainable nitrogen management and resilient food systems.\"\nFurthermore, we provide robust evidence for this divergence. For instance, we show that for low-income\ncountries, farm size expansion is a crucial buffer. In stark contrast, for high-income countries, we find that\nfurther farm size expansion can be detrimental, and resilience depends on technology.\nThis finding is now explicitly discussed in the Results and discussion section (Lines 305-310):\n\"This nuanced finding refines the average positive role of farm size by delineating its limits. For high-income\neconomies, sustainable nitrogen management in the face of continued urbanization will depend less on\nfurther farm consolidation and more on technology-driven intensification, emphasizing mechanization and\nirrigation efficiency to balance productivity, resource use, and environmental resilience.\"\n3. We quantify critical policy trade-offs: Our analysis provides novel, quantitative evidence of the specific\ntrade-offs policymakers face. This is a key theme of our revised discussion.\nFor example, we now articulate this trade-off for upper-middle-income countries in the Results and\ndiscussion section (Lines 287-293):\n\"This divergence underscores a key policy challenge for middle-income economies: distinct modernization\nstrategies carry distinct trade-offs. Policies prioritizing irrigation for environmental benefits may sacrifice\nnitrogen yield, while emphasizing consolidation may undermine efficiency. Effective navigation of this\ntransitional stage requires balanced modernization portfolios, such as coupling irrigation investments with\nmechanization to sustain nitrogen productivity while maintaining environmental gains.\"\nThis quantification of competing pressures on policy levers, supported by the empirical results in Fig. 4,\nis a new and actionable insight.\n2. Regarding indicator selection, the R4-3: Thank you for raising this critical methodological point. You have identified a potential issue of\ncalculations of farm size, \"endogeneity by definition\" that could create a spurious correlation between our urbanization variable and\nmechanization ownership, and the mediating variables. To resolve this issue at its root, we have fundamentally revised the definitions of\nirrigated area per capita all use rural all three agricultural modernization indicators, completely breaking the definitional link with the \"rural\npopulation as the denominator. population\" denominator. This substantial revision not only solves the endogeneity problem but also\nSince urbanization level is also enhances the theoretical soundness of our indicators. These changes are now detailed in the \"Data sources\ndetermined by the proportion of and processing\" subsection of the Methods.\nrural population, these indicators\nmay have endogenous relationships 1. For 'farm size': From rural population to agricultural labor\nwith urbanization itself. Does this We have replaced the denominator \"rural population\" with \"agricultural labor force,\" which is a more precise\nimply that, all else being equal, measure of the active workforce engaged in farming.\nincreases in urbanization level are\nnecessarily associated with Revision in the Methods section (Lines 430-435):\nincreases in farm size, \"Farm size was estimated as the cropland area per agricultural laborer. This metric provides a more direct\nmechanization, and irrigation? measure of agricultural operational scale by linking land resources to the active agricultural workforce,\nHow do the authors account for this thereby avoiding potential endogeneity that could arise from using total rural population as a denominator.\npotential endogeneity and its Data for the agricultural labor force were obtained from the FAOSTAT database.\"\ninfluence on the results?\n2. For 'machinery' and 'irrigation': From per capita to intensity-based metrics\nTo completely eliminate any population-based denominator, we have shifted from \"per capita\" metrics to\n\"intensity-based\" metrics (i.e., per unit of cropland area). This provides a more direct measure of\ntechnological adoption on the land itself.\nRevision for Machinery in the Methods section (Lines 456-458):\n\"Machinery intensity... was measured as total agricultural machinery power per unit of cropland area (CV\nha⁻¹). This approach directly quantifies the density of capital investment in machinery on cropland,\nindependent of demographic variables.\"\nRevision for Irrigation in the Methods section (Lines 468-471):\n\"Similarly, irrigation intensity was defined as the proportion of the 'actually irrigated area under full control\nirrigation' relative to the total cropland area (%). This intensity-based indicator reflects the extent of water\nmanagement technology adoption across cropland.\"\nThrough these systematic revisions, we have ensured that our key mediating variables are now conceptually\nand mathematically independent of the urbanization rate's definition. We are pleased to report that after re-\nrunning all analyses with these more rigorous indicators, our core findings on the nonlinear and\nheterogeneous effects of urbanization remain highly robust. All figures and tables in the manuscript have\nbeen updated accordingly.\n1. The first two sentences of the R4-4: Thank you for this excellent suggestion. We agree that the opening sentences of our manuscript must\nabstract involve multiple concepts. establish a clear and logical framework. We have revised both the abstract and the introduction to make this\nSpecifically, the first sentence logical relationship transparent.\nstates that urbanization impacts\nland use and population 1. Clarifying the logical flow in the Abstract and Introduction:\ndistribution, while the second We have rewritten the opening of the Abstract to establish a more direct and logical pathway. Instead of\nemphasizes a particular focus on listing concepts, we now frame them as an interconnected process.\nmechanization and irrigation. The\nlogical relationships among these The revised opening of the Abstract now reads (Lines 15-17):\nkeywords are not clear. I\nrecommend refining these \"Urbanization reshapes agricultural systems by altering labor and land-use, presenting a critical interplay\nsentences to make them more with agricultural modernization processes such as farm size expansion, machinery, and irrigation.\"\nfocused and concise, and explicitly\ndefining the dimensions of This revised sentence clarifies the logical flow: Urbanization → alters labor and land-use → which creates\nurbanization addressed in this an interplay with modernization (farm size, machinery, irrigation). This directly addresses your concern about\nstudy. The same issue exists in the the unclear relationship between these keywords. A similar logical structure is now established at the\nfirst three sentences of the main beginning of the Introduction.\ntext. In other words, in lines 62-63,\nwhen stating that urbanization 2. Specifying the socioeconomic dimensions of urbanization:\n“accounts for... socioeconomic To address your second point, we have now explicitly defined the key socioeconomic dimensions through\ndimensions,” it would be helpful to which we analyze urbanization's impact. We specify that our study focuses on how urbanization, via rural\nspecify which dimensions are being labor migration and land consolidation, drives changes in agricultural structure.\nconsidered.\nThis is now specified in the Introduction (Lines 55-58):\n\"Specifically, as urbanization coincides with rural labor migration and land consolidation, how do the\nresulting changes in farm size, machinery intensity, and irrigation intensity act as moderating variables that\ncan either amplify or buffer the outcomes associated with urbanization at different stages...?\"\nBy explicitly naming \"rural labor migration and land consolidation\" as the core mechanisms, we have\nclarified the specific dimensions being considered, making the link between the macro trend of urbanization\nand our variables of interest (farm size, etc.) much more direct and theoretically grounded.\n2. The results section of the abstract R4-5: Thank you for this valuable suggestion. We have substantially revised the abstract to integrate specific,\nwould be strengthened by including impactful quantitative data, directly addressing your request for details on significant turning points and\nkey quantitative findings, changes in NUE. Below is the revised abstract.\nespecially regarding changes in\nNUE or significant thresholds We have used bold text to highlight the key quantitative findings that have been added in response to your\nidentified in the study. comment:\nRevised Abstract with key quantitative findings highlighted (Lines 15-32) :\n\" Urbanization reshapes agricultural systems by altering labor and land-use, presenting a critical interplay\nwith agricultural modernization processes such as farm size expansion, machinery, and irrigation. This\ninterplay drives nonlinear patterns in cropland nitrogen use. Using a 61-year panel dataset across 139\ncountries, we show that the relationship between urbanization and agricultural nitrogen outcomes is\nprofoundly nonlinear and contingent on economic development and modernization levels. In low-income\ncountries, urbanization initially intensifies fertilizer use but suppresses nitrogen yield and efficiency, with\nlarger farm size buffering early losses. As countries transition to upper-middle-income status, agricultural\nmodernization reconfigures these associations, enhancing efficiency but introducing trade-offs between\nenvironmental gains and yield growth. In high-income countries, modernization mitigates the adverse\nimpacts of advanced urbanization, though diminishing returns to scale emerge. At a 90% urbanization level,\na high degree of modernization reverses the impact on nitrogen use efficiency from a negative 4% to a\npositive 12% compared to low modernization. These results indicate there is no 'one-size-fits-all' pathway\nfor aligning urbanization with nitrogen sustainability. Navigating divergent pathways by strategically\nintegrating land consolidation, mechanization, and precision irrigation can transform urbanization from a\nchallenge into a catalyst for sustainable nitrogen management and resilient food systems.”\n3. Lines 53-56: What specific R4-6: Thank you for this excellent point. To make this rationale clear from the outset, we have substantially\ncharacteristics of urbanization in revised the Introduction. Instead of simply reviewing the literature, we now explicitly lay out the specific\nlow-income countries (as defined characteristics and mechanisms we hypothesize to be different between low- and high-income countries,\nby the conceptual framework providing a strong theoretical motivation for our stratified analysis.\nabove) lead to potentially different\nimpacts on nitrogen use compared We have rewritten the final paragraphs of the Introduction to detail these expected differences.\nto high-income countries? A Excerpt from the revised Introduction (Lines 53-71):\ndiscussion of this point would help \"A critical knowledge gap exists in understanding this conditional heterogeneity: how does the association\nreaders better understand the evolve with socioeconomic development in a country, and how is this evolution linked to the pace of its own\nrationale for conducting a global agricultural modernization? Specifically, as urbanization coincides with rural labor migration and land\nanalysis and the importance of consolidation, how do the resulting changes in farm size, machinery intensity, and irrigation intensity act as\npaying special attention to the moderating variables that can either amplify or buffer the outcomes associated with urbanization at different\ntransition from low- to high- stages?\nincome contexts.\nTo fill this gap, this study aims to disentangle these complex, interactive relationships. We hypothesize that\nthe association of urbanization… critically contingent upon the stage of economic development and its\nprevailing levels of farm size, machinery intensity, and irrigation intensity in a country. In low-income\ncountries, the primary challenge is whether farm size can grow fast enough to compensate for labor loss\nand create a viable scale for machinery intensity. In contrast, in high-income countries, the challenges\nshift to managing the trade-offs of a highly optimized system, questioning whether further advances in\nmachinery intensity and irrigation intensity can provide resilience against diminishing returns.\"\nAs this revised text shows, we now explicitly argue that the \"specific characteristics\" leading to different\nimpacts are:\n1) The prevailing agricultural structures (e.g., smallholder vs. large-scale).\n2) The nature of the developmental challenge (e.g., overcoming initial productivity loss vs. managing\ndiminishing returns).\nOur new introduction frames the entire study as a direct test of these hypotheses. Furthermore, our Results\nsection now systematically presents the empirical evidence confirming this initial framework (e.g.,\ndemonstrating a U-shaped yield curve for low-income countries versus an inverted U-shaped curve for high-\nincome countries).\n4. Lines 89-91 state that the decline in R4-7: Thank you for this insightful question. Upon re-examining the data and relevant literature, we confirm\nEastern Europe potentially stems that the most significant decline occurred primarily during the 1990s, following the dissolution of the Soviet\nfrom its naturally high soil fertility Union. As you rightly suggested, this downturn was predominantly driven by profound socioeconomic\nand extensive black soils. I am transformations.\nunsure about the specific timeframe\nduring which this decline occurred. We have substantially revised the corresponding paragraph in the Results and discussion section to reflect\nCould this decline also be related to this more accurate and comprehensive explanation.\nsocioeconomic developments, such\nas economic transformation, Revised text in the manuscript (Lines 96-101):\nagricultural decline, or \"A notable deviation occurred in Eastern Europe, where nitrogen fertilizer application fell sharply in the\nenvironmental policies? 1990s after the dissolution of the Soviet Union (Fig. 1e, h). Despite favorable agricultural conditions\nprovided by extensive black soils and high natural fertility, the collapse of state-subsidized systems, rising\ninput costs, and disrupted production and supply chains during the economic transition led to a substantial\nreduction in nitrogen fertilizer application.”\nThis revision now explicitly states the timeframe (the 1990s) and correctly attributes the decline to the major\nsocioeconomic shifts you mentioned, such as the collapse of state-subsidized systems and supply chain\ndisruptions during the economic transition. We have supported this revised explanation with appropriate\ncitations to the literature on transition economies.\nReferences:\n1. Liefert, W. M., & Swinnen, J. F. M. Changes in Agricultural Markets in Transition Economies.\nAgriculture Information Bulletin No. 709. USDA, Economic Research Service. (2002)\n2. Macours, K. & Swinnen, J.F.M. Impact of Initial Conditions and Reform Policies on Agricultural\nPerformance in Central and Eastern Europe, the Former Soviet Union, and East Asia. Am. J. Agr. Econ. 82,\n1149-1155 (2000).\n5. I suggest that the authors include R4-8: This is an excellent suggestion. We have now revised the Results and discussion section to include\nthe analysis of typical countries typical country cases that exemplify the patterns we describe.\nwithin the results section to support\nand illustrate the main conclusions. 1. Example for yield stagnation in Central Asia:\nFor example, in lines 101-103, Following your advice, we have replaced the vague term \"some countries\" with specific examples and\nwhen referring to \"some countries,\" provided a more nuanced explanation, linking their experience to the broader post-Soviet socioeconomic\nit would be helpful to specify which transition.\ncountries are involved and describe\ntheir experiences with unfavorable Revised text in the manuscript (Lines 110-115):\nclimatic conditions, soil \"…Central Asia, particularly countries such as Kazakhstan and Uzbekistan, experienced stagnation and\ndegradation, or limited access to even sharp declines in nitrogen yield during the 1990s (Fig. 1j–l). The breakdown of Soviet-era agricultural\ncritical agricultural inputs. This systems, difficulties in transitioning to market economies, and restricted access to essential inputs severely\nwould help readers gain a clearer constrained productivity. Environmental stressors such as arid climates, soil degradation, and long-standing\nunderstanding of the issues. irrigation challenges further amplified these declines.\"\n2. Examples for divergent trends in NUE:\nSimilarly, when discussing divergent NUE trends, we have now included specific country examples to make\nthe comparison more vivid and grounded.\nRevised text in the manuscript (Lines 130-137):\n\"Some nations achieved higher nitrogen yields while reducing nitrogen fertilizer application. For example,\nDenmark and the Netherlands have substantially improved NUE through precision farming and optimized\nnutrient management. Conversely, many rapidly industrializing countries such as China and India\nexperienced declining NUE during periods of input-driven intensification, as nitrogen fertilizer application\nexpanded faster than nitrogen yield gains.\"\n6. Please check the numerical data R4-9: Thank you for your careful reading and for flagging this for our attention. The sentence in the original\npresented in line 197. manuscript (line 197) contained numerical data related to a specific threshold identified by our previous\nsegmented regression approach. As detailed in our responses to other major comments (e.g., R2-2), we have\nundertaken a fundamental upgrade of our entire empirical strategy. We have replaced the segmented\nregression framework with a more flexible and robust continuous interaction model. This new approach does\nnot rely on identifying discrete thresholds.\nAs a result of this methodological shift, the entire section discussing specific threshold values has been\nremoved and rewritten. Consequently, the original sentence and the numerical data you questioned no longer\nexist in the revised manuscript.\n7. During the study period, some R4-10: Thank you for this insightful question. We have made two major revisions: first, by clarifying our\ncountries experienced significant categorization method and its rationale in the Methods section, and second, by using the concept of\nincome level changes (e.g., \"transitioning economies\" to frame our entire Policy Implications section, as you astutely suggested.\ntransitioning from low- to middle-\nincome). How are these countries 1. Clarification of the country categorization method:\ncategorized in the temporal We have now explicitly stated our methodological choice to use a fixed income-group classification and\nanalysis? Given that the authors explained its rationale.\nmention the transition from low- to New text in the Methods section (Lines 575-584):\nhigh-income contexts in the \"To systematically examine heterogeneity across development stages, we partitioned the full sample into\nintroduction, could the authors use three sub-groups based on World Bank income classifications... This classification is based on the World\ntypical examples of countries that Bank's income group designations for the most recent year of our dataset (2021). This approach ensures\nhave undergone such transitions to that each country remains in the same group throughout the analysis, allowing us to interpret the results\nfurther elucidate their findings? as the average developmental trajectory for countries that have achieved a certain income status. We\nacknowledge that some countries transitioned between groups during the 60-year study period, and we\nexplore the implications for these transitioning economies in our Policy Implications section.\"\nThis clarification explains that our goal is to compare the structural patterns of countries that are currently at\ndifferent development stages, a standard and interpretable approach for this type of analysis.\n2. Analysis of transitioning countries in the Policy Implications section:\nFollowing your excellent suggestion, we have completely restructured the Policy Implications section to\nbe a dynamic discussion framed around the journey of transitioning economies. We now use our cross-\nsectional findings to offer longitudinal insights, using typical transitioning countries as concrete examples.\nExcerpt from the revised \"Policy Implications\" section (Lines 337-349):\n\"As economies advance toward upper-middle-income status, the policy focus must shift from foundational\ncapacity building to balancing environmental sustainability with continued productivity gains. Our analysis\nidentifies irrigation as a powerful lever for improving nitrogen efficiency and reducing nitrogen fertilizer\napplication, while machinery intensity emerges as a key determinant of late-stage nitrogen yield growth. The\noptimal balance between these levers depends heavily on each country’s farm structure. For smallholder-\ndominated systems such as China and Vietnam, policy evolution has followed a two-stage trajectory, from\nearly market liberalization and infrastructure expansion to contemporary efforts targeting nitrogen surpluses\nand water scarcity. Similar economies must integrate environmental and food security objectives within a\nsingle policy framework, ensuring that ecological limits are internalized into growth strategies. Conversely,\ncountries such as Brazil and Mexico, characterized by dualistic agricultural structures dominated by large\ncommercial farms, face a different challenge: ensuring the sustainability of already modernized systems.\"\nFurther excerpt on countries transitioning to high-income status (Lines 356-367):\n\"In high-income countries and those transitioning into this category, the central policy concern shifts from\noptimization to resilience. Our results show that at very high levels of urbanization, even advanced\nagricultural systems encounter diminishing returns, with potential declines in nitrogen yield and NUE. At\nthis stage, expanding farm size is no longer a viable pathway to sustainability; further consolidation can\ngenerate negative externalities, including biodiversity loss, soil degradation, and social inequities. Japan\nexemplifies this inflection point: a highly mechanized, small-scale system now constrained by severe labor\nshortages and declining marginal returns, prompting a strategic pivot toward technological innovation and\nresilience building. Similarly, Poland’s post-socialist agricultural transformation, guided by the EU’s\nCommon Agricultural Policy (CAP), reflects a decisive policy reorientation toward environmental\nsustainability and system stability. For these advanced economies, the policy priority should be technology-\ndriven intensification rather than scale expansion…\"\n8. Please discuss the limitations of R4-11: Thank you for this essential suggestion. We have added a new, comprehensive section\nthis study and suggest directions for titled \"Limitations and future research\" at the end of the manuscript (Lines 471-499). This section\nfuture research. transparently discusses the key limitations of our study, many of which reflect the insightful points you raised\nin your other comments, and proposes several avenues for future inquiry.\nThe new section is structured around five key points:\nSummary of the new \"Limitations and future research\" section (Lines 378-412):\n1. Country-level Data: \"First, our analysis relies on country-level panel data... [which] inevitably masks\nsignificant sub-national heterogeneity...\" (Suggests future research using grid-level data).\n2. Proxy Variables: \"Second, due to data limitations... we employed proxy variables for key concepts such\nas 'farm size'... it remains an approximation...\" (Suggests future studies using detailed agricultural census\ndata).\n3. Endogeneity: \"Third, while our dynamic panel models account for key econometric challenges, they\ncannot fully eliminate all potential endogeneity concerns...\" (Suggests future research using quasi-\nexperimental methods).\n4. Dynamic Soil Conditions: \"Fourth, our model does not explicitly control for time-varying soil\nconditions... Constructing a consistent, multi-decade global panel dataset for such dynamic soil properties\nremains a major challenge...\" (Suggests future research at smaller regional scales).\n5. Focus on the Nitrogen Cycle: \"Finally, our study focuses on the nitrogen cycle... Future research should\nextend this framework to analyze... other critical dimensions of sustainability, such as water use, greenhouse\ngas emissions, biodiversity, and rural livelihoods...\"\nAs summarized above, this new section addresses limitations related to data granularity, variable\nmeasurement, model identification, and the scope of our analysis. For each point, we have proposed a\nconcrete direction for future research. Notably, we have integrated the discussions on proxy variables,\nendogeneity, and omitted environmental variables (soil conditions) that arose from your other valuable\ncomments, framing them as broader limitations and forward-looking research questions.\nReferee #5\nComment Response and Change\nI co-reviewed this manuscript with one R5-1: Thank you for letting us know about the co-review process. We appreciate the thorough attention that\nof the reviewers who provided the your collaborative review has brought to our manuscript.\nlisted reports. This is part of the Nature\nCommunications initiative to facilitate\ntraining in peer review and to provide\nappropriate recognition for Early\nCareer Researchers who co-review\nmanuscripts."}, {"round": 2, "reviewer_comments": {"reviewer_2": "Comments on the Revised Manuscript “Interplay of urbanization and agricultural modernization shapes nitrogen use in global croplands”\nJournal: Nature Communications\nThe revised manuscript has substantially improved upon the original version, particularly in clarifying terminology, redefining key variables, strengthening methodological rigor, and addressing the appropriateness of the statistical framework. The replacement of the original cross-sectional SEM with dynamic panel path analysis and interaction models represents a major methodological upgrade. Nevertheless, several conceptual, methodological, and interpretational issues remain.\nThese issues do not invalidate the analysis but require further clarification and refinement before the work can reach the level of robustness and transparency. Below I detail both the strengths and remaining weaknesses.\n1. The revised manuscript represents a major methodological improvement: whereas the first version relied on a conventional SEM framework unsuited for 61-year panel data, the authors now employ a system of dynamic two-way fixed-effects models with lagged dependent variables, effectively addressing concerns about autocorrelation, unobserved heterogeneity, and non-stationarity. However, several issues remain. First, because each path is estimated separately, the analysis does not account for potential cross-equation dependence or correlated errors, which may bias derived indirect effects; a brief acknowledgement of this limitation would improve transparency. Second, the manuscript should more clearly explain the interpretational role of the lagged dependent variables, which materially affect both the magnitude and meaning of the estimated coefficients. Third, although a robustness check excluding countries with major structural breaks has been added, its implications are not adequately summarized in the main text; this is important because several key nonlinear patterns, such as the U-shaped yield response, are sensitive to sample composition.\n2. The revised manuscript offers a clearer and better-justified treatment of the farm-size proxy: the authors replace the previously inconsistent use of per-capita cropland area with cropland area per agricultural laborer, provide a precise definition, and move the empirical validation into both the main text and the supplement, where they report a strong correlation with actual farm size. Terminology is now consistently applied, and the rationale for preferring this proxy over census-based measures is more clearly articulated. Nonetheless, several conceptual concerns remain. Although correlated with observed farm size, cropland per laborer fundamentally captures labor availability rather than the true operational scale of holdings, and this distinction should be explicitly acknowledged as a limitation. In addition, the theoretical framework implicitly assumes that farm size, mechanization, and irrigation evolve jointly, yet in many regions with fragmented land tenure, such as parts of South Asia, mechanization may increase independently of consolidation; the manuscript would benefit from a more nuanced discussion of such divergences.\n3. While the revised introduction is more coherent and focused, the overall novelty of the study still requires clearer articulation. Prior global research has examined urbanization, farm structure, and nitrogen trends largely in isolation, and the distinctive contribution of this work lies in integrating these domains within a nonlinear, moderating framework, yet this message is at times obscured by extended descriptive background on global trends. The authors should sharpen the narrative by explicitly emphasizing why capturing joint nonlinearities yields insights that neither linear global models nor regional case studies can provide, and by positioning this integrative approach as the central conceptual innovation of the manuscript.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "I appreciate the authors' substantial revisions in response to my previous comments. These changes have significantly enhanced the clarity and impact of the manuscript. I have only a few minor comments remaining:\n1. Regarding my previous comment 1, I note clear improvements in the revised summary of existing literature and the identified knowledge gap. To further strengthen the paper, I suggest that the authors consider referencing relevant findings or discussions in the Results and Discussion section. It would be beneficial to highlight how this study responds to, or extends, previous work or discussions on the relationship between urbanization and nitrogen use.\n2. Regarding the selection of agricultural modernization indicators, I appreciate both the substantial revisions and the rationale provided. I recommend including a brief discussion, perhaps in the second point of the limitation section, on the potential biases associated with these proxy variables and how they may affect the interpretation of farm size expansion, mechanization, and irrigation. For example, what relevant aspects might not be fully captured, and how could this influence the results? Additionally, please note the repetition of \"Machinery intensity\" in Line 456.\n3. I encourage the authors to reconsider the title of the second Results section, \"Average transmission pathways of urbanization,\" as its meaning is not entirely clear. If the authors prefer to retain this title, I suggest providing a brief explanation at the beginning of the section to clarify what is meant by \"Average transmission pathways.\"", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to Referees’ Comments\nReferee #2\nComment Response and Change\nThe revised manuscript has substantially improved upon the original version, R2-1: We are very grateful for your encouraging feedback.\nparticularly in clarifying terminology, redefining key variables, strengthening\nmethodological rigor, and addressing the appropriateness of the statistical framework. In this current revision, we have carefully addressed each of\nThe replacement of the original cross-sectional SEM with dynamic panel path analysis the remaining conceptual, methodological, and\nand interaction models represents a major methodological upgrade. Nevertheless, several interpretational points you raised. We have provided further\nconceptual, methodological, and interpretational issues remain. These issues do not clarifications and expanded our discussion on model\ninvalidate the analysis but require further clarification and refinement before the work limitations and proxy variables to ensure the manuscript\ncan reach the level of robustness and transparency. Below I detail both the strengths and reaches the desired level of robustness and transparency.\nremaining weaknesses.\n1. The revised manuscript represents a major R2-2: Thank you for this rigorous technical insight. In the interest of transparency and\nmethodological improvement: whereas the first methodological rigor, we have added an explicit acknowledgement of this limitation in the\nversion relied on a conventional SEM framework\n\"Limitations and future research\" section.\nunsuited for 61-year panel data, the authors now\nemploy a system of dynamic two-way fixed-effects\nChange: We have added the following sentences to the third point of the \"Limitations and\nmodels with lagged dependent variables, effectively\nfuture research\" section (Line 434):\naddressing concerns about autocorrelation,\nunobserved heterogeneity, and non-stationarity. \"...Additionally, while we estimated each transmission path independently to ensure the\nHowever, several issues remain. First, because each robustness of individual coefficients within our complex interaction framework, this approach\npath is estimated separately, the analysis does not does not account for potential cross-equation dependence or correlated errors. If unobserved\naccount for potential cross-equation dependence or factors simultaneously influence multiple agricultural modernization mediators and nitrogen\ncorrelated errors, which may bias derived indirect outcomes, the derived indirect effects could be subject to minor biases. Future studies could\neffects; a brief acknowledgement of this limitation explore simultaneous equation modeling to further refine the precision of these pathways.\"\nwould improve transparency.\nSecond, the manuscript R2-3: Thank you for highlighting this critical econometric nuance. To ensure maximum transparency and\nshould more clearly explain methodological rigor, we have added clarification in the Results and Discussion, and Methods.\nthe interpretational role of\nthe lagged dependent\nChange 1: We have expanded the explanation in the \"Average transmission pathways of urbanization\" subsection of\nvariables, which materially\nthe Results and Discussion (Line 147):\naffect both the magnitude\nand meaning of the \"The direct association between urbanization and nitrogen yield and use efficiency evolves with economic development.\nestimated coefficients. All reported coefficients represent short-run marginal effects, capturing changes within a single year. Consistently high\nself-lagged values (0.6 to 0.9, Supplementary Tables S2 to S3) reflect the persistence of farming habits and soil nutrient\ncarry-over. Cross-lagged terms further account for feedbacks where past production outcomes influence current input\ndecisions. By absorbing these historical effects, our model isolates the immediate net impact of urbanization.\"\nChange 2: We have expanded the explanation in the \"Dynamic panel interaction model\" subsection of the Methods\n(Line 580):\n\"… the inclusion of this autoregressive term directly dictates the interpretational role of our estimates. Because the\nlagged dependent variable absorbs the historical accumulation of effects, the estimated coefficients for urbanization\nand modernization must be interpreted strictly as short-run (contemporaneous) marginal effects. They quantify the\nimmediate shift in nitrogen outcomes conditional on the system's state in the previous year, rather than the cumulative\nlong-run equilibrium effects, thereby yielding a more conservative and precise measure of contemporary dynamics...\"\nThird, although a robustness R2-4: Thank you for pointing out this critical omission. We have integrated a summary of the robustness check\ncheck excluding countries\nimplications directly into the paragraph that introduces the nonlinear curves (in the \"Nonlinear effects of\nwith major structural breaks\nurbanization\" section). This addition reassures the reader upfront that the subsequent U-shaped curves are driven by\nhas been added, its\ngenuine structural mechanisms rather than skewed by outlier countries experiencing state dissolutions or geopolitical\nimplications are not\nshocks.\nadequately summarized in\nthe main text; this is\nimportant because several Change: We have added the following sentences to the end of the introductory paragraph of the \"Nonlinear effects of\nkey nonlinear patterns, such urbanization\" section (Line 215):\nas the U-shaped yield\nresponse, are sensitive to \"These nonlinear patterns remain stable even after excluding countries with major structural breaks, confirming that\nsample composition.\ntrajectories are driven by fundamental economic mechanisms rather than anomalous geopolitical shocks\n(Supplementary Note S3).\"\n2. The revised manuscript offers a clearer and R2-5: We sincerely thank the reviewer for this highly perceptive conceptual point.\nbetter-justified treatment of the farm-size proxy:\nthe authors replace the previously inconsistent\nTo be fully transparent about this conceptual distinction, we have explicitly acknowledged this\nuse of per-capita cropland area with cropland\nlimitation in the revised manuscript. We noted that relying on this proxy means our results capture\narea per agricultural laborer, provide a precise\nbroad structural land-to-labor shifts rather than precise operational consolidation dynamics.\ndefinition, and move the empirical validation\ninto both the main text and the supplement,\nChange: We have substantially expanded the second point in the \"Limitations and future research\"\nwhere they report a strong correlation with\nsection (Line 413):\nactual farm size. Terminology is now\nconsistently applied, and the rationale for\n“...Second, due to data limitations across the long time series, we employed proxy variables for\npreferring this proxy over census-based\nagricultural modernization, which introduces potential interpretational biases. For instance, our\nmeasures is more clearly articulated.\nproxy for 'farm size' (cropland area per agricultural laborer) fundamentally captures labor\nNonetheless, several conceptual concerns\navailability and land-to-labor endowments rather than the true operational scale or distribution\nremain. Although correlated with observed farm\nof individual holdings. Consequently, it may not fully capture the nuanced dynamics of land\nsize, cropland per laborer fundamentally\nconsolidation or tenure fragmentation...\"\ncaptures labor availability rather than the true\noperational scale of holdings, and this\ndistinction should be explicitly acknowledged as\na limitation.\nIn addition, the theoretical R2-6: Thank you for highlighting this critical and insightful nuance. We agree that while our global macro-level\nframework implicitly\nframework captures the average structural trends of co-evolution, the actual pathways of agricultural modernization can\nassumes that farm size,\nbe highly context-specific.\nmechanization, and\nirrigation evolve jointly,\nTo address this limitation and enrich our manuscript, we have added a nuanced discussion in the Policy\nyet in many regions with\nImplications section. We explicitly acknowledge this divergence and suggest that for regions like South Asia, policy\nfragmented land tenure,\nframeworks should prioritize service-oriented mechanization models rather than relying solely on land consolidation.\nsuch as parts of South Asia,\nmechanization may\nChange: We have added the following nuanced discussion to the Policy implications section (Line 359):\nincrease independently of \"Modernization components can also diverge. In fragmented systems like South Asia, mechanization often decouples\nconsolidation; the\nfrom consolidation via machinery rental and custom hiring services (Rahman et al., 2021). Policymakers should support\nmanuscript would benefit\nthese service-oriented models through cooperatives and credit, rather than relying exclusively on land consolidation (Ma\nfrom a more nuanced\net al., 2023; Balana et al.,2022; Jiang et al., 2024.)\"\ndiscussion of such\ndivergences.\n1. Rahman, M.M., et al. Farm mechanization in Bangladesh: A review of the status, roles, policy, and potentials. J.\nAgric. Food Res. 6, 100225 (2021).\n2.. Ma, W., Marini, M.A. & Rahut, D.B. Farmers' organizations and sustainable development: An introduction. Ann.\nPublic Coop. Econ. 94, 683-700 (2023).\n3. Balana, B.B. & Oyeyemi, M.A. Agricultural credit constraints in smallholder farming in developing countries:\nEvidence from Nigeria. World Development Sustainability 1, 100012 (2022).\n4. Jiang, M., Li, J. & Mi, Y. Farmers' cooperatives and smallholder farmers' access to credit: Evidence from China. J.\nAsian Econ. 92, 101746 (2024).\n3. While the revised R2-7: Thank you for this highly insightful and constructive comment.\nintroduction is more\ncoherent and focused, the\nFollowing your suggestion, we have condensed the descriptive background in the Introduction (Line 33) to bring the\noverall novelty of the study\ncore research gap into focus more quickly:\nstill requires clearer\narticulation. Prior global\nChange 1:\nresearch has examined\nurbanization, farm \"Urbanization, one of the most significant global trends of the 21st century, profoundly reshapes land-use patterns and\nstructure, and nitrogen agricultural system1-3. A critical nexus in this transformation is the agricultural nitrogen (N) cycle. While N-driven\ntrends largely in isolation, agricultural intensification has been instrumental in global food security4,5, its inefficient management has triggered\nand the distinctive severe, cascading environmental challenges6-8. Consequently, understanding the complex drivers shaping nitrogen\ncontribution of this work\nmanagement is paramount for navigating the trade-offs between productivity and sustainability in an increasingly\nlies in integrating these\nurbanized world.\"\ndomains within a\nnonlinear, moderating\nFurthermore, we restructured the study aims to highlight our central conceptual innovation. By contrasting our\nframework, yet this\nframework with linear global models that mask threshold effects and regional studies that lack generalizability across\nmessage is at times\ndevelopment stages, we have revised the Introduction to emphasize our income-specific mechanistic hypotheses (Line\nobscured by extended\n58):\ndescriptive background on\nglobal trends. The authors\nChange 2:\nshould sharpen the\n\"To fill this gap, the central conceptual innovation of this study lies in integrating urbanization, agricultural\nnarrative by explicitly\nmodernization, and nitrogen trends into a unified, nonlinear moderating framework. Prior research has largely\nemphasizing why\nexamined these domains in isolation. While linear global models are insightful, they often focus on estimating a net\ncapturing joint\naverage effect, thereby masking critical threshold effects and structural phase transitions. Conversely, regional case\nnonlinearities yields\nstudies offer rich granular details but often lack the macro-level perspective necessary to generalize how these\ninsights that neither linear\nrelationships evolve across distinct stages of economic development. By capturing these joint nonlinearities, our\nglobal models nor regional\napproach yields macro-systemic insights that neither linear models nor regional studies can provide alone.\ncase studies can provide,\nand by positioning this\nWe hypothesize that the association of urbanization with nitrogen fertilizer application, nitrogen yield, and nitrogen\nintegrative approach as the\nuse efficiency (NUE) is not uniform, but is profoundly nonlinear and critically contingent upon the stage of economic\ncentral conceptual\ndevelopment and the prevailing levels of farm size, machinery intensity, and irrigation intensity. In low-income\ninnovation of the\ncountries, the primary challenge is whether farm size can grow fast enough to compensate for labor loss and create a\nmanuscript.\nviable scale for machinery intensity15, 16. In contrast, in high-income countries, the challenges shift to managing the\ntrade-offs of a highly optimized system, questioning whether further advances in machinery intensity and irrigation\nintensity can provide resilience against diminishing returns17. To test these hypotheses, we compiled a comprehensive\npanel dataset spanning 139 countries from 1961 to 2021...\"\nReferee #3\nComment Response and Change\nI co-reviewed this manuscript with one of the reviewers who provided R3-1: Thank you for letting us know about the co-review process. We\nthe listed reports. This is part of the Nature Communications initiative appreciate the thorough attention that your collaborative review has brought\nto facilitate training in peer review and to provide appropriate to our manuscript.\nrecognition for Early Career Researchers who co-review manuscripts.\nReferee #4\nComment Response and Change\nI appreciate the authors' substantial revisions in response to my previous R4-1: Thank you for the positive feedback. We have addressed each of your\ncomments. These changes have significantly enhanced the clarity and remaining minor comments with great care, particularly in strengthening the\nimpact of the manuscript. I have only a few minor comments remaining:\ndialogue between our empirical results and the existing literature, which\nfurther consolidates the theoretical contributions of this study.\n1. Regarding my previous R4-2: We appreciate your constructive suggestion. Following your excellent advice, we have structured our literature\ncomment 1, I note clear dialogue across the three income groups. We carefully aligned these dialogues with the core conclusions\nimprovements in the revised\nhighlighted in our Abstract, referencing high-quality, recent empirical studies to demonstrate how our global\nsummary of existing literature\nframework validates regional challenges in developing contexts while revealing structural shifts in advanced\nand the identified knowledge\neconomies.\ngap. To further strengthen the\npaper, I suggest that the authors\nconsider referencing relevant Specifically, we have woven these dialogues into the following three sections:\nfindings or discussions in the\nResults and Discussion section. 1. Low- & Lower-Middle-Income Countries: Urbanization suppresses nitrogen yield, with larger farm size\nIt would be beneficial to\nbuffering early losses\nhighlight how this study\nTo explain the initial yield suppression, we linked our findings to regional studies in Sub-Saharan Africa,\nresponds to, or extends,\ndemonstrating how scale expansion acts as a necessary buffer against capital constraints.\nprevious work or discussions on\nthe relationship between\nChange 1: We have integrated the following literature dialogue into the discussion of Low- and Lower-middle-\nurbanization and nitrogen use.\nincome countries (Line 224):\n\"...For nitrogen yield, urbanization exerts a prolonged suppressive effect under low modernization. This aligns with\nrecent evidence of widespread yield stagnation among African smallholders (Wollburg et al., 2024; Lobell &\nLee, 2025). This stagnation is driven by rural outmigration and capital constraints that hinder mechanization (Patil et\nal., 2026). Our framework reveals that larger farm size acts as a critical structural buffer against these early\nlosses. Scaling up helps overcome capital barriers, as productivity rises significantly only after farms reach a\nviable size for mechanization (Muyanga & Jayne, 2019). Consistent with this technological threshold, the impact\nof urbanization on nitrogen yield reverses dramatically under high modernization: the turning point of the U-\nshaped curve moves from 52% to 34% urbanization...\"\n2. Upper-Middle-Income Countries: Modernization reconfigures associations, introducing trade-offs between\nenvironmental gains and yield growth\nTo contextualize the observed input intensification, we connected our results to the \"input substitution\" phenomenon\nin rapidly urbanizing Asia, highlighting that consolidation and mechanization must act synergistically to optimize\nefficiency and reduce fertilizer overuse.\nChange 2: We have added the following discussion to the Upper-middle-income context (Line 299) to provide\ndeeper context for the environmental benefits:\n\"...Machinery intensity and irrigation intensity enhance the nitrogen fertilizer application-reducing effect of\nurbanization, producing clear environmental benefits. This reflects the input substitution phase where urban labor\npull drives severe fertilizer overuse to compensate for lacking fixed inputs (Steinhuebel-Rasheed et al., 2025;\nWu et al., 2018). However, our results demonstrate that this trade-off can be optimized. Consistent with recent\nempirical studies, land consolidation must be coupled with mechanization to achieve genuine efficiency gains\nand positive regional spillovers (Abdulwahid et al., 2024; Sheng et al., 2019). This synergy provides the necessary\nscale to significantly reduce fertilizer intensity (Zhen et al., 2024; Li et al., 2023). However, they also intensify the\nnegative association of urbanization with nitrogen yield...\"\n3. High-Income Countries: Modernization mitigates adverse impacts, while further farm size expansion leads\nto diminishing returns\nTo explain why farm size becomes a negative moderator, we used our global model to extend the literature on the\n\"limits of scale,\" referencing recent evidence on yield plateaus and ecological trade-offs in developed nations.\nChange 3: We have updated and streamlined the discussion on High-income countries (Line 318):\n\"In high-income countries, modernization dynamics fundamentally shift. Farm size expansion, once a driver of\nefficiency, now acts as a negative moderator that exacerbates the decline in both nitrogen yield and NUE (Fig.\n4). This reversal indicates diminishing or even negative returns to scale. It aligns with recent evidence that highly\noptimized and large-scale operations face yield stagnation and limited environmental benefits (Gerber et al., 2024;\nLarsen & Noack, 2021; Sirami et al., 2019). In contrast to further scale expansion, technological modernization\nthrough machinery intensity and irrigation intensity emerges as the principal source of resilience. \"\n1. Abdulwahid, N., Bakari, L., Hussein, A., Kawa, S. K., Lavoe, F., Mwisomba, T., ... & Wineman, A. (2024).\nSpillover effects of medium-and large-scale farms on smallholder farmers in Tanzania: Evidence from the\nNational Sample Census of Agriculture 2019/20. World Development Perspectives, 34, 100590.\n2. Gerber, J. S., Ray, D. K., Makowski, D., Butler, E. E., Mueller, N. D., West, P. C., ... & Sloat, L. (2024).\nGlobal spatially explicit yield gap time trends reveal regions at risk of future crop yield stagnation. Nature\nFood, 5(2), 125-135.\n3. Larsen, A. E., & Noack, F. (2021). Impact of farm size on biodiversity: Modeling the effect of farm size and\nfarm specialization after controlling for environmental variables. Journal of Applied Ecology, 58(12), 2806-\n2816.\n4. Li, Y., Huan, M., Jiao, X., Chi, L., & Ma, J. (2023). The impact of labor migration on chemical fertilizer use of\nwheat smallholders in China- mediation analysis of socialized service. Journal of Cleaner Production, 394,\n136366.\n5. Lobell, D. B., & Lee, R. J. (2025). Crop productivity in southern Africa is stagnant despite moderate climate\ntrends. Nature Food, 6, 762-765.\n6. Muyanga, M., & Jayne, T. S. (2019). Revisiting the farm size-productivity relationship based on a wide range\nof farm sizes: Evidence from Kenya. American Journal of Agricultural Economics, 101(4), 1140-1163.\n7. Patil, P., Schuenemann, F., & Strzepek, K. (2026). Labour Dynamics, Factor Substitution and Capital Adoption\nin Sub-Saharan African Agriculture. Journal of Agricultural Economics, 77(1), 183-204.\n8. Sheng, Y., Ding, J., & Huang, J. (2019). The relationship between farm size and productivity in agriculture:\nEvidence from maize production in northern China. American Journal of Agricultural Economics, 101(3), 790-\n806.\n9. Sirami, C., Gross, N., Baillod, A. B., Bertrand, C., Carrié, R., Hass, A., ... & Fahrig, L. (2019). Increasing crop\nheterogeneity enhances multitrophic diversity across agricultural regions. Proceedings of the National\nAcademy of Sciences, 116(33), 16442-16447.\n10. Steinhuebel-Rasheed, L., Christiaensen, L., Minten, B., Swinnen, J., & Vandercasteelen, J. (2025). Cities and\nAgricultural Development. Annual Review of Resource Economics, 17(1), 401-421.\n11. Wollburg, P., Bentze, T., Lu, Y., Udry, C., & Gollin, D. (2024). Crop yields fail to rise in smallholder farming\nsystems in sub-Saharan Africa. Proceedings of the National Academy of Sciences, 121(21), e2312519121.\n12. Wu, Y., Xi, X., Tang, X., Luo, D., Gu, B., Lam, S. K., ... & Chen, D. (2018). Policy distortions, farm size, and\nthe overuse of agricultural chemicals in China. Proceedings of the National Academy of Sciences, 115(27),\n7010-7015.\n13. Zhen, S., Wang, W., Peng, Y., & Liu, Y. (2024). Does cross-policy mix work on agriculture total factor\nproductivity: analyzing by the impact of agricultural mechanization and large-scale operations in China.\nEnvironment, Development and Sustainability, 1-2\n2. Regarding the selection of R4-3: We greatly appreciate your constructive feedback. We have expanded the Limitations section to explicitly\nagricultural modernization address the uncaptured aspects of each proxy and their impact on result interpretation.\nindicators, I appreciate both the\nsubstantial revisions and the\nrationale provided. I\nrecommend including a brief\ndiscussion, perhaps in the\nsecond point of the limitation\nsection, on the potential biases\nassociated with these proxy\nvariables and how they may\naffect the interpretation of farm\nsize expansion, mechanization,\nand irrigation. For example,\nwhat relevant aspects might not\nbe fully captured, and how\ncould this influence the results?\nAdditionally, please note the\nrepetition of \"Machinery\nintensity\" in Line 456.\n3. I encourage the authors to R4-4: We appreciate the reviewer’s helpful suggestion. Following your constructive advice, we have chosen to retain\nreconsider the title of the the title to contrast with the subsequent \"Nonlinear effects\" section, but we have added a clear and explicit\nsecond Results section,\nexplanation at the very beginning of this section.\n\"Average transmission\npathways of urbanization,\"\nas\nits meaning is not entirely clear.\nIf the authors prefer to retain\nthis title, I suggest providing a\nbrief explanation at the\nbeginning of the section to\nclarify what is meant by\n\"Average transmission\npathways.\"\nReferee #5\nComment Response and Change\nI co-reviewed this manuscript with one of the reviewers who provided R5-1: Thank you for letting us know about the co-review process. We\nthe listed reports. This is part of the Nature Communications initiative appreciate the thorough attention that your collaborative review has brought\nto facilitate training in peer review and to provide appropriate to our manuscript.\nrecognition for Early Career Researchers who co-review manuscripts."}]

Decision letter