Biofertilizers can enhance nitrogen use efficiency of sugarcane.

Fertilizers are costly inputs into crop systems. To compensate for inefficiencies and losses from soil, farmers apply on average double the amount of nitrogen (N) fertilizer acquired by crops. We explored if N efficiency improves with biofertilizers formulated with organic waste, mineral N or plant growth-promoting rhizobacteria (PGPR). We compared treatments receiving mineral N fertilizer or biofertilizers at industry-recommended (100%) or lower (60%) N rates at two commercial sugarcane farms. Biofertilizer at the 60% N-rate generated promising results at one farm with significantly higher biomass and sugar yield than the no-N control, which matched the 100% mineral N treatment. This yield difference was accompanied by a shift in microbial diversity and composition. Correlation analysis confirmed that shifts in microbial communities were strongly linked to soil mineral N levels, as well as crop productivity and yield. Microbial co-occurrence networks further revealed that biofertilizer, including treatments with an added PGPR, can enhance bacterial associations, especially in the context of complex fungal networks. Collectively, the results confirm that biofertilizers have quantifiable effects on soil microbial communities in a crop system setting, which underscores the opportunities for biofertilizers to promote N use efficiency and the circular N economy.

Introduction

Maximizing crop yields with minimum input and environmental harm is a key goal for agriculture. However, conventional mineral fertilizers, including N fertilizers, are generally applied at higher rates than strictly needed by crops due to the inherent inefficiencies and risk of nutrient loss from soil (Udvardi et al., 2021). This overapplication has led to excess nutrient loads in the biosphere, wide‐ranging pollution and soil deterioration, including reduced microbial biodiversity (Umrani and Jain, 2010; Thomas and Singh, 2019). Despite the evidence of the profound negative impacts of excessive fertilizer application on crops and their soil systems, environmental and human health, numerous barriers complicate solving this problem. For example, Zhang et al. (2020) found that 86% of farmers with high environmental awareness who overapply mineral fertilizers nevertheless noted a low likelihood of reducing their fertilizer use. Our study addresses the need to advance fertilizer alternatives in line with several of the Sustainable Development Goals to reliably deliver nutrients to crops, avoid pollution and enable the circular nutrient economy.

Organic fertilizers (e.g. manures, composts, we use the term ‘biofertilizer’ here) can provide benefits beyond crop nutrition, but also carry the risk of under‐ or over‐supplying nutrients to crops. From a sustainability standpoint, biofertilizers are attractive as they are typically derived from waste streams (Delgado et al., 2016). Biofertilisers ameliorate soil with organic matter, which is an advantage over mineral fertilizers (Manlay et al., 2007). However, drawbacks relating to the rate of N release from organic substrates and subsequent mineralisation kinetic rates have called into question the sole reliance on biofertilizers for satisfying a crop's N requirements (Delgado et al., 2016). This risk of diminished crop yield is the most important factor for farmers who are less concerned about the longer‐term sustainability of farming than the immediate yields (Zhang et al., 2020). Here we address the need to develop practices that ‘increase agricultural N use efficiency (NUE) while sustaining or building soil organic matter and soil fertility’ (Udvardi et al., 2021).

The concept of bioengineering soils is debated in the context of nutrient use and soil fertility which involves augmenting the microbial communities of soils and crops. This can be achieved in several ways including adding plant growth‐promoting rhizobacteria (PGPR) to crop systems, where they can exert benefits on plants such as increased resilience against biotic (pest, pathogens) and abiotic stresses (e.g. heat, drought, toxicity) (Goswami et al., 2016; Vejan et al., 2016; Singh et al., 2020). PGPR that produce phytohormones and other biochemicals can increase nutrient mobilization for immediate crop use and consequently improve agricultural productivity (Adesemoye et al., 2008; Lugtenberg and Kamilova, 2009; Abbasi et al., 2011; Qiu et al., 2019). In addition, PGPR performing secondary mechanisms of action (e.g. defence against phytopathogens, inhibition of plant stress responses) are harnessed as biocontrol agents (Naseby et al., 2000; Yang et al., 2014; Arya et al., 2018) to augment the rhizosphere microbiome by reducing the presence and/or effectiveness of pathogens in crops, including sugarcane and other crops (Fu et al., 2017; Hamonts et al., 2018; Araujo et al., 2019; Elsayed et al., 2020). There are solid evidences that PGPR biofertilizers benefit plants in experimental conditions and some of these aforementioned key mechanisms are well understood (Taulé et al., 2012; dos Santos et al., 2020), but the effects of commercial biofertilizers have rarely been examined in crop field situations (or if they have, this information is often not in the public domain). Therefore, such steps into effective real‐world application demand further research and development (Finkel et al., 2017).

Here, we use a multi‐factorial design with sugarcane grown at two commercial farms to examine if biofertilizer (organic waste amended with mineral fertilizer and added PGPR) benefits crops or crop systems. Previously we found that pasture grass receiving fertilizer combinations (50:50 mixture of mineral and organic fertilizer) and supplemented with PGPR matched crop growth with mineral fertilizer alone (Paungfoo‐Lonhienne et al., 2019). Importantly, biofertilizer reduced N leaching by 95% compared to mineral fertilizer, which supports the concept of a hybrid organo‐mineral‐PGPR fertilization regime that achieves similar yields as mineral fertilizers but reduces the environmental harm (Paungfoo‐Lonhienne et al., 2019). Similarly, increasing the crop nutritional value of biofertilizer and reducing the use of mineral fertilizers benefitted soil chemical–biological properties and crop yields (Kaur et al., 2005; Chivenge et al., 2011). However, the way in which such blended biofertilizers affect crops in specific situations demands research so that farmers can access information to make judicious decisions. We explored the effect of biofertilizer, with or without added PGPR, in sugarcane. Sugarcane is a major tropical crop with high nutrient needs, and is one of the fastest growing crops (producing up to 240 t of biomass per hectare per year). Sugarcane is often grown in high rainfall coastal regions, which contributes to the leaching of N into the deeper soil, ground and surface water, causing damage to local aquatic ecosystems which are becoming a growing concern. In response, regulators are increasingly restricting the application of mineral fertilizers on farms, leading to the need for alternative N fertilization regimes that can generate the desired yields with a reduced pollution footprint (Robinson et al., 2011; Udvardi et al., 2021). We hypothesised that the combination of fast (mineral) and slow (organic) N supply in biofertilizer improves NUE and that the addition of PGPR will have demonstrable effects on the composition of soil microbial communities and crop performance.

