The field experiment was initiated in 2013 at the Agriculture and Agri-Food Canada (AAFC)’s Beaverlodge Research Farm (55°12’ N 119°24’ W), located in northwestern Alberta, Canada. The site features a slope gradient of 2.5% and average elevation of 745 m. The predominant soil at the site was a dark gray Luvisol developed from lacustro-till and glacial-lacustrine deposits (Broersma et al., 1997). Soil analysis of initial soil samples taken in 2013 at depths of 0-15, 15-30, and 30–60 cm depths showed organic matter content of 5.5, 4.3, and 3.6%; pH values of 5.43, 6.13, and 6.43; and bulk densities of 1.3, 1.5, and 1.6 g cm⁻³, respectively (Khanal et al., 2021). Soils at the surface layer (0–15 cm) were loamy in texture (31.3% sand, 46.0% silt, and 22.7% clay), and the sub-surface layer (15–30 cm) had a clay loam texture (28.1% sand, 36.1% silt, and 35.8% clay). The experimental region has a continental climate with a mean annual precipitation of 435 mm based on the 30-year average of weather data (1991–2020) (Anonymous, 2024). Typically, 50–60% of the annual precipitation occurs during the primary growing season (May to September). However, a seasonal shift in precipitation distribution during transitional months like spring and fall, higher frequency of extreme rainfall events, and longer intervals between rain events are more noticeable, leading to crop failure in recent years (Yang et al., 2019). The mean monthly maximum temperatures varied from 17 to 22 °C between May and September (Anonymous, 2024). The monthly distribution of precipitation and temperature during the growing seasons from 2013 to 2024 are presented in Figure 1.

The experiment was conducted in a split-plot design with four replicates. The main-plots were assigned to eight different cropping sequences, while the sub-plots nested within the main-plots received three different nitrogen (N) levels (0, 45, and 90 kg N ha^-1), randomly assigned during the establishment year. All sub-plots were 2.5 m in width and 8.0 m in length (20 m²), and a 0.6 m buffer area was maintained between main plots. Nine species including the seed crops of two perennial forage legumes (alsike clover [Trifolium hybridum L.], and red clover [Trifolium pratense L.]), three perennial grasses (creeping red fescue [Festuca rubra L.], meadow bromegrass [Bromus riparius Rehm.], and timothy [Phleum pratense L.]), and four annual crops (canola [Brassica napus L.], wheat [Triticum aestivum L.], pea [Pisum sativum L.], and barley [Hordeum vulgare L.]) constituted different phases of eight different cropping sequence treatments. More detailed information concerning years and crops grown in each cropping sequence can be found in Table 1 and Supplementary Figure 1. Six of the eight cropping sequence treatments were diversified with perennial forage legumes and grasses to compare with annual-based traditional sequences (i.e. canola-wheat only or in combination with barley and/or pea). During the seeding of annual crops, N fertilizer treatments were applied as band placement, while for perennial grasses, it was broadcast in the fall over the production years. All perennial forage seed crops underwent an initial establishment year characterized by vegetative growth, followed by one year of economic seed production for legumes and two years for grasses.

The experiment was conducted under a no-tillage management system since its establishment in 2013. The test crops were seeded each year with a 2.4 m wide air seeder (Valmar Cross Slot IP Limited) with 30-cm row spacing (8 rows). Basal fertilizers for phosphorus (13.1 kg P ha^-1), potassium (16.6 kg K ha^-1), and sulfur (20 kg S ha^-1) were band-placed at sowing time. The broad-spectrum herbicides glyphosate (Roundup WeatherMAX, 3.2 L ha⁻¹) was applied twice, once in the fall and once in the spring, to effectively terminate perennial forage crops after their production years and to establish subsequent crops (Supplementary Table 1). Pests were monitored throughout the growing season; however, no infestations warranted the application of pesticides. The perennial crop residues were mowed post-harvest and uniformly distributed across the plot to maintain homogeneity and facilitate reproductive tiller production.

The crops were sampled at physiological maturity from the middle two rows by taking 4 m row lengths from each treatment sub-plot for seed and aboveground biomass yield, and subsequent post-harvest analyses. The plant samples were air-dried to a constant weight to measure total aboveground biomass yields on an area basis. The samples were then threshed mechanically using an ALMACO LPR thresher (ALMACO Nevada, USA) to measure grain and straw yields.