Results

Effects of biofertilizer ‘EcoNPK’ and PGPR ‘SOS3’ on sugarcane growth and sugar yield

seven treatments were designed (Table 1) Experiments were performed at two commercial sugarcane farms and seven treatments were designed (Table 1) to evaluate the effect of organic fertilizer. Sugarcane biomass ranged from 61.7 to 107.8 t ha−1 in Ingham and 71.9 to 106.9 t ha−1 in Innisfail (Fig. 1). Sugar yield ranged from 8.1 to 14.9 and 9.7 to 15.7 t ha−1 in Ingham and Innisfail respectively. Biomass and sugar yields were significantly higher in Innisfail than in Ingham (P < 0.01). In Ingham, all fertilized and control treatments had statistically similar (P > 0.05) biomass and sugar yields. In Innisfail, treatments 120NB and 72NE had significantly (P < 0.05) higher biomass and sugar yield than the zero N treatment (0N) (Table S2).

Table 1

The experimental design comprised of two N rates, with the industry recommended rate (100%, 120 kg N ha−1) and a 40% reduced rate (60%, 72 kg N ha−1).

Treatment	Description	Treatment	Conventional fertilizer	SOS3 PGPR	EcoNPK
			(kg N ha−1)	(kg ha−1)	(kg N ha−1)
1	No‐N	0N
2	Mineral N (100%)	120N (100%N)	120
3	Mineral N (100%) + PGPR	120NB (100%N)	120	94
4	Mineral N (60%)	72N (60%N)	72
5	Mineral N (60%) + PGPR	72NB (60%N)	72	94
6	Mineral N (30%) + EcoNPK (30%)	72NE (60%N)	36		36
7	Mineral N (30%) + EcoNPK(30%) + PGPR	72NEB (60%N)	36	94	36

Organic–inorganic fertilizer EcoNPK was applied at final N input equal to mineral fertilizer, while SOS3 PGPR product (107 cells g−1) was applied at the same time in the corresponding treatments.

Fig. 1

Cane and sugar yield from Ingham (A and C) and Innisfail (B and D). No significant difference was found in either cane or sugar yield in Ingham (P > 0.05), while significant differences were found between 0N and 120NB, 0N and 72NE, 120NB and 72NB from cane yield, 0N and 120NB, 0N and 72NE from sugar yield in Innisfail (P < 0.05). Colour bars indicate different fertilizing treatments.

Effect of EcoNPK and SOS3 on soil and plant physiochemical parameters

To evaluate the effects of the various fertilizer treatments on soil N relations, we analysed the soluble soil mineral N to 1 m soil depth over time to quantify the release of immediately crop‐available N and its potential transport through the soil profile. Overall, the main difference occurred between the no‐N control and the fertilized treatments. In Ingham, soluble ammonium and nitrate concentrations changed over time following fertilization and differed between time points and soil depths (Fig. 2). Ammonium concentrations in topsoil (0–20 cm) steadily decreased over time in 120N and 72N treatments. EcoNPK treatments (72NE and 72NEB) had increased ammonium concentrations (P < 0.1 or 0.05) between T1 (1 month) and T3 (3 months) after fertilization which then decreased (Fig. 2A). Topsoil nitrate increased in all fertilized treatments (120N, 72N, 72NE, 72NEB, Fig. 2A) in the first 5 months before declining sharply at T4 (6 months). Ammonium in the zero N control was significantly lower than all other treatments at T1 and T2 (P < 0.05), and the 120N treatment was significantly higher than 72 NE (P < 0.05) and elevated over 72NEB (P < 0.1, Fig. 1A). At T3, topsoil (0–20 cm) ammonium of the control was significantly (P < 0.05) lower than 120N, 72N and 72NE, while all other depths and sampling times were similar (Fig. 2B). Nitrate in the control was significantly lower than 120N in T1, 120N and 72NEB in T2, and lower than all treatments in T3 topsoil (0–20 cm, P < 0.05). Nitrate in the control was also significantly lower than 72NEB at T3 at 40–60 cm depth, and 120N at T4 at 20–40 cm depth (P < 0.05, Fig. 2B). Furthermore, nitrate was significantly lower in 72N than in 72NEB at T3 at 40–60 cm depth (P < 0.05, Fig. 2B).

Fig. 2

Soil mineralised nitrogen at different time points and at different depths.

A. Nitrogen content (mg kg−1, square‐rooted) of ammonium nitrogen (NH4‐N) and nitrate–nitrogen (NO3‐N) in topsoil (0–20 cm) at 1 month (T1), 3 months (T2), 5 months (T3), 6 months (T4) and 12 months (T5) after fertilization.

B. Nitrogen content (mg kg−1, square‐rooted) of ammonium nitrogen (NH4‐N) and nitrate nitrogen (NO3‐N) at 0–20, 20–40, 40–60 and 60–100 cm of soil depth at 5 (T3), 6 (T4) and 12 months (T5) after fertilization. Level of significance: *P < 0.05, **P < 0.01, ***P < 0.001. Colours indicate different fertilizing treatments.

In Innisfail, ammonium concentrations decreased from T1 to T2 across all treatments. At T3, ammonium in 72N and control treatments increased. EcoNPK groups (72NE and 72NEB) remained at the same level while 120N continued to decrease (Fig. 2A). Similar to Ingham, ammonium decreased in all treatments at Innisfail after T3 until the end of the experiment. Nitrate in all fertilized treatments increased from T1 to T2 and then decreased until the end of the experiment, except for the control which increased at T3 followed by a decrease (Fig. 2A). Overall, the no‐N control had significantly lower ammonium than all other treatments in T1, as well as 120N and 72N in T2 (P < 0.05), while 72NE and 72NEB had slightly higher ammonium than the control in T2 (P < 0.1). Differences in nitrate concentrations only occurred between the control and 120N (P < 0.05), and 72N (P < 0.1) during T1, whereas the control had significantly lower nitrate concentrations than all other treatments during T2 (P < 0.05, Fig. 2A). Soil depth did not feature different ammonium or nitrate levels across treatments. All statistical results were presented in Table S3.

Microbial assemblages in response to EcoNPK and SOS3

After rarefying to equal sampling depth per sample, a total number of 4436 unique bacterial operational taxonomic units (OTUs) and 1940 unique fungal OTUs were analysed in this study. In Ingham, slight significant differences were found between 0N and 120N in the bacterial Shannon index as well as 0N and 72NEB in fungal Shannon index (P < 0.1), while both bacterial and fungal communities had similar Chao1 or Simpson indices (P > 0.05, Fig. S2). In Innisfail, no significant differences were found in either bacterial and fungal communities in Chao1, Shannon or Simpson indices (P > 0.05, Fig. S2).

To evaluate the effects of PGPR SOS3 and organic fertilizer EcoNPK on microbial communities, samples were pooled according to their treatments (application of SOS3 or EcoNPK), and the resulting impacts on microbial compositions and abundances were analysed. Overall, there was no clear pattern of variations in bacterial communities under SOS3 or EcoNPK treatments on either site (Fig. S3A–D). However, the addition of SOS3 slightly increased Basidiomycota occurrence in fungal communities compared to the control, which was consistent on both sites (Fig. S3E and F). There was no observable pattern in fungal community changes in EcoNPK treatments (Fig. S3G and H). PERMANOVA results and non‐metric multidimensional scaling (NMDS) plots showed that at both sites, time was the major factor that significantly impacted bacterial and fungal abundances and compositions (P < 0.05, Table 2a; Fig. S4). In Ingham, SOS3 and EcoNPK had no significant effects on overall bacterial or fungal communities (P > 0.05, Table 2). In Innisfail, SOS3 treatments had no significant effect on bacterial or fungal communities between fertilized treatments (P > 0.05, except C ≠ B in month 3, fungal community), while significant difference occurred between EcoNPK and mineral fertilizer on bacterial communities (C ≠ E, P < 0.05) and fungal communities (V ≠ E, P < 0.05 in month 1 and 6, P < 0.1 in month 3, Table 2; Fig. S4).

Table 2

Pairwise PERMANOVA of bacterial (16S) and fungal (ITS) communities based on Bray–Curtis and measures of square‐root transformed relative abundances under different treatments (C = control, S = SOS3, V = conventional fertilizer, E = EcoNPK, NSD = no significant difference).

Location	Ingham				Innisfail
Treatments	SOS3		EcoNPK		SOS3		EcoNPK
Source	16S	ITS	16S	ITS	16S	ITS	16S	ITS
Time	<0.001	<0.001	<0.001	<0.001	<0.001	<0.001	<0.001	<0.001
Treatments	0.399	0.476	0.169	0.320	0.375	0.720	0.035	<0.001
Time × Treat	0.925	0.759	0.996	0.613	0.422	0.009	0.069	<0.001
Pairwise analysis
T1	NSD	NSD	NSD	NSD	NSD	NSD	V ≠ E	C, V ≠ E
T2						C ≠ S		C ≠ E
T4						NSD		V ≠ E

Significant P values (P < 0.05) are highlighted in bold.

Using linear discriminant analysis effect size (LEfSe) analysis, we identified a suite of bacterial and fungal OTUs as distinct microbial markers, indicating the most prominent microbial taxa present under SOS3 and EcoNPK treatments. In the bacterial communities, six OTUs were associated with SOS3 treatments in Ingham in the phyla Actinobacteria, Bacteroidetes, Chloroflexi and Proteobacteria, while 12 OTUs from phyla Actinobacteria, Firmicutes, Thermodesulfobacteria and Proteobacteria were identified in EcoNPK treatments in Ingham (Fig. 3A). In Innisfail, seven OTUs were associated with SOS3 treatments, with five OTUs in the phylum Actinobacteria and two in the phylum Proteobacteria, and 11 OTUs were associated with EcoNPK treatments, in phyla Actinobacteria, Chloroflexi, Firmicutes and Proteobacteria (Fig. 3B). In the fungal communities, three OTUs from the order Capnodiales were markers in Ingham EcoNPK treatments (Fig. 3C), and seven OTUs from the orders Saccharomycetales, Hypocreales, Microascales and Ustilaginales in Innisfail EcoNPK treatments (Fig. 3D). No unique fungal marker was detected in SOS3 treatments (Fig. 3C and D). A detailed list of microbial markers is presented in Table S4.

Fig. 3

Cladogram of bacterial (16S, A and B) and fungal (ITS, C and D) communities based on LEfSe analysis in Ingham (A and C) and Innisfail (B and D). Microbial markers at different taxonomic levels were highlighted with colours based on treatments (SOS3: blue, EcoNPK: green, conventional fertilizer: purple, Control: red). Full details of microbial taxa labelled with alphabet and numbers were provided in Table S4.

Soil physiochemical parameters in correlation with microbial communities

To further explore the variation of microbial communities in the context of soil factors, we calculated the correlation between microbial richness (Chao1) and environmental parameters including soluble N, pH, electrical conductivity (EC) and moisture using a Spearman correlation. In Ingham, soluble N and EC significantly increased bacterial richness (R 2 = 0.287, 0.096 and 0.138 respectively, Table 3; Fig. S5), while soil pH and moisture significantly reduced fungal richness (R 2 = 0.185 and 0.014, Table 3; Fig. S5). In Innisfail, the increase of soil ammonium, pH and moisture significantly increased bacterial richness (R 2 = 0.147, 0.291 and 0.132, Table 3; Fig. S5), while fungal richness was unaffected.