A random sub-sample of cleaned seed and threshed straw from each sub-plot was ground using an electric grinder (Marathon Electric, Wausau, Wisconsin 54401, USA). Nitrogen concentration in the plant samples was determined using the dry combustion method described by Nelson and Sommers (1982). Briefly, 10–20 mg of ground plant materials were encapsulated in tin foils with 6 × 6 × 12 mm dimensions. The processed samples were then analyzed by a Vario EL Cube (Elementar Analysensysteme GmbH, Germany), where the combustion and reduction temperatures were maintained at 1150 °C and 850 °C, respectively. The tissue N concentration data were used to calculate N uptake and apparent fertilizer N recovery under experimental management conditions (see below).

For uniform comparison among the sequences, the annual and perennial forage crop seed yield was expressed as Canola Equivalent Yield (CEY) in kg ha^-1. The CEY was calculated using the following formula (Liu et al., 2019; Khanal et al., 2021).

The prices of both canola and non-canola seeds were obtained from Agriculture Financial Services Corporation (Agriculture Financial Services Corporation, 2024). The seed commodity prices varied every month, and the yearly average price was used in the calculation (Supplementary Table 1). The cumulative CEY for each treatment subplot was calculated by summing the CEY of individual crops grown under the cropping system.

Gross revenues were calculated based on the seed yields of all crops grown in each treatment sequence from 2013 to 2024. For each crop year, seed yields were multiplied by the annual average sale prices of seeds of the respective crops, which were summed over the years to obtain the total cumulative gross revenues. Gross margins were calculated as the difference between total revenue and variable costs across treatment sub-plots, including seed and seeding, fertilizer and fertilizer application, herbicide and herbicide application, mowing, harvesting, and storage. The variable inputs and operational costs were valued in Canadian Dollars ($CDN) and calculated for each treatment by adding all expenses on a per-hectare basis using the current price ($ ha⁻¹) (Supplementary Table 2). Other costs, such as equipment repairs, depreciation, land rent, common operating costs, insurance and property taxes, were not considered for the analyses as they did not vary across cropping sequences.

This study focused exclusively on seed yield as the measure of crop productivity and did not include forage biomass as an economic output. While this may appear to limit a complete assessment of overall crop output, the decision aligns with the research objective centered on seed production. The study location, a region recognized as a major exporter of forage seed, has limited demand for hay or biomass use, making seed yield the most relevant and economically significant metric. The CEYs, gross revenues, and gross margins were calculated based on the average annual prices of variable inputs and outputs; therefore, they do not account for how the profitability of crop sequences may vary under price elasticity.

For uniform comparison among the cropping sequences, NUE was calculated in terms of aboveground biomass yield and agronomic NUE for CEY. Accordingly, the aboveground biomass and the agronomic NUE for CEY were calculated using the following formula described by Khanal et al. (2021), which represents the amount of additional aboveground biomass and CEY seed harvested per kilogram of N fertilizer applied.