Table 3

Spearman correlation between microbial richness and soil physiochemical parameters in (a) Ingham and (b) Innisfail (EC = Electrical conductivity).

	16S			ITS
	R 2	P	rho	R 2	P	rho
a) Ingham
Ammonium	0.287	<0.001	0.346	0.024	0.971	−0.196
Nitrate	0.096	0.002	0.528	0.032	0.095	−0.004
pH	0.0002	0.985	−0.002	0.185	<0.001	−0.424
EC	0.138	<0.001	0.370	0.013	0.133	−0.176
Moisture	0.006	0.409	−0.094	0.014	0.042	−0.237
b) Innisfail
Ammonium	0.147	<0.001	0.364	0.013	0.345	−0.104
Nitrate	0.110	0.187	−0.149	0.008	0.671	0.047
pH	0.291	<0.001	0.545	0.005	0.376	−0.098
EC	0.005	0.274	0.124	0.0001	0.919	−0.011
Moisture	0.132	<0.001	0.410	0.002	0.451	−0.083

Significant P‐values (P < 0.05) are highlighted with bold.

DistLM analysis of soil microbial communities revealed that in Ingham, most of the soil parameters except TN were correlated with both bacterial and fungal compositions and abundances, explaining 34.94% and 23.47% of total variation respectively (Table 4a). More specifically, ammonium, nitrate, EC and pH had high contributions (9.67%, 7.48%, 7.94% and 5.01% respectively) to bacterial communities, while soil moisture, ammonium, nitrate and pH contributed substantial variation (6.34%, 6.29%, 2.86% and 3.62% respectively) to fungal communities. TC significantly drove the variation of fungal communities (P = 0.03), but the contribution was relatively low (1.82%). Together, most of the parameters significantly contributed to soil microbial community variations except TN and TC (Table 4a).

Table 4

DistLM analysis of multivariate data based on Bray–Curtis dissimilarity and environmental parameters of (a) bacterial and (b) fungal communities for specified variables.

Variable	R 2	F	P	Var.%	Variable	R 2	F	P	Var.%
a) Ingham					b) Innisfail
16S					16S
Moisture	0.024	1.918	0.014	2.40	Moisture	0.043	3.508	<0.001	4.30
Ammonium	0.096	8.406	<0.001	9.67	Ammonium	0.044	3.847	<0.001	4.70
Nitrate	0.013	1.152	<0.001	7.48	Nitrate	0.032	2.778	<0.001	3.44
EC	0.012	0.994	<0.001	7.94	EC	0.019	3.612	<0.001	4.43
pH	0.035	3.232	<0.001	5.01	pH	0.028	3.480	<0.001	4.27
TC	0.011	0.984	0.355	1.30	TC	0.014	1.025	0.347	1.30
TN	0.011	1.012	0.583	1.14	TN	0.009	0.905	0.611	1.15
ITS					ITS
Moisture	0.063	4.875	<0.001	6.34	Moisture	0.025	2.010	0.028	2.51
Ammonium	0.063	5.109	<0.001	6.29	Ammonium	0.119	10.660	<0.001	11.86
Nitrate	0.029	2.369	<0.001	2.86	Nitrate	0.047	4.384	<0.001	4.67
EC	0.014	1.166	0.145	1.40	EC	0.026	2.405	0.007	2.52
pH	0.036	3.097	<0.001	3.62	pH	0.013	1.278	0.155	1.33
TC	0.018	1.566	0.003	1.82	TC	0.009	0.876	0.627	0.91
TN	0.012	0.983	0.495	1.14	TN	0.008	0.807	0.706	0.84

Var.% = variations explained by each parameter, EC = electrical conductivity. Significant P‐values (P < 0.05) are highlighted in bold.

In Innisfail, all parameters explained 23.59% and 24.64% of total bacterial and fungal community variations respectively (Table 4a). Similar to Ingham, most of the parameters significantly contributed to variations in both communities except TC and TN. Soil moisture, ammonium, nitrate, EC and pH explained 4.30%, 4.70%, 3.44%, 4.43% and 4.27% variation in bacterial communities respectively, while ammonium was the most important factor driving variation of fungal communities (11.86%) followed by nitrate (4.67%). pH did not significantly affect fungal communities (Table 4b). Collectively, soil moisture, ammonium and nitrate are the most important factors that crucially influenced both bacterial and fungal communities. EC and pH significantly drove bacterial communities but only affect fungal communities on specific sites. Generally, TC and TN did not provide high contributions to soil microbial communities due to the high soil background of TC and TN.

In the dbRDA analysis, a shift of the microbial communities was observed over time, in which variations were mainly affected by pH, EC, NH4 + and NO3 − (Fig. 4). In Ingham, pH, EC and moisture were correlated closely with bacterial and fungal communities at T1, while NH4 + and NO3 − had strong correlations with bacterial and fungal communities at T2 and T4 (Fig. 4A and B). In Innisfail, bacterial and fungal communities at T1 were closely correlated with NH4 +, followed by pH and EC, while NO3 − and EC had a stronger correlation with fungal communities at T2 (Fig. 4C and D).

Fig. 4

Distance‐based RDA (dbRDA) analysis of bacterial (A and C) and fungal (B and D) communities in correlation with ammonium nitrogen (NH4‐N), nitrate–nitrogen (NO3‐N), pH, EC and soil moisture at 1 month (T1, orange), 3 months (T2, green) and 6 (T4, blue) months after fertilization. Different shapes indicate nitrogen input in the fertilization, and solid/hallowed symbols indicate additional input of SOS3 (solid: SOS3, hallowed: SOS3‐free).

Bacterial–fungal networks

In the co‐occurrence network analysis, bacterial and fungal OTUs with strong correlations (|r| > 0.8) were extracted for a network analysis (Fig. 5). Bacterial taxa with strong correlations dominated in Ingham, while fungal species with strong correlations dominated in Innisfail (Fig. 5). Higher average degrees were generally found in SOS3‐free treatments at both sites when comparing under the same fertilization strategies (i.e. 120N vs. 120NB, 72N vs. 72NB and 72NE vs. 72NEB) except 72NE versus 72NEB in Ingham (Table 5; Fig. 5A). Additionally, higher positive rates in correlations were found in SOS3 treatments when comparing the same fertilization regimes except 120N versus 120NB in Ingham (Table 5; Fig. 5A). In intra‐ and inter‐kingdom correlations, higher proportions of bacterial–bacterial correlations were commonly found in Ingham, especially in fertilized treatments (all treatments except 0N, Table 5a). However, higher proportions of fungal–fungal correlations were commonly found in Innisfail (Table 5b, Fig. 5B). Relatively lower proportions of fungal–fungal correlations were found in the EcoNPK treatments in Innisfail (72NE and 72NEB), compared to other treatments (Table 5b). In Innisfail soil, a core fungal OTU group was consistently present in the first five treatments (0N, 120N, 120NB, 72N and 72NB). These core OTUs were tightly and positively connected in EcoNPK‐free treatments, while such association was disrupted in EcoNPK treatments (nodes in red, Fig. 5B). This pattern was even stronger in EcoNPK + SOS3 treatment, where a group of bacteria was established to replace the fungal association (72NEB, Fig. 5B). Although there was no strong fungal association observed at Ingham, a similar pattern in the bacterial network was also detected with 72NEB (Fig. 5B).

Fig. 5

Co‐occurrence network analysis of bacterial and fungal communities in Ingham (A) and Innisfail (B) under different fertilization treatments. Colours of nodes indicate bacterial (blue) and fungal (non‐core: brown, core: red) OTUs, and colours of edges indicate positive (green) and negative (red) correlations. The size of nodes indicates the weight of the corresponding OTU (numbers of edges connected), and the weight of the edges indicates the degree of the correlation.

Table 5

Bacterial and fungal co‐occurrence networks in the seven treatments.

	1	2	3	4	5	6	7
Treatment	0N	120N	120NB	72N	72NB	72NE	72NEB
a) Ingham
Nodes	193	260	143	114	103	248	250
Edges	289	483	213	167	135	355	492
Avg. degree	3.00	3.72	2.98	2.93	2.62	2.86	3.94
Positive percentage	77.16	69.77	68.54	73.05	82.22	61.69	70.33
B‐B percentage	53.29	73.29	70.42	73.65	74.81	59.15	79.67
B‐F percentage	33.22	22.57	16.43	16.77	11.11	29.58	14.84
F‐F percentage	13.49	4.14	13.15	9.58	14.07	11.27	5.49
b) Innisfail
Nodes	115	96	102	78	93	193	114
Edges	315	283	185	280	189	429	217
Avg. degree	5.48	5.90	3.63	7.18	4.07	4.45	3.81
Positive percentage	89.52	85.16	99.46	74.29	96.3	70.86	75.58
B‐B percentage	22.54	20.49	18.38	5.36	15.34	27.51	46.54
B‐F percentage	12.70	15.55	3.78	13.21	3.17	29.84	10.14
F‐F percentage	64.76	63.96	77.84	81.43	81.48	42.66	43.32

Nodes (OTUs) and edges (correlations) with only strong correlations (|r| > 0.8) were kept in the results. B‐B = bacterial–bacterial correlations; B–F = bacterial‐fungal correlations; F–F = fungal‐fungal correlations.

Discussion

Biofertilizer can partially replace mineral nitrogen fertilizer

Nitrogen is one of the most important factors that drive growth and productivity of non‐leguminous plants (Below, 2001; Plett et al., 2020), and poor N fertilizer use efficiency in sugarcane production is one of the main hurdles that must be overcome for sugarcane to reach its full potential as a bioeconomy crop (e.g. Grant et al., 2016; Rodriguez et al., 2019; Wang et al., 2019). Biofertilisers have drawn attention as a sustainable alternative to mineral fertilizers in sugarcane and other crops (del Carmen Rivera‐Cruz et al., 2008; Molla et al., 2012; Nunes Oliveira et al., 2017; Yadav and Sarkar, 2019). As most plant‐associated microorganisms have intimate relationships with their hosts (Bonfante and Anca, 2009), biofertilizers containing effective PGPR can increase nutrient use efficiency by mineralizing, solubilizing, mobilizing and supplying N to plants (Maeder et al., 2002; Qiu et al., 2012; Arif et al., 2017; Pereira et al., 2020). PGPR that produce beneficial phytohormones, fix nitrogen and solubilize phosphorus can further enhance nutrient delivery to hosts (dos Santos et al., 2020; Grover et al., 2021; Verma et al., 2021). Nutritional benefits of microbes may also include being themselves consumed by the roots, through a supposed mechanism currently under investigation termed ‘rhizophagy’ (Paungfoo‐Lonhienne et al., 2010; White et al., 2018).

Irrespective of microbial taxa and mechanisms, in field growth conditions, PGPR increased the yields of cereals, root crops, legumes and vegetables (Ali et al., 1998; Tena et al., 2016; Islas‐Valdez et al., 2017; Mukhongo et al., 2017). In contrast, yields of field‐grown sugarcane were unaffected by PGPR, although the rhizosphere fungal community was altered (Berg et al., 2019).