Where CEY is Canola Equivalent Yield, Nk is the applied N rate (45 or 90 kg ha^-1, and N0 is without N.

The N-uptake (kg ha⁻¹) by the crops grown in 2023 and 2024 was calculated by multiplying the dried sample biomass (straw and grain) by the corresponding N concentration. To examine apparent fertilizer N recovery, N uptake in non-nitrogen fertilized plot was subtracted from the N uptake in a plot at given nitrogen fertilizer rate and then divided by the corresponding fertilizer rate (Nk) (Craswell and Godwin, 1984), as shown in the equations below:

where tissue N concentration and aboveground biomass were expressed in percentage and kg ha^-1, respectively.

where Nup k is N uptake in a plot at given N fertilizer rate, Nup control is N uptake in non-N fertilized plot and Nk is the applied N rate (45 or 90 kg ha^-1).

All data analyses were performed using the SAS OnDemand for Academics (SAS Institute, Inc., Cary, NC) in compliance with the standard procedure of split-plot design. The assumptions of normality (PROC UNIVARIATE) and homogeneity of variances (Levene’s test) were examined prior to analysis. Data transformations were not required to meet normality assumptions. The data were then analyzed using the PROC GLIMMIX procedure, and the Kenward-Roger approximation method was used to determine degrees of freedom for the treatment factors. Cropping sequences and N rates were treated as fixed effects, whereas blocks and their interactions with cropping sequences and N were considered random effects.

Mathematically, the analytical model is expressed as:

Indices:

Blocks: i = 1,…., r

Cropping sequences (whole-plot levels): j = 1,…., a

Nitrogen: k = 1,…., n

Effects:

Fixed: μ (overall mean), S_j (cropping sequence), N_k (nitrogen), (SN)_jk (interaction)

Random: B_i ~N (0, σ^2B) (Block), BS_ij ~ N(0, σ^2BC), (whole plot error: block x cropping sequence), ϵ_ijk ~ N(0, σ²) (sub-plot residual)

Error strata and tests:

Cropping sequence (c): test S_j against the whole-plot error (BS_ij).

Nitrogen (N): test N_k against the sub-plot residual ϵ_ijk,

Interaction (S x N): test (SN)_jk against the sub-plot residual ϵ_ijk.

Wherever treatment differences were found, Tukey’s Honestly Significant Difference (HSD) procedure was used to separate means at a 5% significance level.

The cumulative CEY varied significantly for the interaction between cropping sequences and N treatments (P <.0001) over the production years 2013-2024 (Figure 2). The legume to vernalizing grass rotation (S₄), with intermittent reoccurrence of red clover leading to meadow bromegrass phases, resulted in higher CEY than other cropping sequences across three N levels. In this cropping sequence, red clover was grown biennially in two alternating phases for four years, and meadow bromegrass was grown for three consecutive years over the twelve-year production period. The highest CEY in this sequence was driven by the higher yield and premium price of meadow bromegrass seed in the production years of 2023 and 2024. The seed yield of meadow bromegrass was 2420 and 910 kg ha^-1 in 2023 and 2024, respectively, which was significantly higher than the average productivity of meadow bromegrass in the Peace region of western Canada, as reported by Azooz and Johnson (2014). Moreover, the seed price of meadow bromegrass was nearly eight times higher than that of canola ($4.82 vs. $0.74 per kg in 2023 and $4.79 vs. $0.59 per kg in 2024) in those years (Agriculture Financial Services Corporation, 2023, 2024).

Notably, the seed price of perennial forages like meadow bromegrass increased substantially compared to canola in recent years. For example, the seed price of meadow bromegrass rose by 69% from $2.83 per kg in 2018 to $4.79 per kg in 2024. In contrast, the seed price of canola showed a modest increase of 28% from $0.46 per kg in 2018 to $0.59 per kg in 2024 (Agriculture Financial Services Corporation, 2018, 2024). The substantial price increase indicates an opportunity to leverage demand in emerging markets, which is a key consideration for crop planning (Edgerton, 2009). In contrast, lower CEY was recorded in the alsike clover (a perennial legume) to timothy (a non-vernalizing grass) rotation (S₅), irrespective of N fertility levels. The lower seed productivity of alsike clover (80 kg ha^-1 in 2018) resulted in lower CEY in S₅. Our visual observations indicate that excessive moisture during alsike clover’s flowering to maturity stages critically affected seed productivity in 2018. The physiological legacy of drier conditions in the preceding fall and winter could be the reason behind the shorter plants and smaller heads of timothy during the production periods of 2023 and 2024. Timothy displays reproductive asynchrony due to weather variability, resulting in reduced seed yield. Earlier in 2021, the Peace River region experienced a heatwave that significantly affected crop productivity, including timothy (Pandey, 2023). Garwood and Sinclair (1979) reported that timothy is sensitive to soil moisture deficit at the deeper soil horizons. The CEY of fescue-dominated sequence (S₇) was higher than that of the timothy- and clover-integrated cropping sequences, while remaining comparable to CEYs of annual crop-based sequences (S₁ and S₂). As fescue species are recognized for their eco-physiological benefits, such as improved soil structure, enhanced nutrient cycling, increased soil organic matter, and augmented water infiltration and retention (Hartwig and Ammon, 2002; Athar et al., 2022; Skersiene et al., 2023), their inclusion in crop diversification could be a sustainable approach for the study region.

The gross revenue was significantly influenced by cropping sequences (P <.0001), N rates (P = 0.005), and their interactions (P < 0.02) (Figure 3). The highest gross revenue was recorded from the legume to vernalizing grass rotation (S₄), with intermittent recurrence of red clover leading to meadow bromegrass phases, compared to annual-based and other diversified cropping sequences. As explained in the preceding sub-section, the economic advantage of this rotation was primarily influenced by crop prices and yields of different crops.

The higher seed price contributed to the significant increase in the gross revenue of legume-to-vernalizing grass rotation cropping sequence. Conversely, the cropping sequence with the lowest gross revenue was recorded from the legume-to-non-vernalizing grass rotation (S₅), with intermittent recurrence of alsike clover leading to timothy phases, across all N fertility levels (Figure 3). This was primarily due to the lower productivity of both timothy and alsike clover during their production phases. The average seed yield of timothy was less than 200 kg ha^-1, while that of alsike clover was 80 kg ha^-1, both substantially lower than previously reported in the region (Yoder, 2004; John and Ogle, 2008).

Understanding crop adaptative behavior and economic dynamics is essential for developing resilient crop rotational strategies amid climate variability and shifting market conditions. Bansal et al. (2024) conducted a comparative analysis of diversified and conventional crop rotation systems in South Dakota, USA, and found that higher economic returns in diversified cropping systems are associated with increased crop yields. The crop-specific sensitivity to weather anomalies could have resulted in seasonal yield variability of crops in our study. The productivity gains typically increase expected revenue per unit area as they directly translate to higher output and potential income, but higher supply can also lower crop prices (Pouliot, 2019). Therefore, the commodity prices can be a key driver of on-farm profitability (Canales et al., 2021). The gross revenues of different forage-based sequences (S₆, S₇, and S₈) were comparable to those of annual cropping systems (S₁ and S₂), irrespective of different N application rates. This result suggests that well-managed perennial-integrated cropping systems can be viable alternatives to recurrent annual crop rotations enhancing the resilience of farm income.

Gross margin is considered here as a primary indicator of the economic viability of diversified cropping systems, as it accounts for the differential costs across different cropping sequences. Overall, the highest gross margins were recorded from the legume (red clover) to vernalizing grass (meadow bromegrass) rotation (S₄) with an average of $21,143 ha^-1, and the lowest in the legume (alsike clover) to non-vernalizing grass (timothy) rotation (S₅) with an average of $1,688 ha^-1 (Figure 4). The legume to vernalizing grass rotation (S₄) consistently exhibited a higher gross margin across all N rates. High yields and strong market prices for meadow bromegrass in 2023 and 2024 drove the profitability of the S₄ sequence, underscoring the importance of understanding crop adaptability and market signals.

Abundant literature substantiates the economic viability of diversified cropping systems. Sánchez et al. (2022) performed a meta-analysis of 119 peer-reviewed articles and found that diversified systems generate profits comparable to or greater than simplified farming systems. Alcon et al. (2024) studied the medium- and long-term economic performance of diversified cropping systems and found higher gross margins from diversification than monocrop farming practices. Although the magnitude of profit margins depends on crops, diversification can promote better financial outcomes by mitigating risks related to extreme weather, pest outbreaks, yield loss, and price volatility (Gaudin et al., 2015; Mortensen and Smith, 2020; Khanal et al., 2021; Mihrete and Mihretu, 2025). The thoughtful selection of crops with complementary growth characteristics and resource requirements can enhance profit margins by reducing input costs (Khanal et al., 2021; Mihrete and Mihretu, 2025). Overall, strategic crop planning can enhance efficiency, reduce input costs, and improve profit margins in diversified farming systems.

Farm profitability is driven by both higher yields and cost savings from reducing or appropriately substituting inputs such as fertilizers and pesticide costs. Higher crop yields achieved through excessive use of inputs beyond optimal levels can lead to diminishing returns, where additional inputs do not proportionally increase yields but raise costs (Reader et al., 2018; Pittman, 2020; Liang, 2023). As expected, we observed a significant improvement in gross margin with the increase of N rates from 0 to 45 kg N ha^-1 for most of the diversified treatment sequences. However, a further increase to 90 kg N ha^-1 did not yield significant economic benefits. These findings highlight the potential to fine-tune N management in cropping systems that include annual or perennial legumes. Analyzing over 34 years of experimental results in the USA, Breza et al. (2023) concluded that cropping systems complexity enhanced organic nutrient recycling, but N fertilization counteracted the beneficial effects of crop diversification. Moreover, the excess N may trigger environmental issues like nitrate leaching (De Laporte et al., 2021), leading to soil degradation (Awadelkareem et al., 2023), water pollution (Akinnawo, 2023), and increased greenhouse gas emissions (Menegat et al., 2022) and impairment of beneficial soil microbiome functionality (Beltran-Garcia et al., 2021). Therefore, rational N fertilization is essential to improve productivity and reduce environmental footprint of diversified cropping systems.

During 2013–2024, the effects of cropping sequences (P <.0001), N rates (P <.0001), and their interaction (P = 0.0002) were significant for cumulative aboveground biomass yield (Figure 5). The conventional and diversified annual cropping systems (S₁ and S₂) consistently produced higher aboveground biomass yields across the N application rates compared with forage-based sequences (S₃ to S₈). Among perennial forage-based systems (S₃ to S₈), the legume-to-vernalizing grass rotation (S₄) resulted in the highest aboveground biomass yield, whereas the turf grass-to-legume rotation (S₇) had the lowest aboveground biomass yield across all N rates. High aboveground biomass yield in the cropping sequence with intermittent inclusion of red clover (legume) and meadow bromegrass (vernalizing grass) was concurrent with an exceptionally high seed yield of meadow bromegrass under favorable weather and crop-growing conditions in 2023 and 2024. In contrast, the lower aboveground biomass yield of creeping red fescue led to reduced overall aboveground biomass productivity in the turfgrass-to-legume rotation (S₇) across all N rates. It is well established that perennial grasses have high root-to-shoot ratios because of their extensive, deep root systems and contribute to soil health and other ecosystem services (Beard, 1973; Bolinder et al., 2002; Thiagarajan et al., 2018). The perennial growth habit, characterized by greater allocation to belowground biomass can enhance carbon sequestration and help mitigate climate change risks (Franzluebbers et al., 2014; Ojija, 2024; Skersiene et al., 2024).

Enhanced aboveground biomass productivity in annual cropping systems is primarily attributed to their higher responses to N fertilization (Larsen et al., 2024). This responsiveness is due to their efficient uptake and utilization of available N, leading to accelerated growth rates and higher yields (Omara et al., 2019; Liu Z. et al., 2022; Wen et al., 2023). In contrast, forage-based sequences with perennial grasses and legumes often exhibit slower initial growth and allocate a larger proportion of resources to root development, which can result in slower aboveground growth compared with annual-based sequences (Bolinder et al., 2002; Hakala et al., 2009; Thiagarajan et al., 2018; Vico and Brunsell, 2018). It is also noteworthy that aboveground biomass yield measurements of perennials were only performed during the seed production phases (one year after crop establishment), because the plant growth in both perennial forage and seed production systems is insufficient to justify harvesting during the establishment year in this region. The exclusion of perennial biomass in the establishment year partly contributed to lower aboveground biomass yields observed in perennial forage-based diversified sequences. We also observed that aboveground biomass yield increased substantially with the application of 45 kg N ha⁻¹, and no further yield improvement was realized at higher N level (90 kg N ha⁻¹). The perennial crops, therefore, contribute to sustainable agroecosystems because of lower N demand and greater belowground biomass inputs.

The mean NUE indices for the cropping sequence treatments from 2013 to 2024 are summarized in Table 2. The highest agronomic NUE was observed in the perennial legume-to-vernalizing grass rotation (S₄) at both N application rates of 45 and 90 kg N ha^-1, followed by the conventional annual rotation (S₁) at 45 kg N ha^-1. While annual crops generally exhibited higher agronomic NUE than perennials, the yield-specific N efficiency (measured as CEY-based agronomic NUE) can be higher in perennial forage-based sequences (Weih et al., 2011; Khanal et al., 2021). The higher agronomic NUE in the perennial legume-to-vernalizing grass cropping sequence, mainly due to the higher seed production of meadow bromegrass in 2023 and 2024, showed greater efficiency in converting applied N into grain yield through improved uptake and assimilation (Fageria and Baligar, 2005; Govindasamy et al., 2023). In contrast, cropping sequences S₅, S₆, and S₈ exhibited the lowest agronomic NUE across both N application rates. The magnitude of variation in agronomic NUE among cropping sequences is largely influenced by the seed yields of different crops, which vary considerably among crop species. Numerous studies have substantiated that agronomic N use efficiency was influenced by crop types and varieties, growth characteristics, growing-season weather, and agronomic management practices (Ivić et al., 2021; Khanal et al., 2021; Govindasamy et al., 2023).

Similar to agronomic NUE, aboveground biomass NUE varied greatly across cropping sequence treatments (Table 2). The highest aboveground biomass NUE was observed in conventional annual rotation (S₁) at both N levels, indicating a strong ability to utilize soil-applied N for aboveground biomass and grain production. Annual crops often exhibit higher aboveground biomass NUE than perennials because of their rapid growth rates, shorter life cycles, and intensive nutrient uptake patterns (Weih et al., 2018). Annual crops usually allocate more resources to aboveground biomass and seed production, whereas perennial grasses tend to allocate more resources to root development (Tilman et al., 2006; Means et al., 2022). These differences explain why annual crops display higher aboveground biomass NUE than perennials. Previous studies suggest that annual crops depend heavily on external fertilizer inputs, whereas perennials rely on vigorous root systems and microbial symbiosis for nutrient retention and recycling (Mueller et al., 2012).

The N concentration and uptake by straw and grain were influenced by treatment interactions, as shown in Tables 3 and 4. As anticipated, higher N concentrations were recorded in the plant tissues of legume species compared with annual cereals and perennial grasses in both 2023 and 2024. Previous research reported that legumes have evolved mechanisms to maintain a high N concentration in both shoots and seeds (McKey, 1994; Schenk et al., 1995). Naturally, the N-fixing symbiosis with rhizobia contributes to the accumulation of more N and the maintenance of elevated tissue N concentrations in legume plants (Zahran, 1999). However, despite this higher N concentration, legumes often exhibit lower rates of photosynthesis and plant growth per unit of N in the plant relative to cereals (Evans, 1989; Del Pozo et al., 2000). Conversely, cereals typically have a more efficient use of N in photosynthesis, as they invest more N into photosynthetic enzymes like Ribulose-1,5-bisphosphate carboxylase/oxygenase (RUBISCO) and have higher rates of C assimilation per unit of leaf tissue N (Makino, 2011).

While tissue N concentration increased with higher N rates, N uptake varied significantly among cropping sequence treatments because of differences in crop species. The observed variations in straw and grain N uptake were attributed to differences in aboveground biomass production among crops grown under the tested treatment sequences (Figure 5, Tables 3, 4). For example, thegrain N concentration of alsike clover was highest in 2024, but the uptake was among the lowest due to low seed yield. Conversely, the highest grain N uptake was recorded in wheat during 2023 and 2024, which was attributed to its higher seed yield compared to other test crops. With an external N supply, annual crops often achieve greater biomass accumulation and higher yields compared to perennials under similar environmental conditions (Fageria and Baligar, 2005; Govindasamy et al., 2023). The results suggest that annual grain crops tend to have greater total N uptake, likely due to their ability to utilize available N efficiently (Guinet et al., 2020). A similar result was observed in ANR, as annual crops recorded the highest recovery and demonstrated superior fertilizer N uptake. Unlike cereals, perennial legumes are less responsive to applied N fertilizers since they meet their N requirements mostly through symbiotic N fixation (Šidlauskaitė and Kadžiulienė, 2023; Tang et al., 2024). Furthermore, the ANR calculation considered only the fertilizer-derived N uptake, excluding any contribution of atmospheric N₂ fixation, which was not estimated in this study. The capacity of legumes for N-fixation can make their uptake of fertilizer N less responsive, resulting in a lower or negative N recovery. Similarly, legumes can leave more N in the soil than is removed by uptake, resulting in negative ANR values (Guinet et al., 2020). Excessive N application to legumes can even suppress root nodulation and N fixation (Sharma and Bali, 2017), suggesting a complex interaction of ecophysiological feedback mechanisms.