Here, mineral or biofertilizers (with/without added PGPR) generated similar amounts of soluble ammonium and nitrate in the top meter of soil, confirming that both fertilizers deliver similar amounts of readily crop‐available N to sugarcane over the growing season. Fertilizer type, or reducing the N supply by 40%, had only negligible effects on the detectable ammonium and nitrate, while omitting N fertilizer resulted in a noticeable decline in soil N (but this did not impact crop yield). That crops were N sufficient irrespective of N application rate is unsurprising as sugarcane draws on soil N stores with a considerable portion of soil organic N mobilized over a crop season (Allen et al., 2019). Similar soluble N levels were observed with mineral or biofertilizer and supplied N to pasture grass (Paungfoo‐Lonhienne et al., 2019). Thus, the biofertilizer tested here underpins the global aim of transitioning to a circular N economy with a proportion of N retrieved from wastes. To what extent biofertilizer can increase NUE and/or build soil organic matter and soil fertility requires longer‐term research.

Biofertilizer and PGPR can augment soil microbial communities

Application of biofertilizer significantly changed soil microbial communities with fungal communities altered at certain cropping stages at the Innisfail site. This aligns with previous studies demonstrating that biofertilizers can alter soil microbiota (Paungfoo‐Lonhienne et al., 2015; Dong et al., 2019) through interactions with plant roots (Azizoglu, 2019; Trivedi et al., 2020) and competition with potential pathogens (Zhang et al., 2018). Similar to previous research, we find that the impact of biofertilizers on soil microbial communities is site‐specific and depends on soil and climate properties, as well as local microbial biodiversity (Ma et al., 2016; Delgado‐Baquerizo et al., 2018; Xue et al., 2018). In the top 20 cm of soil (where microbial communities were analysed in our study), the two soils differed considerably in texture and organic matter content. Ingham sandy loam soil had considerably less soil organic matter (0.5%C, 0.04%N) than Innisfail clay loam (3.2%C, 0.15%N). While soil microbial diversity was similar across all treatments, we detected variation in microbial abundance and composition between treatments and regions. The identified microbial markers between sites and biofertilizer application confirm that changes of soil microbiota are region and/or soil specific.

Similarly, soil microbiota were region‐specific and possibly associated with soil physiochemical properties among other factors (Ng et al., 2014; Ma et al., 2016; Yeoh et al., 2016). While soil pH is often reported as a primary driver of soil microbial communities (Fierer and Jackson, 2006; Lauber et al., 2009; Hong et al., 2020), other factors, including mineralised N (Paungfoo‐Lonhienne et al., 2015; Zhou et al., 2017), soil moisture (Brockett et al., 2012) and salinity (Ren et al., 2018) were closely related with soil microbial assemblages. We found that Spearman's correlation, dbRDA and DistLM analyses confirmed that, despite the variations in microbial responses to soil parameters, ammonium had a consistent positive correlation with bacterial richness, highlighting that soil N status and microbial diversity are linked. Most of the five parameters investigated (soil moisture, ammonium, nitrate, EC, pH) contributed to the variation of microbial communities at both sites, suggesting a close relationship of soil physiochemical parameters and microbial dynamics. The decrease of soluble N (especially ammonium) at both sites explained the significant difference in microbial communities at different time points, supporting similar trends of shifts in soil microbiota in two contrasting soils. Similarly, the application of biofertilizer with lower N input can also achieve alteration of the soil microbiota due to changes in the soil physiochemical properties.

Engineering soil microbiomes is receiving much attention in the quest for sustainable cropping (Mueller and Sachs, 2015; Qiu et al., 2019). While the abundance of inoculated PGPR can decrease over time, lasting effects can include that microbial networks are enhanced by establishing core and hub microorganisms that in turn promote microbial functions, crop fitness and growth (Agler et al., 2016; Lemanceau et al., 2017; Toju et al., 2018; Zhuang et al., 2021). Application of PGPR can contribute to sustainable cropping (Ma et al., 2016; Li et al., 2017; Sun et al., 2020), although efficacy has to be confirmed for a given crop system including sugarcane (Berg et al., 2019).

Biofertilizer can establish stronger bacterial associations in soil

Microbial interactions in soils critically influence plant performance and productivity (Trivedi et al., 2017; Trivedi et al., 2020). Positive interactions are usually considered as mutual cooperation between organisms, while negative interactions can be competitors that exert antagonistic mechanisms (Abhilash et al., 2012; Hassani et al., 2018; Araujo et al., 2019). While network analysis cannot truly reflect microbial interactions, it highlights associations based on co‐occurrence of microbial species and as such enables insight into microbial associations at an inter‐kingdom scale (Li and Wu, 2018; Zheng et al., 2018). The general concept is that positive associations often occur on an intra‐kingdom basis, whereas negative associations are more commonly observed on an inter‐kingdom scale (Sweet et al., 2019). Our results showed that the PGPR inoculant (a Paraburkholderia species), added to mineral N and biofertilizer, can increase the rate of positive associations (especially bacterial–bacterial associations) in most of the scenarios compared to treatments without PGPR addition. This suggests that the PGPR strain used here can enhance bacterial networks in the soil microbiome. Augmenting positive correlations might link to ecological functions of microbial communities, for example carbon and nitrogen cycling (Hoppe et al., 2014; Banerjee et al., 2016). It may have positive effects on soil biological function and promote beneficial microbial associations in the rhizosphere as observed with other PGPR in previous research (Ling et al., 2016; Fan et al., 2017; Chen et al., 2018), but longer term research with quantifiable response variables is needed to support this notion.

PGPR‐strengthened microbial networks generated more associations. We found that the biofertilizer established a stronger bacterial network in replacement of the fungal network, which was particularly prominent in Innisfail samples, possibly due to strong fungal associations pre‐existing in the soil. The higher soil organic matter content in the Innisfail soil may be a key factor contributing to these stronger networks. Although fungal communities can be beneficial or neutral for crops (Bonfante and Genre, 2010), fungi are more often reported as crop pathogens (Deacon, 2005). Biofertilisers can increase plant fitness and protect crops from infections of fungal pathogens by establishing stronger bacterial associations with better resilience to environmental stresses (Siddiqui, 2005). Whether this is the case in our study – and which organisms are involved – requires further investigation.

The input of organic matter can increase fungal–bacterial competition and alter community structure in the soil microbiome to improve microbially driven ecosystem functions (Ling et al., 2016; Zheng et al., 2018; Gu et al., 2019) and boost plant growth and resilience to biotic and abiotic stresses. Sugarcane soil with lower mineral N loads favoured beneficial fungi compared to soil with a higher N load (Paungfoo‐Lonhienne et al., 2015). Together, these observations deserve further investigation. The findings here certainly confirm that biofertilizer and PGPR affect soil and crop responses, and that they can benefit soil function and sugarcane production in certain situations.

Conclusions

Our two main findings are that a biofertilizer with a 40% reduced N load formulated with 50:50 N (recycled waste:mineral fertilizer) was superior to mineral fertilizers at one of the two studied sites. The second finding was that this biofertilizer combined with PGPR enhanced bacterial and fungal associations with sugarcane which may benefit crops as this occurred at the same site where the biofertilizer appeared to benefit yield. Both findings warrant further investigations with long‐term field testing, including at sites where the soil microbiome negatively impacts crop yield. The next step for this research should explore how NUE can improve with biofertilizer application to reduce N losses from soil. With so much at stake for crop production and environmental integrity, efforts should be heightened to optimize N supply for sugarcane as an critical tropical bioeconomy crop.

Experimental procedures

Experimental design

The research was performed at two commercial sugarcane farms in Australia's Wet Tropics (Queensland, Australia), one near Ingham (−18.65067° N, 146.15472° E) on a sandy loam soil, and the other near Innisfail (−17.52448° N, 146.02815° E) on a clay loam soil. The crop at both locations was rainfed (2144 and 3759 mm average annual rainfall at Ingham and Innisfail respectively, Bureau of Meteorology) and consisted of plant cane (the start of a multi‐year ratoon cycle) of commercial variety Q208. The blocks had previously been fallowed without addition of organic amendments for 6 months following a former sugarcane crop. Before sugarcane was planted, soil was collected on both sites at four depths (0–20, 20–40, 40–60, 60–100 cm) to record the initial soil physiochemical parameters.

At each site, plots within the field were established in a randomized block design with each plot consisting of six rows of sugarcane (20 m long, 1.8 m row spacing) and four replicated plots per treatment. To evaluate the effect of the organic fertilizer (EcoNPK, containing 4.0% N, 1.6% P, 1.9% K with C:N 5.3:1.0) and microbial products SOS3 PGPR bacterium (Paraburkholderia sp., NCBI ID: 1926494; both products from Sustainable Organic Solutions Pty., Brisbane, Australia), seven treatments were designed (Table 1). All fertilizers and the microbial products were applied at planting to an open furrow close to the cane setts. Nitrogen was added as conventional urea (120 or 72 kg N per hectare for 100%N or 60%N treatments) or as an N‐equimolar combination of urea and poultry manure‐based organic EcoNPK fertilizer. Superphosphate, Muriate of Potash and SuPerfect were applied in all plots as a source of P, K and S respectively. All treatments received the same amount of P, K and S (30 kg P, 100 kg K, 25 kg S per hectare).

Sample collection

After 12 months of growth (August 2019 to August 2020 in Ingham, September 2019 to September 2020 in Innisfail), sugarcane was hand‐harvested by cutting 10 (Ingham) and 30 (Innisfail) stalks from the middle two rows of all replicate plots of all treatments. Stalks were weighed and processed at industry service provider HCPSL Macknade to obtain three yield measures: tonnes of cane, sugar yield and commercial cane sugar yield (CCS, i.e. recoverable sugar content). Yield per hectare was estimated by multiplying the weight of harvested stalks by the total number of stalks in each plot in line with long‐term industry practice.

To determine soil mineral N (ammonium, NH4 +; nitrate, NO3 −) dynamics, soil samples were collected 1 (T1), 3 (T2), 5 (T3), 6 (T4) and 12 months (T5) after fertilizer application. In the first two time points, soil from 0 to 20 cm depth was collected from three positions on the fertilizer band and pooled. To evaluate leaching of mineral N, soil cores to a depth of 100 cm were sampled from the fertilizer band at the latter three sampling time points. These deep cores were collected using a manual lever and occurred after heavy rainfall events (>65 mm over 3–4 day period or >40 mm over 24 h). Three cores were sampled per plot and bulked into four depths (0–20, 20–40, 40–60 and 60–100 cm) before being transported to the laboratory in insulated boxes with ice blocks and stored at 4°C until further processing.

For the microbial community analysis, soil samples were collected from the topsoil (0–20 cm) on the first day (T0), 1 month (T1), 3 months (T2) and 6 months (T4) after fertilizer application. Briefly, approximately 100 g of soil was collected 1–10 cm from roots and aseptically transferred into sample bags and transported to the laboratory using the method described above. Samples were stored at −80°C until further processing.

For leaf chemical analysis, leaf samples were collected 6 and 12 months after planting. Ten leaves (the third leaf from the top of the stalk) were collected from the middle two rows of each replicate plot.

Determination of soil and plant physiochemical parameters

The concentration of mineral N (ammonium, nitrate) in soil samples was determined via spectroscopy after extracting a 20 g subsample of mixed field‐moist soil with 2 M potassium chloride following the protocol of Rayment and Lyons (2011). Soil gravimetric moisture content was determined by drying a sub‐sample of soil at 105°C for 48 h. Total organic C and N contents in soil were determined using fine ground air‐dried sub‐samples with the Dumas combustion method (LECO CN928 Analyser, Michigan, USA). Primary particle size distribution was determined using the hydrometer method (Thorburn and Shaw, 1987). A subsample of field moist soil was dried at 40°C for 48 h, before being sieved to determine pH and EC in a 1:5 soil water suspension using calibrated electrodes (Rayment and Lyons, 2011). All soil variables were expressed on a dry soil basis.

Leaf samples were dried at 65°C for 48 h, before grinding to <1 mm. The total N content of the ground leaf samples was determined by combustion (LECO CN928 Analyser).

DNA extraction and Illumina MiSeq sequencing

Frozen stored soil (~250 mg dry weight) was weighed and soil microbial DNA was extracted using DNeasy PowerSoil Kit (Qiagen, Hilden, Germany), following the manufacturer's instructions. Extracted DNA was quality checked by NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA), DNA quantity was checked by Qubit Fluorometer (Thermo Fisher Scientific) and PCR‐checked to confirm the amplifiability. Amplicons targeting the V3–V4 region in the 16S rRNA gene for bacterial communities (341F‐805R, Herlemann et al., 2011), and ITS1 region (ITS1F‐ITS2R, Adams et al., 2013) for fungal communities were generated. Sequencing was performed at the Australian Genome Research Facility (AGRF) using Illumina MiSeq 2 × 300 bp paired‐end chemistry. All raw sequence data related to this study are available in the NCBI Sequence Read Archive (PRJNA723115).

Sequences processing

Raw data obtained from NGS facility were processed using Mothur standard operating procedure (Schloss et al., 2009). Briefly, forward and reverse sequences were merged into contigs. Sequences that contained unidentified bases or had greater than eight homopolymers were filtered out. For bacterial sequences, an additional step aligning sequences against Silva 16S rRNA gene database version 132 (Pruesse et al., 2007) was applied, and unaligned sequences were removed. Refined sequences were pre‐clustered (diffs = 1) and chimaera checked using UCHIME (Edgar et al., 2011) and singleton was removed to reduce error (Reeder and Knight, 2009). Bacterial and fungal sequences were taxonomically classified according to the Silva database version 132 and UNITE database version 8 respectively, with 60% cut‐off confidence. Classified sequences were clustered into OTUs at 97% identity where taxonomy was assigned to. Bacterial and fungal OTUs were rarefied to equal sequencing depth per sample. Rarefaction curves of the processed sequences were asymptotic, suggesting good coverage of the microbial diversity present (Fig. S1).

Statistical analysis and data visualization

The OTU matrix was analysed using permutational multivariate analysis of variance (Anderson, 2001a) in PRIMER v. 6 (PRIMER‐E, UK) to compare bacterial and fungal communities under different treatments. Similarity matrices were calculated based on Bray–Curtis distances on square‐root transformed abundance data to compare the composition and abundances of community structure. Main analyses used 9999 permutations of residuals under a reduced model (Anderson, 2001b). Pair‐wise analyses were performed to compare the differences between treatments with significant interactions based on the main analyses. Permutational multivariate dispersion analysis was used to test for homogeneity of multivariate dispersion within groups (Anderson, 2006). Alpha and beta diversity were analysed using R package ‘phyloseq’. Data visualization, including NMDS plots, was generated based on Bray–Curtis distances. Spearman's correlation analysis between microbial richness and soil physiochemical parameters was performed with R. Distance based linear model (DistLM) and distance‐based redundancy analysis (dbRDA) were performed to evaluate the contribution of environmental factor on microbial variances (Anderson, 2008). To compare the relative abundance of each OTU across different treatments, LEfSe was used to determine the unique microbial features from different treatments (Segata et al., 2011), and taxonomic analysis was performed using R packages ‘phyloseq’, ‘dplyr’ and ‘ggplot2’ (Lozupone et al., 2012).

Co‐occurrence network analysis

To infer the bacterial and fungal associations between the OTUs and compare the associations between treatment groups, co‐occurrence networks between the OTUs defined as core microbiome above were performed using FastSpar (Watts et al., 2019), an implementation of network analysis using SparCC algorithm (Friedman and Alm, 2012). Fastspar was run using 50 iterations, 1000 bootstraps and used null models to exclude the false correlations (i.e. controlling the FDR), with statistically significant correlations (P < 0.05) kept for further analysis. Based on the refined datasets, Gephi (0.9.2) was used to visualize networks using the Fruchterman Reingold layout algorithm (Bastian et al., 2009).

Author's Contributions

Conceptualisation: C.P.‐L, L.D.B., W.W. and S.S. Methodology, data collection and original data analysis: Z.Q., J.Y., A.G.G., I.P., R.H., M.I., and S.R. Data presentation, writing, reviewing and editing: Z.Q., C.P.‐L., J.Y., L.D.B., W.W. and S.S. All authors have read and agreed to the published version of the manuscript.

Fig. S1. Rarefaction curve for the sequences of soil microbiomes obtained from (A) bacterial 16S rRNA gene sequencing and (B) fungal ITS region sequencing.

Click here for additional data file.

Fig. S2. Alpha diversity indices of bacterial (A, B, C) and fungal (D, E, F) communities from Ingham and Innisfail samples. Colours of boxplots indicate different fertilization treatments.

Click here for additional data file.

Fig. S3. Taxonomic analysis of bacterial (A, B, C, D) and fungal (E, F, G, H) communities from Ingham and Innisfail samples with SOS3 (C: control, N: SOS3‐free, B: SOS3) and EcoNPK (C: control, V: conventional fertilizer, E: EcoNPK) treatments.

Click here for additional data file.

Fig. S4. NMDS plots of bacterial (A and B) and fungal (C and D) communities from Ingham and Innisfail samples at one day (T0: purple), one month (T1: orange), three months (T2: green) and six months (T4: blue) after fertilization. Shapes of plots indicate different fertilization treatments.

Click here for additional data file.

Fig. S5. Spearman's correlation of bacterial (A and B) and fungal (C and D) richness (Chao1) from Ingham and Innisfail samples correlated with ammonium nitrogen (NH4‐N), nitrate–nitrogen (NO3‐N), pH, electrical conductivity (EC) and soil moisture.

Click here for additional data file.

Table S1. Supporting information.

Click here for additional data file.