Cyanobacteria can also form filaments called hormogonia. Hormogonium lacks heterocysts and can be morphologically distinguished from the latter as the former is capable of locomotion appearing as short motile filaments (Hernández-Muñiz and Stevens, 1987). Hormogonia play a role in selfdispersal and in forming symbiotic association during the early stages of infection by responding to chemical signals produced by prospective host plants (Khayatan et al., 2017). Chemoattractants are said to be involved in stimulating hormogonia formation and directing the cyanobacteria toward the targeted plant tissue that will house the cyanobionts (Meeks et al., 2001; Meeks and Elhai, 2002; Bergman et al., 2007).

Filamentous, heterocyst-forming cyanobacteria are also capable of differentiating into spores called akinetes (Kaplan-Levy et al., 2010; Adams et al., 2013). When environmental conditions become unfavorable (e.g., low salinity, temperature fluctuations, lack of phosphates, low light, insufficient nutrients), active cells transform into a resting state that may last for up to 60 years in which they can be restored to their vegetative cell state when favorable conditions arise (Kaplan-Levy et al., 2010). Although not resistant to intense heat, akinetes can survive cold and desiccation (Adams and Duggan, 1999). Akinetes are rare but were reported to occur in strains of Nostoc isolated from coralloid roots (Grilli Caiola, 1980; Grobbelaar et al., 1987) and more are common in Azolla-Anabaena symbioses (Peters and Perkins, 1993). Additionally, Sukenik et al. (2007) state that resting cells can still perform minimal metabolic activities such as photosynthesis, carbon fixation and protein synthesis. Aside from akinetes, lysing or dying cells called necridia, which allow excision of cells along the filament, also occur in cyanobacteria. Characterized by thickened walls and non-granulated cells, these are commonly referred to as deteriorating cells that were unable to differentiate into heterocysts (Grilli Caiola and Pellegrini, 1979).

#### General Mechanism of Symbiosis

All known plant-cyanobacterial symbioses are acquired from the environment. The only exception is seen in the fern Azolla, where the cyanobiont is an integral part of the host throughout its developmental stages and gets inherited to the next generation (Bergman et al., 1996). Thus, for most symbioses, efficient communication between the host and cyanobacterium must be carried out to ensure successful entry of the symbiont (**Figure 2**). This requires signal molecules induced by the host and/or the symbiont during the initial stages of invasion. The general mechanism based on other host plant-cyanobacterial symbioses, as controlled by a set of regulatory genes, is that the host elicits hormogonium-inducing factors causing the surrounding cyanobacteria to transform into motile hormogonia (Adams and Duggan, 2012; Warshan, 2017). Chemotactic signals permit entry into the partner plant, subsequently demanding the host to produce hormogonium-repressing factors to allow the cyanobacteria to develop heterocysts for nitrogen fixation to occur (Warshan, 2017).

Formation of motile hormogonia is essential for successful migration of a cyanobiont to the internal tissues of its potential plant partner. Certain compounds which may be produced by the host, trigger this phenomenon and initiate cell differentiation (Khayatan et al., 2017). In plants known to accommodate symbiotic cyanobacteria, modified structures can be formed regardless of whether a cyanobiont will be present in its life cycle (Meeks and Elhai, 2002) such as the case in coralloid roots. Unlike Rhizobia, which requires specialized symbiosomes to be formed by its host, cyanobionts of coralloid roots do not need similar structures to survive within its host. However, particular spaces are dedicated to symbionts. Cyanobacteria resides in the cyanobacterial zone in the roots of cycads, in the red stem glands of Gunnera, in dorsal leaf cavities of Azolla water fern, in slime cavities within the thallus of hornworts, in auricles underneath the thallus of liverworts, and in the bladders of the fungi Geosiphon pyriformis (Adams et al., 2013). Among these plant-cyanobacterial symbioses, it is only in Gunnera where the cyanobiont invasion occurs intracellularly (Meeks and Elhai, 2002; Adams et al., 2013). In the coralloid roots, an acidic viscous mucilage in the cortex (Grilli Caiola, 1980), also observed by the authors, is present in the cortex layer probably involved in attracting motile hormogonia filaments as this was the case observed in Gunnera (Nilsson et al., 2006). Mucilage-filled cavities were also observed in hornworts and liverworts and this heat-labile putative signal inducer is yet to be characterized, but was estimated to be a low molecular mass protein around 12 kDa (Adams et al., 2013). In line with this, root extracts from cycad coralloid roots were found to significantly initiate hormogonia formation in strains of Nostoc (Campbell and Meeks, 1989; Bergman et al., 1996; Meeks, 1998) suggesting hosts produce hormogonia-inducing factors (Nilsson et al., 2006). In an in vitro study, hormogoniainducing factor (HIF) was not released in a culture medium with excess nitrogen. Thus, HIF production is assumed to be stimulated when nitrogen levels are depleted (Campbell and Meeks, 1989; Adams et al., 2013). Aside from the HIF signal, plant hosts are believed to release various chemoattractants to entice potential symbionts. These were hypothesized to be sugar-based molecules, since simple sugars were proven to be attractants for hormogonia initiation (Nilsson et al., 2006; Adams et al., 2013). Recently, a bioassay was conducted using crude methanolic extracts of Cycas revoluta coralloid roots that transformed symbiotic Nostoc filaments into motile hormogonia. This led to the successful isolation of a hormogonium-inducing factor from cycads characterized as diacylglycerol 1-palmitoyl-2 linoleoyl-sn-glycerol (Hashidoko et al., 2019).

Once a cyanobiont successfully enters its host, the plant partner must stop releasing hormogonia-inducing factors to halt hormogonia formation. The reason is to direct the cyanobiont to the next stage of symbiosis and to start forming heterocysts

for nitrogen fixation to occur – a task that cannot be performed by hormogonia filaments. Therefore, the release of hormogoniarepressing factors (HRF) hinders HIF activities. Cohen and Meeks (1997) observed this phenomenon using aqueous tissue extracts of Anthoceros wherein activation of two genes blocks hormogonium formation and the expressions are only induced by HRF but not HIF.

Besides inducing and repressing factors, high amounts of phenolic substances are also present in the mucilaginous material embedded in the cyanobacterial zone of coralloid roots and surrounding cortical layers (Lobakova et al., 2004). Phenolic substances are rich sources of antimicrobial compounds and secondary metabolites (Obukowicz et al., 1981). Interestingly, phenolics are also known to participate in cell signaling and might play a role in establishing and maintaining stable cycad-cyanobacterial symbioses (Grilli Caiola, 1980; Obukowicz et al., 1981).

## An Overview of Genes Involved in Cycads-Cyanobacteria Symbiosis

Though still limited, sequencing of representative complete genomes of hosts and cyanobionts led to the identification of genes involved in plant host-cyanobacterial symbiosis (Ran et al., 2010; Li et al., 2018). This section aims to list some studied genes and their functions significant in the establishment of symbiosis as well as other genes they interact with (**Table 1**). Although most of the genes cited here were not studied using cycad coralloid roots, a similar, if not exact, genetic mechanism in cycad-cyanobacterial symbiosis seems to be occurring as the genes discussed in this section were commonly found among plant-cyanobacterial symbioses.

As mentioned previously, HRF are involved in halting hormogonium formation after initial migration of the symbiont into its target organ in the host that varies among plant hosts as mentioned in the previous section. Using Anthoceros tissue extracts, hrmU and hrmA genes were discovered that suppresses hormogonia formation and allow differentiation of cells into heterocysts (Cohen and Meeks, 1997). The expression of these two genes appears to be controlled by the host, which hinders further activities attempted by HIF (Cohen and Meeks, 1997; Meeks et al., 1999). Other open reading frames (ORFs) – hrmI, hrmR, hrmK and hrmE – were identified that were similar to a family coding for transcriptional inhibitor proteins (Campbell et al., 2003). It was proposed that external repressing signals can possibly be detected by Nostoc incapacitating Hrm proteins in successfully binding to hrm operons and this affects transcription of genes that induce hormogonium formation (Meeks et al., 1999; Campbell et al., 2003; Adams et al., 2013).

The ability to form heterocysts is a vital feature for a cyanobiont to establish symbiosis with a host. Without

heterocysts, a cyanobiont cannot fix nitrogen for its host. Mutations in hetR and hetF did not allow one strain of Nostoc to differentiate into heterocysts and thus, those genes were determined to be directly responsible for heterocyst formation (Wong and Meeks, 2001). The main activator for the development of heterocysts is the hetR gene, which works simultaneously with hetF gene coding, which gives rise to a protein that enhances subsequent transcription of hetR (Wolk, 2000; Wong and Meeks, 2002; Zhang et al., 2006). The expression of ntcA gene is necessary for the production of proteins that cyanobacteria use to regulate nitrogen and in addition, is involved in activating the transcription of various genes, including the hetR gene (Herrero et al., 2001). Studies showed that a mutant species of Nostoc with a defective ntcA gene failed to infect its host despite retention of hormogonia-forming capacity and therefore might be linked with the activation of other genes required for a successful symbiotic relationship to be established (Vega-Palas et al., 1992; Herrero et al., 2001; Adams et al., 2013). A study by Leganes et al. (1994) in Anabaena showed that an altered hepA gene disrupts proper formation of the cell envelope in both heterocysts and akinetes. The thickened walls, made up of a polysaccharide layer, are essential for nitrogen fixation activities to be concentrated in the heterocyst cells and prevent oxygen diffusion (Tikhonovich and Provorov, 2007). As a consequence of a non-functional hepA gene, normal development of the cell envelope is impossible and affects the ability of the cyanobiont to fix nitrogen in its host (Leganes et al., 1994). Both free-living and symbiotic bacteria have nif genes not restricted to cyanobacteria (Corbin et al., 1982; Fay, 1992). Those genes are responsible for forming nitrogenase complexes that convert unusable atmospheric dinitrogen to useful forms like ammonia and in this process, the nifD and nifK genes encode a dinitrogenase heterotetramer that contains an active site to reduce dinitrogen atoms (Dos Santos et al., 2012). The nifDK genes work together with it redox partner, an iron protein encoded by nifH gene, to complete the structure and function of nitrogenase complex (Jasniewski et al., 2018). Therefore, these genes are important for maintaining mutual associations with a host because defects in the genes may affect the nitrogen-fixing capacity of the symbiont (Spaink, 1998).

Several genes have been identified, but their significance in symbioses is still poorly understood. A sigma factor gene, sigH, which appears to be directed by HIF, might play a role in the ability of symbionts to increase invasion success in a targeted host (Campbell et al., 1998). Likewise, the ctpH gene, which is in close proximity with the sigH gene, is also controlled by HIF in Nostoc(Adams et al., 2013) and has a function in the photosystem II mechanisms in Synechocystis (Anbudurai et al., 1994). This gene is interesting as it may have varying physiological roles in different symbiotic strains. Another gene, tprN, codes for proteins known to be necessary in various functions, such as regulating the cell cycle, suppressing transcription mechanisms and transporting proteins (Lamb et al., 1995). This gene is also essential for heterocyst maturation process (Campbell et al., 1996). In line with this, Meeks et al. (1999) observed that a silenced tprN gene in N. punctiforme caused no phenotypic change but the infection rate doubled compared to the wild-type strain, and its transcription elevated when exposed to both HIF and HRF exudates of A. punctatus. However, the implications and its involvement in the infection process are still unknown.

### LIMITATIONS IN CYCAD-CYANOBACTERIA RESEARCH

A substantial number of studies have focused on the growth and development of coralloid roots and the diversity of the endosymbiont community within the cyanobacterial zone using cultures and 16S identification methods (Costa et al., 1999; Thajuddin et al., 2010; Yamada et al., 2012). The underlying mechanisms, however, have been left unexplored. Probable reasons challenging the study of symbiosis compared to other hosts with a cyanobacterial partner are discussed in this section.

One of the difficulties in this field of research is studying the in vivo development of cycad – cyanobacterial symbiosis. To determine the mutualistic interactions exhibited by both partners, a significant length of time, which may take years of research, might be needed to confirm and validate hypotheses as development in cycads takes much time. Tissue culture methods are also lacking for most cycads to facilitate in vitro experimental assays that could be used to compare wild-type and symbiont-free conditions. Another reason is the availability of coralloid roots of good quantity and quality. According to McLuckie (1922) cited in Milindasuta (1975) and also based on personal observations, even when coralloid roots are present, a green layer in the cortex may be absent (Nilsson et al., 2006). Moreover, the availability of both root tissue samples and cyanobionts per plant host are often insufficient for experiments, especially when replicates are

needed. Thus, samples are usually pooled using a number of representative hosts. Also, most studies obtain their samples from botanic gardens and not from the wild. Although promising results can still be achieved, gathering samples from natural habitats could add more interesting findings.

Another frustration in this field of research is trying to separate root tissues from the endosymbiont as well as the endosymbiont from the sticky mucilaginous material. Various methods have been proposed (Ferris and Hirsch, 1991; Lindblad et al., 1991; Gehringer et al., 2010), but completely separating them from each requires improvement as most requires manually scraping off the green symbiont from the roots. In determining the functions of cycads and their cyanobionts during symbiosis, the gene expression must be analyzed capturing the symbiotic condition while both partners are together. This eliminates the need to separate them because it would cause changes in gene expression, but poses the more difficult task of how to monitor their association in vivo.

#### CONCLUSION AND FUTURE PERSPECTIVES

Research on coralloid roots is not new, but progress has been quite slow compared to research on Rhizobia and other host-cyanobiont associations. Due to the limitations mentioned above, the research field has been confined to studying the diversity of symbionts obtained from various cycad hosts and revalidating previous hypotheses on how both partners benefit from each other. Cruz-Morales et al. (2017) identified novel biosynthetic gene clusters unique to cycad coralloid rootscyanobacteria symbiosis through genome mining, a research area worth undertaking to understand coevolution and to discover pathways responsible for the synthesis of natural products. With the advancements in genomic research that are becoming more affordable, even for small-scale laboratories, studying cycad host-cyanobiont symbiosis may now move forward at a quicker pace.  In line with this, DNA microarray technology can now be used for analysis of expressed genes of microbes from multiple genomes simultaneously using probes to determine up or downregulated genes. Alternatively, RNA sequencing (RNA-Seq) of cyanobiont genomes approximately 5.4–9.0 Mbp – in obligate cyanobiont Nostoc azollae and facultative cyanobiont Nostoc punctiforme PCC 73102 of Azolla filiculoides, respectively (Ran et al., 2010) and 6.7 Mbp in Nostoc cycadae of Cycas revoluta (Kanesaki et al., 2018) – can be

#### REFERENCES

Adams, D. G., Bergman, B., Nierzwicki-Bauer, S. A., Duggan, P. S., Rai, A. N., and Schubler, A. (2013). "Cyanobacterial-plant symbioses," in The Prokaryotes – Prokaryotic Biology and Symbiotic Associations, eds E. Rosenberg, E. F. DeLong, S. Lory, E. Stackebrandt, and F. Thompson (Berlin: Springer-Verlag), 359–400. doi: 10.1007/978-3-642-30194-0\_17

utilized that could detect new genes and splicing events. Aside from ruling out the need for manual separation of microbe from host tissues, Next generation sequencing (NGS) technology provides high resolution data. Likewise, tissue-specific whole transcriptome profiling might be applied to coralloid root tissues with genome size of about 20–30 Gbp (Wang et al., 2018) for determining expressed genes during symbiosis. Furthermore, gaining bioinformatics skills is necessary to maximize the analysis of the outputs obtained from huge amounts of genomic data.

Cyanobacteria are also valuable sources of metabolites (Haque et al., 2017). Thus, venturing into cycad and/or cyanobiont metabolomics research might help to identify useful products, such biofertilizers, antimicrobials and other natural products of economic importance. It may also facilitate studies on cyanotoxins and biomagnification theories. Likewise, metabolites are also believed to play major roles in cell signaling and communication and could also assist in providing insights on symbiotic pathways and identification of enzymes involved in symbiosis. Many areas of cycad-cyanobacteria symbiosis are still waiting to be explored. With technology rapidly advancing simultaneously with the skills of researchers, restrictions previously deemed overwhelming currently appear to be promising.

### AUTHOR CONTRIBUTIONS

TC and AC designed the study. AC gathered the materials and wrote the manuscript. TC, NL, and JD reviewed the manuscript. All authors read and approved the final manuscript.

### FUNDING

This work was supported by the special financial project from the National Forestry and Grassland Administration for Rescues and Breeding of Wild Rare and Endangered Species (NL 2017, 2018, and 2019), the Shenzhen Urban Management Bureau (NL 201411), the Shenzhen Key Laboratory of Southern Subtropical Plant Diversity, Fairy Lake Botanical Garden (TC 2017–2019), and South China Botanical Garden (JD 2017–2019), Chinese Academy of Sciences.

### ACKNOWLEDGMENTS

We would like to thank Dr. David E. Boufford at Harvard University Herbaria for the help while preparing the manuscript.

Adams, D. G., and Duggan, P. (1999). Heterocyst and akinete differentiation in cyanobacteria: Tansley Review No. 107. New Phytol. 144, 3–33. doi: 10.1046/j. 1469-8137.1999.00505.x

and F. Baluška (Berlin: Springer Berlin Heidelberg), 93–121. doi: 10.1007/978- 3-642-20966-6\_5

Chan, C. X., and Bhattacharya, D. (2010). The origin of plastids. Nat. Educ. 3:84.

Haselkorn, R. (2016). Cyanobacteria. Curr. Biol. 19, 277–278.

in Nostoc ellipsosporum. Mol. Microbiol. 12, 679–684. doi: 10.1111/j.1365-2958. 1994.tb01055.x

Phytopathology, eds Y. T. Dyakov, T. Korpela, and V. G. Dzhavakhiya (Amsterdam: Elsevier), 365–420. doi: 10.1016/b978-044452132-3/50018-3

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Chang, Chen, Li and Duan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Non-Legume Parasponia andersonii Mediates the Fitness of Nitrogen-Fixing Rhizobial Symbionts Under High Nitrogen Conditions

#### Edited by:

Kevin Garcia, North Carolina State University, United States

#### Reviewed by:

Sergio Svistoonoff, Institut de Recherche pour le Développement (IRD), France Oswaldo Valdes-Lopez, National Autonomous University of Mexico, Mexico

\*Correspondence:

Simon E. Dupin<EMAIL_ADDRESS>† These authors have contributed equally to this work

#### Specialty section:

This article was submitted to Plant Microbe Interactions, a section of the journal Frontiers in Plant Science

Received: 18 March 2019 Accepted: 20 December 2019 Published: 07 February 2020

#### Citation:

Dupin SE, Geurts R and Kiers ET (2020) The Non-Legume Parasponia andersonii Mediates the Fitness of Nitrogen-Fixing Rhizobial Symbionts Under High Nitrogen Conditions. Front. Plant Sci. 10:1779. doi: 10.3389/fpls.2019.01779 Simon E. Dupin1,2\*, René Geurts <sup>2</sup>† and E. Toby Kiers <sup>1</sup>†

<sup>1</sup> Department of Ecological Science, Vrije Universiteit Amsterdam, Amsterdam, Netherlands, <sup>2</sup> Department of Plant Sciences, Wageningen University, Wageningen, Netherlands

Organisms rely on symbiotic associations for metabolism, protection, and energy. However, these intimate partnerships can be vulnerable to exploitation. What prevents microbial mutualists from parasitizing their hosts? In legumes, there is evidence that hosts have evolved sophisticated mechanisms to manage their symbiotic rhizobia, but the generality and evolutionary origins of these control mechanisms are under debate. Here, we focused on the symbiosis between Parasponia hosts and N2-fixing rhizobium bacteria. Parasponia is the only non-legume lineage to have evolved a rhizobial symbiosis and thus provides an evolutionary replicate to test how rhizobial exploitation is controlled. A key question is whether Parasponia hosts can prevent colonization of rhizobia under high nitrogen conditions, when the contribution of the symbiont becomes nonessential. We grew Parasponia andersonii inoculated with Bradyrhizobium elkanii under four ammonium nitrate concentrations in a controlled growth chamber. We measured shoot and root dry weight, nodule number, nodule fresh weight, nodule volume. To quantify viable rhizobial populations in planta, we crushed nodules and determined colony forming units (CFU), as a rhizobia fitness proxy. We show that, like legumes and actinorhizal plants, P. andersonii is able to control nodule symbiosis in response to exogenous nitrogen. While the relative host growth benefits of inoculation decreased with nitrogen fertilization, our highest ammonium nitrate concentration (3.75 mM) was sufficient to prevent nodule formation on inoculated roots. Rhizobial populations were highest in nitrogen free medium. While we do not yet know the mechanism, our results suggest that control mechanisms over rhizobia are not exclusive to the legume clade.

Keywords: Parasponia, nitrogen fixing bacteria, plant nutrition, host control, rhizobium fitness, nitrogen fertilizer, non-legume, nodulation

**85**

## INTRODUCTION

Symbiotic partnerships have transformed the Earth's nutrient cycles and facilitated rapid adaptation of species to divergent new niches (Joy, 2013; McFall-Ngai et al., 2013; Oldroyd, 2013). Because almost all organisms rely on symbiotic associations for some form of metabolism, protection, or energy (Bronstein, 2015), this immense symbiotic microbial world has been called "the biggest frontier that biology has presented us in a long time" (Carey, 2015).

Despite their importance, understanding the origins and evolutionary trajectories of symbiotic partnerships remains a major challenge. The problem is that mutual benefit does not guarantee evolutionary stability, and partnerships can be vulnerable to exploitation (Sachs and Simms, 2006; Kiers and Denison, 2008; Sachs et al., 2011; Sachs and Hollowell, 2012). What prevents microbial mutualists from defecting from symbiotic cooperation and parasitizing their hosts? This is a question being asked across a diversity of host organisms, from plants and insects to humans (Werner et al., 2014; Hallam and McCutcheon, 2015; Keeling et al., 2015; Werner et al., 2015). While it is appreciated that symbiosis is a key underlying mechanism behind the complexity of life, we do not have a general understanding of how symbiotic associations are controlled and harnessed (Keeling et al., 2015).

 Host plants employ a range of (non-exclusive) strategies to maximize rhizobial benefits, including: 1) pre-nodule control in which the plant evolved high levels of specificity to achieve species, or even strain-specific selection on rhizobial partners, and 2) control based on rhizobial performance, such that higher performing nodules receive proportionally more resources than poor performing nodules (Kiers and Denison, 2008; Westhoek et al., 2017). Because signaling (i.e., pre-nodule control) can be vulnerable to partners that cheat by evolving the correct signal while providing few resources (Edwards and Yu, 2007), it is thought that some form of basal nodule-level control is necessary to prevent rhizobial exploitation (Denison and Kiers, 2011; Regus et al., 2017a; Kiers et al., 2003; Kiers et al., 2006; Oono et al., 2011).

Control of exploitation is particularly important when plant hosts have direct access to high levels of nitrogen in the soil. Under these conditions, the nitrogen benefits provided by the rhizobial symbiont become redundant with nitrogen provided by the environment. Hosts are therefore expected to evolve mechanisms which prevent nodule formation when grown in high nitrogen soils (Regus et al., 2017b). While such regulatory pathways have been documented in legume species (Streeter and Wong, 1988; Cabeza et al., 2014; Soyano et al., 2014; Nishida et al., 2016; Nishida and Suzaki, 2018), the generality of these control mechanisms are unknown (Heath et al., 2010; Regus et al., 2014; Regus et al., 2015).

Our aim was to study host control mechanisms outside the legumes. We focused on the non-legume Parasponia andersonii (Cannabaceae). The genus Parasponia is composed of pioneer plant species typically found on nitrogen poor slopes of volcanic hills in the Malay Archipelago. Parasponia is the only lineage outside the legume family to be able to form a nodule symbiosis with rhizobium (Trinick, 1973; Trinick and Galbraith, 1980; Trinick and Hadobas, 1988). Recent phylogenomic studies suggest that nodule symbioses with diazotrophic bacteria share a single evolutionary origin (Griesmann et al., 2018; van Velzen et al., 2019). As the Parasponia and legumes diverged >100 million years ago, microbial partner selection strategies were shaped independently in both lineages (van Velzen et al., 2019). Therefore, Parasponia is a unique evolutionary replicate of rhizobium symbiosis that allows us to better understand the evolutionary origins of control mechanisms.

In legumes, such Glycine max, past work has shown that the severity of symbiont control is mediated by the addition of external nitrate—as nitrogen availability increases, the host reduces resources allocated to the symbiont. Because soil nitrate is generally less costly for legumes compared to biologically fixed nitrogen, this leads to an inhibition or severe reduction of legume nodule formation (Streeter and Wong, 1988; Voisin et al., 2002; Wendlandt et al., 2018; Regus et al., 2017b). This process, however, is not well understood in non-legumes. In Parasponia, early work has shown nodules can continue to form in high nitrogen environments (Becking, 1983b), but whether rhizobia remain viable in nodules under high nitrogen conditions is unknown. This is important because in Parasponia nodules, rhizobium are not terminally differentiated (Alunni and Gourion, 2016). This means that the endosymbiont population within nodules can replicate, and will be added to the soil population upon nodule senescence. Can non-legumes control symbionts under conditions when the symbionts become a cost rather than a benefit? Studying the existence of control patterns has been particularly challenging in Parasponia nodules due to the difficulty of growing these tropical trees in greenhouses, and the small size of their nodules compared to most model legumes. Typical metrics, such as growth parameters of individual nodules, poly-3-hydroxybutyrate (PHB) content, and rhizobial fitness measures, have been historically difficult to obtain.

Here we study the effects of increasing nitrogen fertilization on the symbiosis between P. andersonii and the rhizobial symbiont Bradyrhizobium elkanii. Our aim was to ask if fertilization: 1) reduced or eliminate the growth benefits of rhizobial nodulation for P. andersonii, and 2) reduced the fitness benefits for the rhizobial symbiont. If Parasponia has evolved effective mechanisms to control nodulation under high nitrogen, then we expect increasing fertilization to be correlated with negative fitness consequences for rhizobia. If these mechanisms are costly for the host to enact, then we expect to see a host growth depression in the presence—but not absence of a rhizobial symbiont under the high nitrogen treatments.

We grew plants under 0, 0.0375, 0.375, 3.75 mM ammonium nitrate concentrations, either with or without rhizobial inoculation. Our four levels were chosen to represent specific ecological challenges for the host-symbiont, namely: i) when Parasponia depends entirely on its symbiont for nitrogen input (0 and 0.0375 mM), ii) when benefits from inoculation are minimal (0.375 mM), iii) when benefits from inoculation are absent/negative (3.75 mM). After 4 weeks, we quantified shoot and root dry weight, nodule number, and nodule fresh weight. We developed an imaging protocol to measure the projected surface of individual nodule areas and then converted it to volume as a second metric for rhizobia benefit. We crushed nodules for measures of colony forming units (CFU) to determine rhizobial populations per nodule. Together, these metrics (nodule number, biomass, volume, and CFUs) provided us with a proxy for in planta rhizobial fitness that allowed us to better understand host control in non-legumes.

#### MATERIALS AND METHODS

#### Seed Germination

We harvested fresh berries from in-vitro propagated Parasponia andersonii trees genotype W1-14 (Van Velzen et al., 2018) maintained in a tropical greenhouse. We extracted the seeds from the berries by soaking them in water and gentle rubbing against a fine sieve. We surface sterilized all seeds in 4% sodium hypochlorite for 15 min and washed seven times with sterile MQ water. We induced germination by temperature cycle (4 h 28°C, 4 h 7°C) for 12 days. We incubated seed on Schenk and Hildebrandt medium agar plates for 10 days until cotyledons were fully emerged.

#### Experimental Design and Plant Growth

We prepared the growing medium of sterile perlite and sterile river sand. We added 210 g of each perlite mixture to 10 sterile polypropylene containers (OS140box, Duchefa Biochemie) allowing for better gas exchange. Per pot, we placed 4 cm<sup>3</sup> of the river sand mixture to transfer the seedling and avoid root desiccation. We used a factorial design experiment consisting of two rhizobia conditions, four nitrogen levels, and 10 replicates per treatment (2 x 4 x 10 = 80 pots).

#### Inoculation and Nitrogen Treatments

To inoculate P. andersonii, we chose the highly efficient nodulating strain B. elkanii WUR3 (Op den Camp et al., 2012). To prepare the inoculum, we grew a WUR3 pre-culture from a single colony in liquid peptone-salts-yeast (PSY) medium at 28°C, 60 rpm (Regensburger and Hennecke, 1983). One milliliter of the pre-culture OD600 = 0.8 was used to inoculate a 200 ml Erlenmeyer culture. We harvested the culture by centrifugation [10 min at 3,500 x relative centrifugal force (rcf)] at OD600 = 0.8. We suspended the cells in the different EKM (Becking, 1983a) solutions to an OD600 = 0.05. The nitrogen treatments were based on an EKM-medium with four levels (0, 0.0375, 0.375, and 3.75 mM) ammonium nitrate. Nitrogen and rhizobia inoculum were added by saturating the perlite and river sand with the four EKM-medium prior to transferring the seedlings to pots. Non-inoculated controls received EKM solutions (see below) but without the rhizobia culture. We then randomly placed the pots on a growth chamber table under a 16/8 h light cycle, temperature 26/24°C, light intensity 185 µmol/m2/s, and a relative humidity of 90%.

#### Harvest

We harvested plants after 30 days. We carefully washed off the perlite and sand from the root systems. We counted nodules and harvested each one individually. We then used binoculars equipped with a Nikon camera (DS–Fi2) to image each nodule. To obtain nodule volume, we extracted and measured the area and perimeter of the nodules photographed using FIJI (Schindelin et al., 2012). We calculated the corresponding prolate spheroid volume using the best fitted ellipse of each nodule based on a previously developed formula (Nedomová et al., 2014). We weighed nodules and kept them at 4°C in 0.9% NaCl solution until they could be crushed for fitness assays. We separated shoots and roots, dried them 72 h at 60°C, and weighed them.

#### Metrics for Rhizobial Fitness Proxies

To determine a fitness proxy for rhizobia per plant, we surface sterilized all nodules with 96% ethanol for 20 s, 4% sodium hypochlorite for 1 min, and washed seven times with sterile water. We crushed nodules in 150 µl 0.9% sterile saline solution. Fifty microliters of the crushed nodules was diluted in series and both 10,000 and 100,000 dilutions were streaked on PSY plates with sterile glass beads and incubated at 28°C for 7 days. We then counted colonies to determine total rhizobia per plant.

#### Statistical Analysis

We used R version 3.6.0 (2019-04-26) to conduct all statistical tests. In case of heteroscedasticity or non-normality, a decimal logarithm transformation of the data was performed to meet ANOVA assumptions. If ANOVA assumption could not be met with a transformation of the response variables, a nonparametric Kruskal test was conducted. To compare plant biomass among treatments, we tested the decimal logarithm mean plant dry biomass for significant differences with two-way ANOVA and a post-hoc Tukey tests for pairwise comparison with a 95% confidence interval. To test for differences in allocation to above and belowground parts, we compared mean root to shoot ratio with a pair-wise Wilcoxon test. To compare plant biomass as a function of rhizobial inoculation, we compared decimal logarithm relative plant biomass among nitrogen treatments with a one way ANOVA and a Tukey test for multiple group comparison with a 95% confidence interval. To test for differences in nodule formation, we compared mean nodule number and nodule fresh per plant with one-way ANOVA and a Tukey test for multiple group comparison with a 95% confidence interval. For nodule volume we used the Kruskal test and Dunn's post-hoc. Rhizobial fitness components, defined as viable rhizobia per milligram of plant, nodule mass, or volume, were compared among the four nitrogen treatments of inoculated plants with Kruskal test and Dunn's post-hoc. Dunn's test is appropriate for groups with unequal numbers of observations (Zar, 2010), and was corrected for multiple comparisons following Benjamini and Hochberg method (1995) with a 95% confidence interval.

### RESULTS

#### Plant Biomass

We first asked how increasing soil nitrogen affected plant growth patterns in the presence and absence of rhizobial symbionts. We found that nitrogen treatment and rhizobial treatment both had a significant effect on total plant biomass (ANOVA; F3,72 = 330.6, p < 0.001 and F1,72 = 112.5, p < 0.001 respectively), with a significant interaction term between the two variables (ANOVA; F3,72 = 40.2, p < 0.001, Figure 1A). Specifically, we found the presence of rhizobia increased total plant biomass at 0, 0.0375, and 0.375 mM ammonium nitrate levels. At the highest nitrate level (3.75 mM), we found that the presence of rhizobia became a cost (Figure 1A). By quantifying the root to shoot ratio across these nitrogen levels, we found that at the highest nitrate level, the rhizobial cost was related to a decrease in root biomass (Figure 1B).

To assess plant biomass relative to rhizobia inoculation under different nitrogen regime, we calculated the relative growth benefit of inoculation for the plant (Figure 1C). Here, we took the difference of the total dry weight of the size-matched inoculated and non-inoculated control, divided by the total dry weight of non-inoculated control. We found that the relative growth benefits of rhizobial inoculation decreased with nitrogen fertilization (ANOVA; F3,36 = 2,685, p < 0.001).

## Formation of Symbiotic Organs

We then asked how nitrogen fertilization affected the formation of symbiotic organs, namely nodule number, total nodule fresh weight, and nodule volume. We found a significant effect of fertilization on all three parameters: nodule number (one-way ANOVA; F3,35 = 16.5, p < 0.001) nodule weight (one-way ANOVA; F3,36 = 41.08, p < 0.001), and nodule volume (Kruskal-Wallis test; chi-squared = 23.07, df = 3, p-value < 0.001) (Figure 2). This effect was driven largely by the highest nitrogen concentration. In the three lowest nitrogen concentrations, we found no difference in nodule number (~4– 5 per plant, Figure 2A), fresh weight (~6 mg, Figure 2B), or volume (~70 mm3 , Figure 2C). However, when fertilization was increased to 3.75 mM NH4NO3, all nodule parameters were reduced to nearly zero, across all replicates. This demonstrates

FIGURE 2 | Inoculated Parasponia andersonii nodulation under four nitrogen levels. (A) Mean nodule number per plant. (B) Mean nodule fresh weight per plant. (C) Mean nodule volume per plant. Letters represent groups significantly different from each other (ANOVA, Tukey HSD test and Kruskal test, Dunn's test). Error bars show standard error.

that with enough exogenous nitrogen available, P. andersonii is able to prevent nodule organogenesis, similar as reported for legumes and actinorhizal plants.

In legumes, it has been demonstrated that the rhizobia—as well as exogenous fixed-nitrogen—can trigger expression of CLE peptide encoding genes, which can trigger systemic signalling and inhibit nodulation (Okamoto et al., 2009; Reid et al., 2011a). Studies in Lotus japonicusrevealed that whereas LjCLE-RS1 and LjCLE-RS2 are induced upon rhizobium inoculation, only LjCLE-RS2 is induced by application of exogenous nitrogen (Okamoto et al., 2009). The latter gene is a close homolog of the soybean (G. max) CLE peptide encoding geneGmNIC1a that is induced by exogenous nitrate (Reid et al., 2011b; Hastwell et al., 2017). For P. andersonii, CLE genes have been annotated (Table S1). Of these, PanCLE5 and PanCLE9 showed to be close homologs of LjCLE-RS1 and LjCLE-RS2, and their counterparts MtCLE12 and MtCLE13 in Medicago truncatula (Van Velzen et al., 2018). To obtain insight whether any of the P. andersonii CLE genes is induced by exogenous nitrate, we exploited available RNA sequencing (RNA-seq) data (Van Velzen et al., 2018). This revealed that three CLE genes, namely PanCLE2, PanCLE8, and PanCLE9, have an increased expression in inoculated roots grown at relatively high exogenous nitrate levels (5 mM KNO3) (Figure S1). These three genes showed a similar induction in young, non-infected, nodule primordia (Figure S1). This suggests that in P. andersoniirhizobium and exogenous nitrate trigger an overlapping CLE gene repertoire to regulate nodulation.

#### Rhizobial Fitness Proxies

While nodule number and weight can give some rough estimates of rhizobial benefit, a key parameter is to directly quantify rhizobial densities in nodules. We therefore next measured CFUs per mg of plant and nodule biomass or volume (n = 2 for 3.75 mM as only two plants formed nodules). We found that as exogenous nitrogen levels increased, there was a concurrent decrease in CFUs per plant biomass (Kruskal-Wallis test; chisquared = 9.1, df = 3, p-value = 0.028), and CFUs per mg of nodule (Kruskal-Wallis test; chi-squared = 8.2, df = 3, p-value = 0.041). But surprisingly, no significance difference was found for CFUs per nodule volume (Kruskal-Wallis test; chi-squared = 4.9, df = 3, p-value = 0.18). Despite the fact that nodule number and fresh weight remained the same across the lowest three nitrogen levels, we documented a decrease in viable rhizobial density in nodules as nitrogen levels increased (Figure 3).

#### DISCUSSION

Plant species found in the nitrogen fixing clade rely on diazotrophic symbionts to acquire nitrogen in poor soils (McKey, 1994; Crews, 1999; Vitousek et al., 2002). The ability to control resource allocation to these symbionts is likely a key requirement for evolutionary maintenance of the symbiosis (Werner et al., 2015). We found evidence that the non-legume P. andersonii has likewise evolved mechanisms to control rhizobial fitness, despite an independent evolutionary trajectory from the legume lineage of >100 million years ago. Specifically, we found that increasing

forming units per mg plant dry weight. (B) Mean colony forming units per milligram nodule fresh weight. (C) Mean colony forming units per volume of nodule. Letters represent groups significantly different from each other (Kruskal test, Dunn's test). Error bars show standard error.

nitrogen levels lead to a decrease in rhizobial populations within nodules. In contrast, evidence for control of nodule formation (i.e., organogenesis) was only evident from the highest nitrogen level, in which nodule formation was almost completely suppressed. In the lower nitrogen treatments, we found no evidence for differences in nodule size, number, or volume. Together, these data suggest that despite the continued growth and formation of nodules at lower nitrogen levels, the host still can control the success of the rhizobial populations within those nodules.

An open question is whether Parasponia nodulation is somehow less advanced than in legumes (Behm et al., 2014): our results point to similar levels of nodulation control as found in some legume species. Specifically, it has been shown that at high nitrogen levels, legumes can control nodule formation and mass—a process known as autoregulation of nodulation (Cho and Harper, 1991; Sueyoshi et al., 2003; Jeudy et al., 2010). This has also been demonstrated in nodulated actinorhizal plants, which likewise show a reduced investment in nodulation with application of exogenous nitrogen (Kohls and Baker, 1989; Thomas and Berry, 1989; Arnone et al., 1994; Markham and Zekveld, 2007; Wall and Berry, 2007). Our findings are therefore in line with the idea that a shared mechanism to control nodulation investment across the nitrogen fixing clade was in place before the evolutionary diversion of the nodulating lineages.

How does such control operate? New work has shown that legumes employ a nitrate response to inhibit rhizobial symbiosis by upregulating specific transcription factors (Nishida et al., 2018). Autoregulation of nodulation works as a negativefeedbackfrom the key transcriptionfactor NIN targeting CLE peptides, which induces a shoot response, production of cytokinin and inhibition of nodulation (Nishida and Suzaki, 2018). The nitrate-induced inhibition induces the NIN-like transcription factor NRSYM1 targeting the same CLE peptides. Furthermore, in the actinorhizal plantCasuarina glauca, RNA interference (RNAi) NIN knockdown also showed the essential role of NIN in controlling nodule formation (Clavijo et al., 2015). While it is unknown if Parasponia employs these same mechanisms, CLE peptides are expressed in P. andersonii nodules (Figure S1), suggesting it could be the case (Van Velzen et al., 2018).

However, the ability to control resource allocation to nodules may be less important to Parasponia hosts, given the extremely nutrient poor soils in which they are typically found (Achmad and Hadi, 2016). In this ecological niche, we would not expect hosts to be exposed to high nitrogen conditions, nor the concurrent selection pressures against nodulation. There is also evidence that host responses to nitrate may evolve differently across different plant lineages. Studies of the legume Acmispon strigosus found increases in nodule number and size at low nitrate levels, as expected, but also revealed that nodulation suppression was linked with high plant mortality suggesting a high, direct sensitivity to nitrate (Regus et al., 2017b). In L. japonicus, fertilization reduced nodule size and nodule number, but with no apparent cost on plant fitness (Nishida et al., 2018). In actinorhizal plant lineages closely related to P. andersonii, added nitrate reduces and blocks nodule formation of the Frankia symbiosis. In a split root experiment, Casuarina cunninghamiana showed localized control depending on exogenous nitrate concentration (Kohls and Baker, 1989; Arnone et al., 1994). In Parasponia, nodule size and number did not vary with nitrogen concentration. Instead, nodulation was nearly eliminated for the highest level of fertilization.

Our data show that P. andersonii can block nodulation, but that this might entail a cost. We found that at the highest nitrogen level, plant root biomass was reduced in the presence of rhizobia compared to non-inoculated controls. But this growth depression was not observed at lower nitrogen concentrations. Moreover, we observed a linear decrease in the benefit to plant growth conferred by rhizobia as nitrogen increased. Specifically, our data show that resource allocation to root growth was negatively impacted under high N when rhizobia were present (Figure 1B). This cost in high nitrogen context was not linked with allocation to resources to nodules organogenesis because nodule formation was suppressed. Instead, this may be linked to the presence of intercellular bacteria, nutrient uptake, or plant defense mechanisms (Bakker et al., 2018). A similar result was also found in A. strigosus, whereby the authors noted a cost to the presence of rhizobia under high nitrogen, even though the cost was also not linked to nodule formation (Regus et al., 2017b). They point to past work showing the cost of chemically induced plants defense response to pathogens (Heil et al., 2000). Whether our reduction in root growth (Figure 1B) is linked to the costs of upregulating plant defense response is unknown.

Although our data suggest that Parasponia employs a mechanism to control nodulation upon presence of exogenous nitrogen—and that this is likely linked to CLE signaling—we also find that rhizobial density may be likewise regulated within nodules. Specifically, we found that P. andersonii controlled rhizobial colonization levels at a per nodule milligram and per plant milligram basis. While plants had similar numbers of nodules and equivalent nodule weights across three nitrogen levels, we found that the amount of viable rhizobia hosted in host cells varied. Control of rhizobial fitness within nodules has been shown in legumes, most recently in L. japonicus, in which there is evidence that plants can differentially control fitness of effective and ineffective rhizobia within a single nodule (Quides et al., 2017). Similarly, sanction strength—meaning the ability to control rhizobial fitness in individual nodules—was also predicted (West et al., 2002) and shown (Kiers et al., 2006) to increase with addition of external nitrate in soybeans.

While our data suggest that within nodule regulation of rhizobial fitness depends on nitrogen levels, a broader question is whether Parasponia has evolved a similar response to ineffective rhizobia that fail to provide nitrogen. Increasingly, work has shown that allocation to nodules will depend on the quality of the rhizobial partners, such that nodules containing low-quality partners will be sanctioned, and experience a reduction in resources (Regus et al., 2017a; Kiers et al., 2003; Oono et al., 2009; Oono et al., 2011). For example work in Lupinus arboreus has shown that the size of the nodule is linked with the quality of its occupant (Simms et al., 2006). While we used a well-characterized effective strain, future work should aim to understand how rhizobial partner quality and nitrogen levels interact in the Parasponia-rhizobia symbiosis, studied in some legumes (Grillo et al., 2016; Regus et al., 2017a; Regus et al., 2014; Wendlandt et al., 2018).

A second open question is how the physiology of Parasponia nodules affects the potential for hosts to control rhizobial fitness. P. andersonii has indeterminate nodule with a central vasculature, meaning a meristem sustains a continuous growth of the nodule such that cells become colonized with infection threads containing rhizobia (Behm et al., 2014). When the cells are fully colonized, the rhizobia cells are kept in fixation threads and do not differentiate to bacteroids (i.e., swollen bacteria unable of cell division) (Trinick and Galbraith, 1976). Because of this mode of growth, we had expected P. andersonii to reduce nodule growth (meristematic cell division) upon fertilization, yet we observed similar nodule size with an overall lower cell colonization by rhizobia. While the mechanisms is still unknown, this result suggests that P. andersonii can directly reduce rhizobia cell division within its nodules or reduce rhizobia nodule occupancy by inducing cell senescence.

More generally, future work is needed to better characterize the costs and benefits of the symbiosis physiologically. For example, from the host side, measurements of %Ndfa (nitrogen derived from the atmosphere) can help us more accurately understand the contribution of nitrogen fixation under different fertilizer regimes. Likewise, detailed microscopy of P. andersonii nodules could be conducted. Here, electron microscopy would be useful to study the integrity of the fixation threads and endosymbiotic bacteria, whereas light microscopy on replicate nodules could help develop reliable quantifying techniques based on visual inspection [as in (Regus et al., 2017a)]. From the symbiont side, quantification of metrics such as (PHB), could be helpful in understanding rhizobial fitness, specifically how PHB is linked to reproduction and survival during starvation [e.g., (Ratcliff et al., 2008)].

Overall our work suggests that rhizobial controlmechanisms are not exclusive to legumes. While there was evidence that the relative host growth benefits of inoculation decreased with nitrogen fertilization, we found that Parasponia controls rhizobial fitness, likely by mediating rhizobial density, depending on ammonium nitrate availability. A key open question is how these processes operate in the field, where Parasponia evolved on nutrient poor volcanic soils. Given increasing global nutrient inputs, even to pristine ecosystems, more data are needed to understand fitness alignment in the Parasponia—rhizobium symbiosis under changing nutrient conditions.

#### DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to the corresponding author.

### AUTHOR CONTRIBUTIONS

All authors contributed to frame the research question and design of the experiments. SD performed the experimental work. All authors contributed to write the manuscript and participated in the discussion.

#### FUNDING

This work was supported by NWO open competition grant to ETK (819.01.007).

### ACKNOWLEDGMENTS

We thank Victor Caldas for the help with developing the nodule image analysis, and the three referees for their useful feedback and suggestions.

#### REFERENCES

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019. 01779/full#supplementary-material

FIGURE S1 | Expression profile of P. andersonii CLE peptide encoding gene. Expression profile of P. andersonii CLE genes in non inoculated roots, nodule promordia and inoculated roots under low (0.5 mM KNO3) and high nitrate (5.0 mM KNO3) conditions. Expression is given in DESeq2-normalized read counts; error bars represent SE of three biological replicates. Dots represent individual expression levels.

Trinick, M. J. (1973). Symbiosis between Rhizobium and the non-legume, Trema aspera. Nature 244, 459–460. doi: 10.1038/244459a0

Conflict of Interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Dupin, Geurts and Kiers. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Aphids Influence Soil Fungal Communities in Conventional Agricultural Systems

Thomas D. J. Wilkinson<sup>1</sup>† , Jean-Pascal Miranda<sup>1</sup> , Julia Ferrari<sup>1</sup> \*, Sue E. Hartley1,2 and Angela Hodge<sup>1</sup>

<sup>1</sup> Department of Biology, University of York, York, United Kingdom, <sup>2</sup> York Environmental Sustainability Institute, University of York, York, United Kingdom

Arbuscular mycorrhizal fungi (AMF) form symbioses with the roots of most plant species, including cereals. AMF can increase the uptake of nutrients including nitrogen (N) and phosphorus (P), and of silicon (Si) as well as increase host resistance to various stresses. Plants can simultaneously interact with above-ground insect herbivores such as aphids, which can alter the proportion of plant roots colonized by AMF. However, it is unknown if aphids impact the structure of AMF communities colonizing plants or the extent of the extraradical mycelium produced in the soil, both of which can influence the defensive and nutritional benefit a plant derives from the symbiosis. This study investigated the effect of aphids on the plant-AMF interaction in a conventionally managed agricultural system. As plants also interact with other soil fungi, the non-AMF fungal community was also investigated. We hypothesized that aphids would depress plant growth, and reduce intraradical AMF colonization, soil fungal hyphal density and the diversity of AM and non-AM fungal communities. To test the effects of aphids, field plots of barley enclosed with insect proof cages were inoculated with Sitobion avenae or remained uninoculated. AMF specific and total fungal amplicon sequencing assessed root fungal communities 46 days after aphid addition. Aphids did not impact above-ground plant biomass, but did increase the grain N:P ratio. Whilst aphid presence had no impact on AMF intraradical colonization, soil fungal hyphal length density, or AMF community characteristics, there was a trend for the aphid treatment to increase vesicle numbers and the relative abundance of the AMF family Gigasporaceae. Contrary to expectations, the aphid treatment also increased the evenness of the total fungal community. This suggests that aphids can influence soil communities in conventional arable systems, a result that could have implications for multitrophic feedback loops between crop pests and soil organisms across the above-below-ground interface.

Keywords: arbuscular mycorrhizal fungi, community diversity, amplicon sequencing, Hordeum vulgare, Sitobion avenae, rhizosphere, multitrophic interactions

### INTRODUCTION

Arbuscular mycorrhizal fungi (AMF) form obligate symbioses with the roots of c. twothirds of land plant species, including agriculturally important cereals (Smith and Read, 2008; Fitter et al., 2011). Enhancing this symbiosis has been proposed as an important tool for increasing food security and agricultural sustainability (Gosling et al., 2006;

#### Edited by:

Heike Bücking, South Dakota State University, United States

#### Reviewed by:

Jan Jansa, Institute of Microbiology (ASCR), Czechia Sergio Saia, Council for Agricultural Research and Economics, Italy

#### \*Correspondence:

Julia Ferrari<EMAIL_ADDRESS>orcid.org/0000-0001-6519-4254

#### †Present address:

Thomas D. J. Wilkinson, ADAS Gleadthorpe, Mansfield, United Kingdom

#### Specialty section:

This article was submitted to Plant Microbe Interactions, a section of the journal Frontiers in Plant Science

Received: 05 March 2019 Accepted: 24 June 2019 Published: 12 July 2019

#### Citation:

Wilkinson TDJ, Miranda J-P, Ferrari J, Hartley SE and Hodge A (2019) Aphids Influence Soil Fungal Communities in Conventional Agricultural Systems. Front. Plant Sci. 10:895. doi: 10.3389/fpls.2019.00895

**95**

Fitter et al., 2011; Jacott et al., 2017; Thirkell et al., 2017). Whilst the host plant provides a fixed carbon (C) source for AMF, AMF transfer nutrients such as nitrogen (N) and phosphorus (P) to the plant (Hodge et al., 2001; Smith et al., 2009; Hodge and Fitter, 2010; Karasawa et al., 2012). AMF colonization affects multitrophic interactions between above- and below-ground herbivores (Yang et al., 2014) and may also enhance the uptake of silicon (Si), which can alleviate the impact of both biotic and abiotic stress (Dias et al., 2014; Garg and Bhandari, 2016; Frew et al., 2017).

The bottom-up effect of below-ground AMF on the performance of above-ground herbivores such as aphids can range from positive to negative (Gange and West, 1994; Wurst et al., 2004; Ueda et al., 2013; Simon et al., 2017; Wilkinson et al., 2019). These impacts on aphid performance likely occur because of alterations to plant defense and nutrition due to the AMF symbiosis (Wurst et al., 2004; Meir and Hunter, 2018b) and can depend on the level of AMF colonization of the host plant (Tomczak and Müller, 2017; Maurya et al., 2018; Meir and Hunter, 2018a). In turn, aphids may impose top-down effects on AMF colonization via the host plant (Babikova et al., 2014; Meir and Hunter, 2018a). Top-down and bottom-up effects can therefore modulate the outcome of each other, potentially resulting in above-below-ground multitrophic feedback loops (Meir and Hunter, 2018a). Thus, if aphids influence AMF colonization this could impact how AMF affect plant nutrient uptake and tolerance to abiotic stress in multitrophic systems. The AMF extraradical mycelium (ERM) phase is of also of key importance for interactions between plants and other rhizosphere organisms (Perotto and Bonfante, 1997; Jones et al., 2004; Hodge and Fitter, 2013) and can be directly related to AMF derived plant nutrient acquisition (Hodge et al., 2001; Barrett et al., 2011). Additionally, the ERM can be involved in plant defense, carrying signals of aphid attack to neighboring plants connected via ERM underground networks (Babikova et al., 2013). Elucidating how intra- and extraradical structures of AMF respond to top down effects is therefore important in understanding their potential for use in complex agro-ecosystems. However, current knowledge of how AMF respond to aphids sharing the same host plant is limited to the impact on AMF colonization (Babikova et al., 2014; Vannette and Hunter, 2014; Maurya et al., 2018; Meir and Hunter, 2018a).

The identity of the taxa within the AMF community colonizing the host plant can be important in determining the nutrient uptake or defense benefit gained from the symbiosis. In small, artificially selected AMF communities, AMF species identity determines the level of protection the AMF provides for the host plant against biotic stressors (Pozo et al., 2002; Sikes et al., 2009; Malik et al., 2016), and certain AMF species may deliver more or less nutrients to their host plant (Jansa et al., 2008; Leigh et al., 2009; Thirkell et al., 2016). Similarly, soil community transfer experiments suggests the AMF community structure can also be important in determining nutrient acquisition and plant growth responses (Hodge and Fitter, 2013; Williams et al., 2014; Manoharan et al., 2017; Jiang et al., 2018). Large vertebrate grazing can affect AMF communities (Ba et al., 2012; Guo et al., 2016), and insect herbivory can alter below-ground ectomycorrhizal (Gehring and Bennett, 2009), non-mycorrhizal fungi (Kostenko et al., 2012) and rhizosphere bacterial community characteristics (Kong et al., 2016), but the impact of arthropod herbivory on AMF communities is currently unknown.

The impact of herbivory on AMF structures is variable (Barto and Rillig, 2010), and phloem feeding aphids can increase or decrease the intraradical AMF colonization of their host plant (Babikova et al., 2014; Meir and Hunter, 2018a). The C limitation hypothesis proposes that aboveground removal of fixed C by herbivory will result in less C available below-ground (Wallace, 1987), although the subsequent allocation of this limited amount of C between roots and AMF contained within roots is unknown. The reduced C availability might result in changes to the AMF community because only a limited number of AMF species can be supported and thus the number and relative abundance of less competitive species might be reduced (Gange, 2007; Ba et al., 2012). Alternatively, low levels of herbivory could lead to more C allocated below-ground in an attempt for the plant to take up more nutrients for regrowth (Wamberg et al., 2003), which could increase fungal diversity (Ba et al., 2012). There are also examples of herbivory affecting the composition of ectomycorrhizal communities rather than species richness, and thus altering the beta diversity of communities, making communities more distinct (Gehring and Bennett, 2009).

Arbuscular mycorrhizal fungi communities in conventionally managed agricultural systems are often distinct from AMF communities in other settings and often have low diversity due to tillage, chemical fertilizer, pesticide and fungicide regimes (Jansa et al., 2002; Gosling et al., 2006; Wetzel et al., 2014; Hartmann et al., 2015; Manoharan et al., 2017). We selected a conventionally managed system to investigate agriculturally relevant AMF communities that are tolerant to such practices. Apart from AMF, other fungi also associate with plant roots, including other endophytic mutualists (Murphy et al., 2015; Lugtenberg et al., 2016). As some of these other fungi and their community composition can influence aphid performance (Hartley and Gange, 2009; Battaglia et al., 2013; Kos et al., 2015), the AMF community must be placed in the context of any changes in the wider fungal community.

Here, we investigate the impact of aphids on the belowground fungal community with a focus on AMF given their key role as ecosystem engineers. Specifically, we tested the following hypotheses: (1) As aphids will depress plant growth and nutrient status, the AMF will benefit less from the association with the plant, and consequently AMF structures, both internal and external to the root, will be reduced. (2) Aphids will cause a reduction in both the alpha diversity and evenness of the soil communities, which results in a distinct soil fungal community composition (increased beta diversity).

#### MATERIALS AND METHODS

#### Site Selection

fpls-10-00895 July 11, 2019 Time: 17:37 # 3

Spring barley (Hordeum vulgare L., cultivar: Planet) was drilled into silty clay loam over chalk (Towthorpe, North Yorkshire, SE 91086 62387, GPS: 54.049412, −0.610267, elevation above sea level: 100 m) on March 15th, 2017. Average soil characteristics of the whole field were sampled on May 12th, 2017 and analyzed by NRM Laboratories (Berkshire, United Kingdom; **Table 1**). The field had been conventionally cropped with wheat and oats for the previous 5 years and was treated with conventional agrochemical inputs throughout the duration of this study (**Table 1**).

#### Aphid Treatments

On the April 21st, 2017, lidless and bottomless PVC boxes (40 cm × 40 cm × 25 cm) were inserted 2–3 cm below the surface of the soil around sections of developing three leaf stage seedlings, averaging 26 ± 1.7 SE plants per box. This shallow insertion, so as not to disturb plant roots, allowed an aphid impermeable seal to form between the soil and a cage structure. The interior of each box thus formed an experimental "plot." Experimental plots were assigned to "+Aphid" or "-Aphid" treatments and arranged in a randomized block design. As the site lay on a North Western slope, plots were set out in two rows perpendicular to the slope, in a North East direction. The location of the plots in the North East and North West direction were coded as the NE and NW coordinate (respectively) of each plot within the field site. This attempted to account for any locational environmental gradients within the site, and are referred to as the NE and NW plot location hereafter (**Figure 1**).

One week later (April 28th, 2017) cages were constructed to cover all plots. The frames of the cage consisted of wooden posts inserted 20 cm into the soil and were attached to the interior of the PVC box. The cage extended 90 cm above the soil and was covered with polypropylene horticultural fleece (**Figure 2**), which intercepted c. 14% of the photosynthetically active radiation. English grain aphids (Sitobion avenae) (a single genotype, originally supplied by Koppert, Holland) were cultured on barley plants (cultivar Quench) at 20◦C. From these cultures, ten 4th instar adults were taken at random and added to each +Aphid plot. All experimental plots (including −Aphid) were sealed with cages. S. avenae populations usually peak in the late summer months (Blackman and Eastop, 2000), however, aphids and cages were added to the crop in the current study earlier than this to stop the natural ingress of aphids into the plots. Initially, eight replicates of each treatment were set up, although one +Aphid replicate was discarded during the study due to damage to the cage caused by high wind speeds and so N = 8 for −Aphid treatments and N = 7 for +Aphid treatments. Cages remained over all treatments for the duration of the experiment, and may have reduced direct contact of some of the agrochemical inputs (see **Table 1**) with the plots after this date.

#### Harvest

Plots were harvested on the June 13th, 2017, 46 days post-aphid addition. The number of aphids on five tillers chosen at random within each plot were recorded, and plots were dug out to 20 cm soil depth before storage at 4◦C overnight. Three soil cores at each plot location were taken with a 2 cm diameter cheese corer between 20 and 30 cm soil depth. These technical replicates were processed separately to assess for the fungal hyphal length density (HLD) from each core (Hodge, 2001) but then the resulting data was pooled to create a biological replicate. After storage, the aphids were washed from any above-ground material and plant roots were separated from the soil and washed whilst still attached to the above-ground biomass. Only those roots visibly attached to a plant were stored at −20◦C for DNA extraction or in 40% ethanol for the staining of fungal structures. The above-ground plant material was oven dried at 70◦C for at least 96 h and the total number of plants, tillers and fertile tillers were recorded. The above-ground plant material for each plot was separated into a combined stem and leaf fraction (henceforth referred to as stem material) as well as a separate grain fraction. It should be noted that due to the size of the stem and leaf material, the stem and leaf fraction was weighed to the nearest gram. These fractions were homogenized in a kitchen blender (Igenix ig 8330, Ipswich, United Kingdom) before ball milling to a fine powder (Retsch MM400, Retsch GmbH, Haan, Germany). The resulting material was then analyzed for C and N ratios via a Elementar Vario El Cube (Elementar UK, Ltd., Stockport, United Kingdom) and pelleted for X-ray fluorescence (XRF) analysis (Thermo Fisher ScientificTM portable X-ray fluorescence analyzer) to determine P and Si concentrations as described by Reidinger et al. (2012).

The roots stored in 40% ethanol were stained for fungal structures via the acetic acid–ink staining method (Vierheilig et al., 1998), modified as in Wilkinson et al. (2019). Structures were assessed under a Nikon 50i eclipse microscope (Nikon UK Ltd., Surrey, United Kingdom) under 200 X magnification. As fungi other than AMF colonize plant roots in natural systems, a method that calculates both the most conservative estimation of AMF root length colonization (RLC) (RLC min) and least conservative (RLC max) was employed (Brundrett et al., 1994).

#### DNA Extraction, PCRs, and Sequencing

Frozen root material was lyophilized for 36 h and ball milled to a fine powder. The DNA was extracted from this material using a DNeasy PowerPlant Pro Kit (QIAGEN N.V, Venlo, Netherlands) according to the manufacturer's instructions with the exception that, in order to increase the DNA yield, the DNA solution was eluted twice through the membrane in the final step. DNA concentrations were assessed (NanoDropT<sup>M</sup> 8000 Spectrophotometer Thermo Fisher Scientific) and diluted to 20 ng µl <sup>−</sup><sup>1</sup> before the PCR analysis.

Two regions (amplicons) of fungal ribosomal DNA were amplified via PCR (**Table 2**); an amplicon that captures the diversity of the entire fungal community (total fungi), and an amplicon that captures AMF specific diversity at a higher resolution and species coverage (AMF specific). After initial amplification via primary PCRs, secondary PCRs attached illumina sequencing barcodes. For all PCRs the reaction consisted of 0.5 µl DNA, 0.1 µl of forward and reverse primers (20 mM) and 12.5 µl BioMix Red (Bioline, London, United Kingdom) made up to 25 µl reaction volume with molecular grade dH2O.

TABLE 1 | Fieldsite soil chemical analyses and agrochemical inputs used throughout the study.

PCRs were carried out using a T100TM Thermal Cycler (BioRAD, Hercules, CA, United States). The PCR products of the secondary PCRs were purified using a QIAquick PCR Purification Kit (Qiagen) and the purified concentrations were measured using a Qubit <sup>R</sup> 3.0 fluorometer (Thermo Fisher ScientificTM, Waltham, MA, United States). To mix the products of the two amplicons at equimolar concentrations to reduce sequencing depth bias during simultaneous DNA sequencing, these products were lyophilized overnight and re-suspended in molecular grade H2O to achieve desired concentrations. Quality control and library preparation was carried out by the University of York Bioscience Technology Facility, and the resulting samples were sequenced using an illumina MiSeq system (illumina, San Diego, CA, United States) at 2 ∼ 300 bp: Briefly, unique barcode sequences (Nextera XT index primers, illumina) were added onto amplicons tagged with illumina adapter sequences via

PCR. Amplicons were then purified and pooled at equimolar ratios and then diluted and denatured. Samples were spiked with PhiX library spikes (illumina) for added sequence variety to enhance the distinguishing of fluorescent signals of clusters during sequencing. Samples were run using the MiSeq 600 cycle kit (illumina).

#### Bioinformatic Analysis

The raw forward and reverse reads were merged together resulting in a total of c. 1.5 million reads which were processed using QIIME2<sup>1</sup> (Caporaso et al., 2010). Reads were stripped of their primer and barcoding sequences and untrimmed reads were discarded (0.9%). Read quality of the merged reads were estimated using eestats2 and reads from AMF specific amplicons were truncated to 270 bp from the front end due to estimated sequence quality drop off at this point, whilst total fungal amplicons were not truncated. Reads were then dereplicated and clustered into operational taxonomic units (OTUs) using Usearch10 (Edgar, 2010) with 97% similarity. The resulting 517 total fungi amplicon OTUs were BLASTed against the UNITE ITS (Kõljalg et al., 2013) database and eight non-fungal OTUs were removed, resulting in 509 total fungi OTUs, with the total read number per sample ranging from 53,333 to 129,338. Samples were normalized to 70,000 reads per sample using Usearch10's "norm" function. Identical virtual taxon accessions according to the UNITE database were merged together yielding the 155 total fungi amplicon OTUs used in subsequent analysis.

A total of 27 OTUs were identified for the AMF specific amplicon. These were BLASTed against the maarjAM database (Öpik et al., 2010) accessed on 13th June 2018 and five OTUs with less than 96% coverage or similarity to taxa in the database were discarded. This resulted in a total read number from 12,035 to 21,706 per sample. Reads per sample were normalized to 16,500. AMF virtual taxonomic (VT) identities were assigned to the OTUs according to the greatest BLAST coverage and similarity; where OTUs could not be assigned to a single VT, VTs were labeled as unassigned. A phylogenetic tree of the sequences was built to identify identical VT accessions. Identical VTs were merged together resulting in 15 OTUs used in subsequent analysis.

<sup>1</sup>https://qiime2.org

TABLE 2 | Primer sets and PCR conditions used in nested PCRs for "total fungi" and "AMF specific" amplicon sequencing.

#### Statistical Analysis

fpls-10-00895 July 11, 2019 Time: 17:37 # 6

Alpha diversity metrics [OTU/VT richness, Shannon's Index (e), Peilou's Evenness and Simpson's Diversity] were calculated using Usearch10's (Edgar, 2010) alpha diversity command and all subsequent analyses were carried out in R version 3.3.2 (October 31, 2016) (R Core Team, 2016). To test the effect of the aphid treatment on AMF and total fungi alpha diversity metrics, AMF family relative abundance, and plant and AMF structures, aphid presence was used as an explanatory variable in a linear model, using the R packages "lme4," "lmerTest," "lsmeans," and "car." Plot location, measured as the NE and NW distance of each plot from the NW corner of the field site, was used as a covariate (see **Table 3**). To identify fungal OTUs/VTs whose presence predicts aphid treatment we employed indicator species analysis using the "indicspecies" package. The identification and estimation of abundances of AMF is more accurate using the AMF specific amplicon (Berruti et al., 2017), thus total fungi OTUs corresponding to AMF taxa were excluded from the indicator species analysis of the total fungi amplicon. The effect of aphid presence on community composition between samples as measured via Bray-Curtis dissimilarity (Beta diversity) was analyzed using a PERMANOVA via the "Adonis" function in the "Vegan" R package. This was visualized via nonmetric multidimensional scaling (NMDS). Relationships between community composition and plant biomass and nutrition, and AMF structures as well as plot location were tested by applying the "envfit" function to the NMDS.

#### RESULTS

#### Effect of Aphid Presence on Plant Nutrition and Biomass, and AMF Structures

The mean number of aphids per tiller ± SE in "+Aphid" plots was 24.9 ± 8.2, and ranged from between 1.8 and 75, whilst no aphids were present on any tillers investigated for "- Aphid" plots. Whilst the majority of plant biomass and nutrition, and AMF physiological traits were not influenced by aphid treatment (**Table 3**), aphid treatment significantly increased the stem N:P ratio by 5.4% (**Table 3**), and there was a near significant increase of 8.9% (P = 0.052) of the grain N:P ratio. Grain P concentration and N:P ratio was associated with the position of the plot in the NE direction of the site (**Table 3**). While aphid treatment did not significantly impact any AMF structures, there was a trend that the vesicle frequency almost doubled in the roots of plants hosting aphids (P = 0.097). The frequency of vesicles was associated with the position of the plot in the NW direction (**Table 3**).

#### Effect of Aphid Presence on Total Fungal and AMF Communities in the Root

Within the entire fungal community in plant roots across both aphid treatments, 153 "total fungi" OTUs were identified from nine fungal phyla and 2 OTUs that could not be assigned at the phylum level. The highest abundance of sequences were assigned to Ascomycota (90.2%), followed by Basidiomycota (5.3%), unclassified fungi (3.5%), Glomeromycota (0.83%), and Chytridiomycota (0.08%). In contrast, sequences from Rozellomycota, Mortierellomycota, Entomophthoromycota, Mucoromycota, and Zoopagomycota contributed less than a combined 0.1% of sequence abundance. Within the AMF specific amplicon, 12 OTUs were assigned to VTs whilst three OTUs could not be assigned to a singular VT (see **Supplementary Table S1**). These VTs belonged to seven AMF families: Glomeraceae (6), Paraglomeraceae (1), Diversisporaceae (2), Ambisporaceae (1), Gigasporaceae (1), Archaeosporaceae (3), and Acaulosporaceae (1).

Aphid presence did not affect the species richness within the entire fungal community, but did increase its evenness (**Table 4**). Aphid presence also did not affect any AMF specific alpha diversity metrics. The Simpson's diversity of the AMF specific community was linked to plot location in the NE direction of the field site (**Table 4**).

Aphid presence had no effect on the relative abundance of AMF reads within the entire fungal community. Within the AMF specific amplicon the relative abundance of the Gigasporaceae family tended to increase when aphids were present (**Table 5**). The relative abundance of AMF reads within the entire fungal community, and the relative abundance of the Gigasporaceae and Ambisporaceae families present in the AMF specific community were associated with the location of the plot in the NE direction of the site (**Table 6**).

All but two AMF specific amplicon VTs were found in both aphid and no aphid treatments, however, the exceptions were not strong indicators of aphid presence or absence (**Table 6**). Several Ascomycota taxa were strong predictors of the absence of aphids for the total fungi amplicon, whilst a member of the family Cystobasidiaceae indicated aphid presence.

The total fungi amplicon community composition between plots, measured as Bray-Curtis dissimilarity (beta diversity) was not significantly affected by aphid presence (PERMANOVA: F1,<sup>13</sup> = 1.75, P = 0.131). However, community composition correlated with plot location in the NE direction of the site (R <sup>2</sup> = 0.46, P = 0.031) and grain P concentration (R <sup>2</sup> = 0.46, P = 0.021; **Figure 3**).

The community composition of AMF specific VTs in the roots of barley plants between plots was not affected by the presence of aphids (F1,<sup>13</sup> = 0.46, P = 0.604). However, the environmental factors of stem Si concentration and location of the plots in the NE direction of the field site were significantly correlated with the community composition (R <sup>2</sup> = 0.39, P = 0.049 and R <sup>2</sup> = 0.40, P = 0.049, respectively; **Figure 4**).

#### DISCUSSION

This study aimed to investigate the impact of aphids on soil fungi in a conventionally managed agricultural system. It was hypothesized that by supressing plant nutrition and growth, aphid feeding would lead to negative impacts on AMF structures and species richness, as well as the evenness

TABLE 3 | Mean (±SE) above-ground plant biomass and nutrient concentrations, and AMF structures of experimental barley plots treated with or without aphids, using the location of the plot in the NE and NW direction of the site as a model covariate.

P-values highlighted in bold are significant at a value of <0.05, P-values highlighted with · are >0.05 and <0.1. <sup>1</sup>No., Number; <sup>2</sup>DW, dry weight (g); <sup>3</sup> [], concentration (mg g−<sup>1</sup> ); <sup>4</sup>RLC, AMF root length colonization where Min, most conservative estimate; Max, least conservative estimate; RLC, arbuscule and vesicle values are % of root length; <sup>5</sup>HLD, hyphal length density (m hyphae g−<sup>1</sup> soil).

TABLE 4 | Mean (±SE) Alpha diversity metrics for "total fungi" amplicon OTUs, and "AMF specific" amplicon VTs from barley plots without or without the presence of aphids, using the location of the plot in the NE and NW direction of the site as a model covariate.

Richness: number of OTUs/VTs. P-values highlighted in bold are significant at a value of <0.05, P-values highlighted with · are >0.05 and <0.1.

Plot location covariate Aphid presence NE NW Relative abundance −Aphid +Aphid F1,<sup>11</sup> P F1,<sup>11</sup> P F1,<sup>11</sup> P AMF 0.53 ± 0.29 1.18 ± 0.63 1.450 0.254 5.447 0.040 0.145 0.710 Acaulosporaceae 0.02 ± 0.02 0 – – – – – – Ambisporaceae 3.21 ± 0.79 2.92 ± 0.71 0.123 0.732 8.111 0.016 0.036 0.854 Archaeosporaceae 1.23 ± 0.41 0.63 ± 0.19 1.828 0.204 0.056 0.817 1.848 0.201 Diversisporaceae 3.72 ± 0.68 5.68 ± 2.52 0.767 0.400 2.241 0.163 0.198 0.665 Gigasporaceae 0.77 ± 0.20 1.48 ± 0.52 4.106 0.067· 13.560 0.004 0.553 0.473 Glomeraceae 85.44 ± 2.18 83.43 ± 2.30 0.420 0.530 0.778 0.397 <0.001 0.990 Paraglomeraceae 5.62 ± 1.32 5.86 ± 1.54 0.016 0.900 0.658 0.434 1.119 0.313

TABLE 5 | Mean (±SE) relative abundances of the AMF sequences within the "total fungi" amplicon, and of AMF family sequences from the "AMF specific amplicon," from experimental barley plots infested with or without aphids, using the location of the plot in the NE and NW direction of the site as a model covariate.

P-values highlighted in bold are significant at a value of <0.05, P-values highlighted with · are >0.05 and <0.1.

of the AMF and other fungi communities, measured here through analysis of an AMF specific amplicon and another less precise, but wider encompassing amplicon (total fungi). It was also proposed that the effects of aphids would impact the compositions of these communities (Gehring and Bennett, 2009), and the relative abundance of AMF taxa within them. In fact, aphid presence had less of an impact on AMF community structure than the location of the host plant in the field, although there was a trend of the abundance of AMF vesicles, and the abundance of the Gigasporaceae family to increase when aphids fed on the host plant. Within the total fungi community, the relative abundance of AMF was also affected by location rather than aphid presence. However, aphid presence increased evenness across the total fungi community.

Only total fungi amplicon OTUs which significantly indicate (bold) the aphid treatments are shown. X indicates the presence of the taxon, while × indicates its absence.

### Effects of Aphids on the Plant Biomass and Nutrition, and AMF Structures

Contrary to our hypothesis, S. avenae had little effect on the above-ground nutrition of barley in the field, although aphid presence tended to increase the above-ground plant N:P ratio (**Table 3**), possibly due to nutrient re-allocation caused by aphid feeding (Sandstrom et al., 2000; Thompson and Goggin, 2006; Nowak and Komor, 2010), or differences between the requirement for N and P by aphids (Tao and Hunter, 2012). Moreover, aphid presence did not reduce above-ground biomass in the present study. Since the aphids used in this experiment were cultured under controlled conditions it is likely that as the aphids did not vector any plant viruses that are a major contributor to aphid related yield loss in cereals (McKirdy et al., 2002). While it is possible that the agro-chemical inputs of fungicides, herbicides and plant growth regulators in this conventionally managed system may have influenced aphid development, the total number of aphids per tiller, and thus, estimated in each "+ Aphid" plot, remained high. However, the aphid population in the current experiment was launched at an earlier date than in non-controlled systems (Blackman and Eastop, 2000) which may have also influenced the results. AMF colonization of barley roots was high compared to that measured in field studies previously (Boyetchko and Tewari, 1995), and in certain glasshouse studies (Grace et al., 2009). Thus, while we originally hypothesized that aphids would have a negative impact on above-ground plant biomass and nutrition which in turn, would reduce both the internal and external phases of the AMF, no negative impact occurred. This may explain why both AMF RLC and hyphal length density were not affected by the presence of aphids in this study. However, there was a trend for aphids to increase the proportion of vesicles in plant roots. Mechanical defoliation has also been shown to influence the proportion of vesicles in plant roots when grown with a native AMF soil community (Garcia and Mendoza, 2012). As vesicles are lipid storage organs in AMF, and AMF derive lipids from the host plant (Keymer and Gutjahr, 2018), this might suggest that more fixed C is available to the AMF via the plant under aphid herbivory.

FIGURE 3 | Total fungi community composition: Non-Metric Multidimensional Scaling (NMDS) of total fungi community distribution based on Bray-Curtis dissimilarity (stress = 0.158) obtained from the abundance of "total fungi amplicon" OTUs under "+Aphid," or "–Aphid" treatments. Total fungi community composition was significantly associated with the grain P concentration (mg g-<sup>1</sup> ) in the plots (Grain [P]) and the location of the plot in the NE direction of the site (NE).

Alternatively, AMF vesicle numbers may have increased in a response to less C flow from the plant in order to aid survival under stressed conditions, similar to the AMF response observed in cool conditions (Hawkes et al., 2008), It is also possible that the frequency of AMF with a propensity to form more vesicles occurred within the Glomeracea community. However, the tendency for aphid presence to increase the relative abundance of the AMF family Gigasporaceae, which do not develop vesicles (Smith and Read, 2008), suggests that this explanation is less likely.

### Effects of Aphids on AMF and Total Fungal Communities

The second hypothesis proposed that the effects of aphid feeding on the host plant would reduce AMF and total fungal species richness and evenness. The low richness of AMF species identified in the current study is similar to that documented in other conventionally managed barley monoculture systems (Manoharan et al., 2017), and as the read depth achieved for the AMF specific amplicon is sufficient to capture AMF diversity (Vasar et al., 2017), it can be assumed that this is an accurate representation. The number of species of AMF were not impacted by aphid presence perhaps as there was no effect of aphid presence on above-ground plant biomass; however, this is similar to the lack of an effect of arthropod feeding on pinyon pine associated ectomycorrhizal communities (Gehring and Bennett, 2009). This was also reflected in the number of species in the entire fungal community in the current study. Moreover, and counter to expectations, aphid presence increased the evenness of the entire fungal community in the present study. Nutrient flows below-ground were not measured in the current study, but aphids can affect below-ground respiration depending on plant growth stage, potentially due to alterations in C availability to soil microbes (Vestergård et al., 2004). Aphids can also alter the profiles of compounds, released from plant roots (Hoysted et al., 2018) and can also change the profiles of sugars found in AMF hyphae sharing the same host plant (Cabral et al., 2018). Moreover, aphids excrete honeydew as a result of their C rich diet of phloem sap which can be utilized as a C source by soil microbiota, thus shaping community structure and biomass (Katayama et al., 2014; Milcu et al., 2015). As more C sources become available in the root, it is possible that niches may enlarge allowing less competitive fungi to compete, reducing the dominance of abundant taxa. However, it should be noted that aphid induced alterations to soil organisms can occur independently of honeydew C inputs (Sinka et al., 2009), and that soil microbes can be influenced by aphid induced changes to plant root exudates in systems where honeydew does not reach the microbe (Kim et al., 2016). Above-ground herbivory generally stimulates the cycling of nutrients by decomposers in the soil (A'Bear et al., 2014), which could partly explain increases in the N:P ratios of plant tissues infested by aphids in the current study.

It was hypothesized that the abundance of AMF taxa within the AMF community would be impacted by aphid presence, and there was a marginal increase in the abundance of Gigasporaceae under aphid infestation of the host plant. A recent metaanalysis revealed that members of this family are tolerant to fertilizer input disturbances, suggesting a role aside from nutrient acquisition, perhaps in plant defense (van der Heyde et al., 2017b). Species indicator analysis may also identify taxa affected by treatments, however, as low abundance taxa typically score as poor indicators, the results of this method may differ from those investigating relative abundances (Longa et al., 2017). None of the AMF VTs were indicators of either treatment, but several total fungi amplicon OTUs were significant indicators of aphid presence or absence. Currently, it is unclear whether these fungi are responding to changes in nutrient availability, whether the plant recruits them in response to aphid feeding to aid with defense (Kong et al., 2016; Pineda et al., 2017), or whether the recruitment of specific soil microbes ultimately benefits the aphid (Kim et al., 2016).

No clear effects of aphid presence were found on community composition, perhaps as a longer period of top down pressure is required to impact this metric. For example, the effects of grazing by large vertebrates on AMF community structure are

strongly linked to the length of the grazing (van der Heyde et al., 2017a). However, as plant communities are removed regularly in cereal systems, and aphid feeding is seasonal (Blackman and Eastop, 2000), only a relatively short window is available for these interactions to occur. Plot location, grain P concentration and AMF abundance as well as community composition were tightly linked in the current study. AMF communities may be associated with environmental and nutritional gradients in soil systems (Bainard et al., 2014; Horn et al., 2014; Guo et al., 2016; Yao et al., 2018) which also likely reflects the spatial differences measured in above-ground plant nutrition and AMF physiology here. Although care was taken to reduce the spread of plots across the slope of the site, these spatial differences highlight the importance of environmental heterogeneity within relatively small distances (less than 10 m) in field sites, and this could have contributed to masking top down effects of aphids on community composition in the current study. Whether the associations between AMF community composition and plant Si are due to AMF uptake of Si (Garg and Bhandari, 2016; Frew et al., 2017), or an artifact of AMF responses to soil pH gradients or water availability (Bainard et al., 2014) requires further study.

#### CONCLUSION

Aphids increased the evenness of the entire fungal community within plant roots, and also tended to increase the level of vesicles and abundance of the AMF family Gigasporaceae. Whether these increases are due to increased C allocation below-ground by plants attempting to increase nutrient uptake, or the active selection of fungal taxa in response to herbivory requires elucidation. Whether these changes in the below-ground soil community feed back into altered aphid performance is currently unclear, but the response of agriculturally relevant fungal communities to top-down effects of herbivory suggests that above-below-ground community feedback could occur in agricultural systems.

#### REFERENCES

### DATA AVAILABILITY

The raw sequence datasets for this study are available in the European Nucleotide Archive accession number PRJEB31469. The datasets for this study are available in the Dryad digital depository (doi: 10.5061/dryad.864pt1f/1).

#### AUTHOR CONTRIBUTIONS

TW, AH, SH, and JF contributed to the experimental design of the study. TW conducted the experiments and analyzed the data with consultation of the bioinformatics analysis from J-PM. All authors contributed to the writing of the manuscript.

### FUNDING

This work was funded by the Biotechnology and Biological Sciences Research Council (BBSRC) United Kingdom, Doctoral Training Partnership grant (BB/J014443/1).

#### ACKNOWLEDGMENTS

We thank the RSK A.D.A.S Ltd. (High Mowthorpe, North Yorkshire, United Kingdom), in particular Dr. Kate Storer for advice and facilitation of the field study reported. We also thank members of the Dr. Thorunn Helgason group, particularly Dr. Susannah Bird and Phillip Brailey for advice and facilitation regarding sequencing and bioinformatics analysis. This work was derived from the work presented in a Ph.D. thesis of Wilkinson (2018).

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.00895/ full#supplementary-material

by abiotic factors and host crop in a semi-arid prairie agroecosystem. FEMS Microbiol. Ecol. 88, 333–344. doi: 10.1111/1574-6941.12300

long-term grazing exclusion on semi-arid steppes in inner Mongolia. Geoderma 269, 79–90. doi: 10.1016/j.geoderma.2016.01.039

spectrometer. New Phytol. 195, 699–706. doi: 10.1111/j.1469-8137.2012.04179. x

and aphid performance. New Phytol. 163, 169–176. doi: 10.1111/j.1469-8137. 2004.01106.x

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Wilkinson, Miranda, Ferrari, Hartley and Hodge. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Soil Amendment With Different Maize Biochars Improves Chickpea Growth Under Different Moisture Levels by Improving Symbiotic Performance With Mesorhizobium ciceri and Soil Biochemical Properties to Varying Degrees

#### Edited by:

Dilfuza Egamberdieva1,2,3

\* † , Li Li<sup>3</sup>†

Sonoko Dorothea Bellingrath-Kimura<sup>1</sup>

Sabine Dagmar Zimmermann, Centre National de la Recherche Scientifique (CNRS), France

#### Reviewed by:

Monica Calvo Polanco, Montpellier SupAgro, France Marta Laranjo, University of Evora, Portugal

#### \*Correspondence:

Dilfuza Egamberdieva<EMAIL_ADDRESS><EMAIL_ADDRESS>
†These authors have contributed equally to this work

#### Specialty section:

This article was submitted to Plant Microbe Interactions, a section of the journal Frontiers in Microbiology

Received: 05 March 2019 Accepted: 07 October 2019 Published: 31 October 2019

#### Citation:

Egamberdieva D, Li L, Ma H, Wirth S and Bellingrath-Kimura SD (2019) Soil Amendment With Different Maize Biochars Improves Chickpea Growth Under Different Moisture Levels by Improving Symbiotic Performance With Mesorhizobium ciceri and Soil Biochemical Properties to Varying Degrees. Front. Microbiol. 10:2423. doi: 10.3389/fmicb.2019.02423 <sup>1</sup> Leibniz Centre for Agricultural Landscape Research, Müncheberg, Germany, <sup>2</sup> Department of Microbiology, Faculty of Biology, National University of Uzbekistan, Tashkent, Uzbekistan, <sup>3</sup> CAS Key Laboratory of Biogeography and Bioresource in Arid Land, Xinjiang Institute of Ecology and Geography, Ürümqi, China

, Stephan Wirth<sup>1</sup> and

, Hua Ma<sup>1</sup>

Chickpea (Cicer arietinum L.) is an important legume originating in the Mediterranean and the Middle East and is now cultivated in several varieties throughout the world due to its high protein and fiber content as well as its potential health benefits. However, production is drastically affected by prevalent water stress in most soybean-growing regions. This study investigates the potential of biochar to affect chickpea-Rhizobium symbiotic performance and soil biological activity in a pot experiment. Two different biochar types were produced from maize using different pyrolysis techniques, i.e., by heating at 600◦C (MBC) and by batch-wise hydrothermal carbonization at 210◦C (HTC), and used as soil amendments. The plant biomass, plant nutrient concentration, nodule numbers, leghemoglobin (Lb) content, soil enzyme activities, and nutrient contents of the grown chickpeas were examined. Our results indicated that plant root and shoot biomass, the acquisition of N, P, K, and Mg, soil nutrient contents, soil alkaline and acid phosphomonoesterases, and proteases were significantly increased by HTC char application in comparison to MBC char under both well-watered and drought conditions. Furthermore, the application of both biochar types caused an increase in nodule number by 52% in well-watered and drought conditions by improving the symbiotic performance of chickpea with Mesorhizobium ciceri. Rhizobial inoculation combined with HTC char showed a positive effect on soil FDA activity, proteases and alkaline phosphomonoesterases under well-watered and drought conditions compared to the control or MBC char-amended soils. This concept, whereby the type of producing biochar plays a central role in the effect of the biochar, conforms to the fact that there is a link between biochar chemical and physical properties and enhanced plant nutrient acquisition, symbiotic performance and stress tolerance.

Keywords: maize biochar, loamy sand, chickpea, rhizobia, nutrient acquisition, soil enzymes

## INTRODUCTION

fmicb-10-02423 October 29, 2019 Time: 15:58 # 2

The long-term use of inorganic nitrogen fertilizers, or the application of such fertilizers at rates higher than the optimum, increases residual inorganic N with adverse effects such as soil degradation or the decline of soil biological health (Gong et al., 2013; Reckling et al., 2016). An alternative N resource to mitigate such drawbacks is biological nitrogen fixation by legumes, which plays a major role in sustaining or improving soil productivity (Santi et al., 2013; Egamberdieva et al., 2018). Legumes are grown in many countries of the world and are considered an important food source for human and animal nutrition (Lüscher et al., 2014; Egamberdieva et al., 2015). However, abiotic stresses such as drought and salinity threaten the growth and yield of legumes and other crops (Bodner et al., 2015; Egamberdieva et al., 2016a, 2017b; Abd Allah et al., 2017). Grain legumes are sensitive to drought stress, and the symbiotic performance of plants with rhizobia is known to be especially hindered by drought, resulting in decreased nodulation and nitrogen fixation (Bouhmouch et al., 2005; Farooq et al., 2017). Inhibition of root nodule formation in legumes can be attributed to the failure of rhizobial colonization in the rhizosphere, which indicates the susceptibility of bacterial proliferation to stress factors (Rehman and Nautiyal, 2002).

Chickpea (Cicer arietinum L.) is an important legume crop in many countries and is considered a functional food source, mostly due to its high protein content (17–31% protein) (Merga and Haji, 2019). Water deficit reduces chickpea growth and yield and leads to low N fixation due to the limitation of rhizobium survival and growth in the soil (Rahbarian et al., 2011). Furthermore, a reduction in photosynthetic pigments, CO<sup>2</sup> assimilation rate, leaf water contents and disturbance in nutrient acquisition by plants were found under drought stress (Vadez et al., 2012; Siddiqui et al., 2015).

Different strategies have been developed to improve the tolerance of crops to drought, e.g., breeding for stress tolerant varieties (Witcombe et al., 2008), genetic engineering (Roy et al., 2011), and application of microbial inoculants (Hashem et al., 2016; Parray et al., 2016). Biochar amendment has also been repeatedly discussed as a technique to increase plant tolerance to various biotic and abiotic stresses (Akhtar et al., 2015). Biochar has been widely used as a soil amendment to improve soil fertility through improved water holding capacity (Yu et al., 2013; Lim et al., 2016), soil cation exchange capacity (Novak et al., 2009), nutrient retention, especially in soils deficient in organic matter content (Chan et al., 2007; Ibrahim et al., 2013), or soil microbial activities (Kolton et al., 2011; Egamberdieva et al., 2016b, 2017a; Soudek et al., 2016). Several reports have demonstrated improved plant growth of soybean (Glycine max L.) (Egamberdieva et al., 2016b), pepper (Capsicum annuum L.), tomato (Lycopersicum esculentum Mill.) (Graber et al., 2010), maize (Islami et al., 2011), and wheat (Alburquerque et al., 2014) after biochar amendment. Moreover, nodule number in the cases of soybean (Mete et al., 2015) and faba bean (Mohamed et al., 2017) significantly increased due to the addition of biochar to the soil, indicating improved symbiotic performance of the plant with rhizobia. Concerning soil biochemical properties, enzyme activities were studied to rate or monitor soil fertility and productivity (Dempster et al., 2012; Song et al., 2018). Several studies reported an increase in soil enzyme activity after biochar application; however, others reported a decrease in activity (Ameloot et al., 2013; Paz-Ferreiro et al., 2014). There is evidence accumulating that the response of plant growth, nutrient acquisition, soil biogeochemical processes and microbial communities to biochar depends on the type of feedstock and the thermal conditions during production (e.g., heating period and the final set point temperature) (Gul and Whalen, 2016).

According to Deenik et al. (2011), the chemical composition of biochar, especially pyrolysis conditions, plays a vital role in soil biological responses to biochar soil amendments. For example, Butnan et al. (2015) observed that plant growth was reduced in a sandy Ultisol amended with eucalyptus wood-derived biochar produced by high (800◦C) temperatures, whereas biochar produced at 350◦C enhanced plant growth. Moreover, the effect of biochar characteristics such as pyrolysis temperature and duration and the application rate of biochar on soil biological activities is highly variable (Paz-Ferreiro et al., 2014; Song et al., 2018).

In another study, straw gasification biochar increased the shoot and root growth of barley (Hordeum vulgare L.) in sandy soil compared to the addition of wood gasification biochar (Hansen et al., 2016). In addition, the interaction of biochar with environmental conditions is essential for determining any contrasting effects, which also depend on the physicochemical properties of biochar (Joseph et al., 2010). Therefore, elucidation of the effect of biochar type on plant growth, development and soil biochemical properties provides important guidance on the selection of feedstock type and production technology, which could be applied under specific environmental conditions. Only limited information is available in the literature regarding the effect of biochar types on the symbiotic performance of legume plants with rhizobia and plant growth under drought stress conditions.

We hypothesized that amending soil with biochar might alleviate drought stress in chickpea by improving its symbiotic interactions with rhizobia in a loamy, sandy soil. Moreover, considering that sandy loam has a pH value below 7, it was unknown whether using different biochar types would promote soil biological activity. We selected biochars produced from maize by two different pyrolysis techniques, i.e., heating at 600◦C (MBC) and batch-wise hydrothermal carbonization at 210◦C (HTC), and amended loamy, sandy soil under drought and well-watered conditions. The aim of the present study was (1) to evaluate the effect of two contrasting biochar types on chickpea growth, nutrient uptake (N, P, and K) and symbiotic performance with Mesorhizobium ciceri and (2) to determine the impact on soil nutrients and soil enzymes linked to carbon, nitrogen and phosphorus cycling.

### MATERIALS AND METHODS

fmicb-10-02423 October 29, 2019 Time: 15:58 # 3

#### Soil, Biochars and Plant Material

The soil used in the study was sandy loam collected from the upper horizon (0–15 cm depth) of an experimental arable field under irrigation operated by the Experimental Field Station of Leibniz Centre for Agricultural Landscape Research (ZALF), Müncheberg, Germany. The soil consists of clay and fine silt (7%), coarse and medium silt (19%), and sand (74%) and is characterized by the following properties: pH, 6.2; organic C content (0.55%), total N content (0.07%), P content (0.03%), K content (1.25%), and Mg content (0.18%).

Two biochar materials were used in this study, supplied from the Leibniz-Institut for Agrartechnik Potsdam-Bornim e.V. (ATB), Germany (Reibe et al., 2015). Both biochar materials were produced from maize using two different pyrolysis techniques: (i) heating at 600◦C for 30 min (MBC) and (ii) batch-wise hydrothermal carbonization at 210◦C and 23 bar for 8 h (HTC). The biochar properties are presented in **Table 1**. Chickpea (Cicer arietinum L.) seeds (var. Xalima) were obtained from the International Centre for Agricultural Research in the Dry Areas (ICARDA), Tashkent, Uzbekistan.

#### Plant Growth Experiment

The experiment was conducted in a plant growth chamber at the Leibniz Centre for Agricultural Landscape Research (ZALF), Germany. The following six treatments were set up under wellwatered and drought conditions, for a total of 12 treatments:

Two different biochar types derived from maize, namely, MBC and HTC, were used as soil amendments. Pots (d = 0.16 m, v = 3016 cm<sup>3</sup> ) were filled with 1000 g of air-dried soil mixed with crushed char (particle size < 3 mm) at a 2% (w/v) concentration. The chickpea seeds were surface-sterilized using 10% v/v NaOCl for 5 min and 70% ethanol for 5 min. Afterward, the seeds were rinsed five times with sterile distilled water and transferred to paper tissue for germination in a dark room at 25◦C for 3– 4 days. The bacterial strain Mesorhizobium ciceri HAMBI1750 was obtained from the Culture Collection of Helsinki University. The strain was grown in yeast mannitol broth (Sigma Aldrich, United States) for 4 days, 2 ml of the bacterial culture was pelleted by centrifugation (10,000 × g for 10 min), and the supernatant was discarded. The cell pellets were washed with 2 ml of PBS (20 mM sodium phosphate and 150 mM NaCl, pH of 7.0) and resuspended in PBS to obtain a final bacterial density of 10<sup>8</sup> CFU ml−<sup>1</sup> . Germinated seeds were dipped in bacterial suspension for 5 min and transferred to pots. Three seeds were sown in each pot, and after 1 week, the seedlings were thinned to one plant per pot.

Each treatment was replicated for three times and was arranged in a randomized complete block design. Only one plant was cultured in each pot for nodule and plant analysis. The soil moisture levels, i.e., drought (40% soil moisture) and well-watered (80% soil moisture) conditions, were monitored using the commercially available UMP-1 BT soil moisture sensor (Umwelt-Geräte-Technik GmbH, Müncheberg, Germany). Plants were grown for 40 days at a temperature of 24◦C/16◦C (day/night) and humidity of 50–60%. Well-watered and drought conditions were controlled at 80 and 40% soil moisture, respectively. At harvest, the chickpea root system was carefully removed, and soil adhering to the roots was collected and considered root-associated soil. The number of nodules (nodule size > 1 mm) was counted for each plant visually. The roots and shoots were oven-dried at 70◦C for 48 h. After determination of the dry weight, the materials were ground using a mill fitted with a 1 mm screen and then subsampled for analysis.

### Plant and Leghemoglobin Content Analyses

After 40 days, the root and shoot dry biomass and nodule numbers were determined. The LB content of nodules was determined by the following method of Wilson and Reisenauer (1963). First, 0.5 mg of crushed and ground nodule tissue was mixed into 3 ml of Drabkin's solution. The supernatant was transferred to a 10 ml tube after centrifugation at 500 g for 15 min. The nodule tissue was extracted twice more with 3 ml of Drabkin's solution and combined with the supernatants. The volume was adjusted to 10 ml with Drabkin's solution, mixed and centrifuged at 20,000<sup>∗</sup> g for 30 min. The assay was standardized with a freshly prepared solution of bovine hemoglobin.

#### Soil Nutrients and Soil Enzyme Measurements

The carbon (C), nitrogen (N), phosphorus (P), and potassium (K) concentrations in plant tissues were analyzed with an inductively coupled plasma optical emission spectrometer (ICP-OES; iCAP

TABLE 1 | Chemical characterization of chars (Reibe et al., 2015).

FM, fresh matter; DM, dry matter; MBC, maize biochar; HTC, hydrochar; EC, electrical conductivity.

6300 Duo). The total carbon (Ct) and nitrogen (Nt) of the soil samples were determined by the dry combustion method (Nelson and Sommers, 1996) using an elemental determinator (TruSpec CNS). P, K, and magnesium (Mg) were analyzed with an ICP-OES (iCAP 6300 Duo). Acid and alkaline phosphomonoesterase activities were assayed according to Tabatabai and Bremner (1969). Briefly, 0.5 g of moist soil was placed in a 15 ml vial, and 2 ml of MUB buffer (pH of 6.5 for the assay of acid phosphatase or pH of 11 for the assay of alkaline phosphatase) and 0.5 ml of p-nitrophenyl phosphate substrate solution (0.05 M) were added. The soil suspension was incubated in a water bath at 37◦C with 300 rpm shaking after the vial was capped. After 1 h of incubation, the vial was removed from the water bath, and 2 ml of NaOH (0.5 M), 0.5 ml of CaCl<sup>2</sup> (0.5 M) and 5 ml of distilled water were added to stop the reaction. One milliliter of soil suspension was centrifuged at 6500 rpm for 5 min. The concentration of p-nitrophenol (p-NP) produced in the assays of acid and alkaline phosphomonoesterase and phosphodiesterase activities was calculated from a p-NP calibration curve after subtracting the absorbance of the control at 400 nm wavelength using a Lambda 2 UV-VIS spectrophotometer (Perkin Elmer) (Acosta-Martínez and Tabatabai, 2011).

Protease activity was assayed using the method described by Ladd and Butler (1972). Briefly, 0.5 g of soil was weighed into a glass vial, and 2.5 ml of phosphate buffer (0.2 M, pH of 7.0) and 0.5 ml of N-benzoyl-L-arginine amide (BAA) substrate solution (0.03 M) were added. The vials were capped and shaken in a water bath at 37◦C with a rotary speed of 300 rpm for 1 h. After the incubation, 2 ml of KCl (2 M) was added to the vials to stop the reaction. One milliliter of the solution was pipetted into a microcentrifuge tube and centrifuged at 6000<sup>∗</sup> g for 10 min. A total of 0.5 ml of the clear supernatant was mixed with 4.5 ml of distilled water, 2.5 ml of sodium salicylate-NaOH and 1 ml of sodium dichloroisocyanide solution. The reaction mixture was incubated at room temperature for 30 min. The absorbance of the colored mixture was measured at 690 nm against a blank reagent solution in a spectrophotometer. Controls were used to follow the procedure described for the assay, with the exception that the BAA substrate was added after the incubation. The ammonium released was calculated by relating the measured absorbance at 690 nm to that of a calibration graph containing 0, 1.0, 1.5, 2.0, and 2.5 µg of NH<sup>4</sup> <sup>+</sup>–N mL−<sup>1</sup> .

The FDA hydrolytic activity assay was performed according to Green et al. (2006). A total of 0.5 mg of soil was added to a 50 ml vial, with the subsequent addition of 25 ml of sodium phosphate (0.06 M; pH of 7.6). Then, 0.25 ml of 4.9 mM FDA substrate solution was added to all vials. The tightly capped vials were mixed and incubated in a water bath at 37◦C for 1 h. A 1 ml soil suspension was transferred to a 1.5 ml tube and centrifuged at 8000 rpm for 5 min. The clear supernatant was measured at 490 nm against a blank reagent solution in a spectrophotometer. Controls were used according to the procedure described for the assay, but 0.25 ml of acetone was added instead of the FDA substrate solution. The concentration of fluorescein released was calculated by reference to a standard curve with 0, 0.001, 0.005, 0.05 and 0.15 mg of fluorescein.

## Statistical Analysis

The data were subjected to univariate analysis using a general linear model in the statistical software package SPSS 22 (SPSS Inc., Chicago, IL, United States). Multiple comparisons between treatments were tested at the p < 0.05 level using Tukey's honestly significant difference (HSD) test.

## RESULTS

#### Plant Shoot and Root Growth Well-Watered Conditions

The responses of chickpea to the two biochar types were different under well-watered conditions. The growth of chickpea was significantly (p < 0.05) higher in soil amended with HTC than in soil without HTC addition (**Figure 1**). The shoot and root biomass of un-inoculated chickpea increased by 36%, whereas the shoot and root growth of inoculated plants increased by 16 and 22%, respectively. However, there were no significant effects of MBC on the growth of un-inoculated or inoculated chickpea. MBC char improved the shoot biomass of un-inoculated and inoculated chickpea by 18 and 13%, respectively, but the effect was not significant. The root biomass of neither un-inoculated nor inoculated chickpea was changed in MBC-amended soil (**Figures 1**, **2**).

#### Drought Conditions

Overall, there were positive and significant effects of HTC on the root and shoot growth of both un-inoculated and inoculated chickpea grown under drought conditions. The soil amendment with MBC char produced more benefits to chickpea growth under drought stress than under the well-watered condition. The shoot growth was increased by 9% and root growth by 24% in both un-inoculated and inoculated plants (**Figure 3**). However, the stimulation was not significant when compared to control plants grown without biochar addition.

The interrelation of biochar × moisture level and biochar × inoculation showed no effects but moisture level × inoculation revealed significant interaction effects on shoot biomass (**Table 2**). Biochar × moisture level showed a significant interaction effect but biochar × inoculation and moisture level × inoculation showed no interaction effect on root biomass. Biochar × moisture level × inoculation showed a significant interaction effect on root biomass but not on shoot biomass.

#### Nodulation and Leghemoglobin Content Well-Watered Conditions

It is notable that both biochars improved the symbiotic performance of chickpea with M. ciceri under well-watered conditions. The nodule numbers of inoculated chickpea were 17 ± 4.01, while HTC addition to soil increased nodule numbers to 52 ± 8.1 per plant, and MBC, by 23 ± 7.2 per plant (**Figure 4A**). The LB content was significantly higher (9.0 mg gFW−<sup>1</sup> ) in the nodules of plants grown in HTCamended soil than in plants grown in control and MBC-amended soil (**Figure 4B**). There was no significant difference in the

FIGURE 1 | The shoot and root growth of chickpea grown in soil amended with HTC (produced from maize by batch-wise hydrothermal carbonization at 210◦C) and MBC (produced from maize through pyrolysis by heating at 600◦C) under well-watered conditions. The plants were un-inoculated (control) or inoculated with M. ciceri. Column means with different letters are significantly different based on Tukey's HSD test at P < 0.05.

nodule LB content between plants grown in MBC-amended soil and control plants.

#### Drought Conditions

The nodule numbers decreased in chickpea grown under drought conditions for all variants (**Figures 4A,B**). The symbiotic performance of chickpea with M. ciceri was inhibited by drought stress, reducing nodule numbers from 17.6 to 3.4 per plant. The biochar addition increased nodule numbers to 25 ± 4.2 in HTC char-amended soil and 14.6 ± 3.1 in MBC char-amended soil, which were higher than the nodule numbers in control plants (3.4 ± 2.51) (**Figure 4A**).

The LB content per g of fresh weight of nodule tissue was considerably lower in the nodules of plants grown

FIGURE 3 | The shoot and root growth of chickpea grown in soil amended with HTC and MBC under drought conditions. The plants were un-inoculated (control) or inoculated with M. ciceri. Column means with different letters are significantly different based on Tukey's HSD test at P < 0.05.

Significance denoted by <sup>∗</sup>p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001.

under drought conditions. The LB content in chickpea nodules was not changed in MBC-amended soil compared to control soil. However, there was a significant positive effect on the LB content in nodules of plants grown in soil with HTC addition.

The Interaction of biochar × moisture level showed no significant effect on LB content (**Table 2**).

#### Plant Nutrient Content Well-Watered Conditions

In general, the concentrations of N, P, K, and Mg in plant tissues were affected by the type of biochar and by plant growth conditions. In HTC char-amended soil, the nutrient concentration of un-inoculated and inoculated plants with M. ciceri was higher than that of plants grown in soil without

biochar addition. Significant (p < 0.05) increases in N, P, K, and Mg concentrations over those in the controls were observed after HTC amendment under the well-watered conditions, with increases of 15, 35, 17, and 51%, respectively. Only K and Mg concentrations increased significantly (20–26%) in shoot tissue. The soil amendment with HTC char positively affected inoculated chickpea as well, increasing mainly N, P, and K uptake (**Table 3**). In MBC-amended soil, only K concentrations in the roots and shoots of un-inoculated and inoculated plants were increased compared to those in control plants grown in soil without biochar addition. The K concentration in the root tissue of un-inoculated and inoculated plants was increased by 123 and 119%, respectively, whereas that in the shoot tissue was increased by 84 and 58%. Among the nutrients, only the N and Mg concentrations in the shoots of un-inoculated plants were significantly (p < 0.05) increased by MBC char (**Table 3**).

#### Drought Conditions

Furthermore, a positive effect of MBC char on the nutrient uptake of plants was observed under drought conditions. The soil amended with MBC showed an increased concentration of N, P, and K contents in plant roots inoculated with M. ciceri by 20, 16, and 79%, compared to control plants, respectively (**Table 3**). The highest concentration of K in roots and shoots was also observed under drought conditions, increasing by 86% for un-inoculated plants and 78% for plants inoculated with M. ciceri. A notable result was observed where the MBC char amendment reduced the Mg concentration in plant tissue under drought conditions. The acquisition of N, P, K, and Mg by roots and shoots was improved under drought conditions by HTC char amendment, and significant stimulation was observed in inoculated plants (**Table 3**). For example, the concentrations of C, N, P, K, and Mg in root tissue increased by up to 27, 33, 28, 23, and 26%, respectively.

#### Soil Carbon and Nutrient Contents Well-Watered Conditions

**Table 4** shows the soil nutrient concentrations under wellwatered and drought conditions with biochar treatment. The highest values of soil organic C, N, and Mg were observed in soil amended with HTC and MBC under both inoculation with M. ciceri and no inoculation. The lowest value was found in soil without biochar treatments. The soil organic C, N, P, K, and Mg concentrations in soil amended with HTC (inoculated with M. ciceri or without inoculation) were increased by up to 66, 29, 14, 44, and 44%, respectively, under well-watered conditions compared to uninoculated soil without biochar. The C, P, and Mg concentrations in the MBC-amended soil were twofold higher and the K concentration was threefold higher than those in the soil without MBC addition and inoculation with M. ciceri (**Table 4**). Our results reveal that the uptake of available P and K was significantly affected by MBC biochar amendment compared to soil without biochar and regardless of inoculation with M. ciceri. There were no significant differences in soil nutrients between treatments with inoculated and noninoculated plants.

#### Drought Conditions

The nutrient concentrations in soil without biochar amendment were similar between drought and well-watered soils, except for the K concentration, which was higher under drought conditions. Under drought conditions, both biochar amendments increased soil organic C as well as N content by up to 84 and 46% compared to soil without biochar, respectively. The P, K, and Mg concentrations in soil were not significantly affected by HTC amendment of soil; only the soil with MBC biochar showed higher contents, regardless of M. ciceri inoculation. Furthermore, the P concentration increased by 75%, that of K, by 186%, and that of Mg, by 22% in soil with inoculated plants and amended with MBC. The nutrient content did not differ significantly between the soil of plants inoculated with M. ciceri and that of uninoculated plants.

The interaction of biochar × moisture level showed significant effects on soil P, K, and Mg contents but no impacts on soil N content (**Table 5**). Biochar × inoculation showed significant interaction effects only on soil K content but no impacts on soil N, P, and Mg content. Moisture level × inoculation showed significant interaction effects on soil Mg content but no effects on either soil N, P, or K contents. There were no interaction effects of biochar × moisture level × inoculation on soil nutrients.

### Soil Enzymes

#### Well-Watered Conditions

Among the biochar treatments, adding HTC char to soil increased FDA activity in well-watered soil with inoculated and uninoculated plants by 37% (**Figure 5A**). The soil FDA activities did not differ between soil amended with MBC biochar and control soil without biochar (**Figure 5A**). The soil acid and alkaline phosphomonoesterase activities under chickpea were also affected by biochar amendment. A significant increase of 70 and 40% in acid phosphomonoesterase activity was observed under HTC and MBC amendments compared to soil without biochar, respectively. Alkaline phosphomonoesterase activity was increased by up to 30% in HTC-amended soil. Both enzyme activities were slightly but not significantly higher for chickpea inoculated with M. ciceri than for uninoculated controls.

Protease activity was slightly higher in HTC-amended soil than in the control and MBC-amended soil; however, the effect was not significant (**Figure 6**). There was no difference in protease activity between MBC-amended soil and soil without biochar.

TABLE 3 | Carbon and nutrient concentrations (%) in chickpea shoot and root tissue grown in soils after application of HTC and MBC chars under well-watered and drought conditions.

Plants were grown for 40 days after the start of the experiment. The different letters indicate significant differences based on Tukey's HSD test at P < 0.05.

TABLE 4 | Soil carbon and nutrient concentrations after application of HTC and MBC chars under well-watered and drought conditions.

Plants were grown for 40 days after the start of the experiment. The different letters indicate significant differences based on Tukey's HSD test at P < 0.05.

TABLE 5 | Interaction effects of biochar, moisture level, and inoculation on soil nutrients and soil enzyme activities.

Significance denoted by <sup>∗</sup>p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001.

#### Drought Conditions

Under drought conditions, soil FDA activity was higher under both biochar types, regardless of bacterial inoculation, than that of soil without biochar (**Figure 5B**). The HTC char amendment increased soil FDA activity by 44% and MBC by 24%, and soil FDA activity was higher with bacterial inoculation than without inoculation. Drought stress inhibited the acid and alkaline phosphomonoesterase activities of soil without biochar, regardless of inoculation with M ciceri (**Figures 7**, **8**). Only HTC showed a positive effect on acid and alkaline phosphomonoesterase activities, with an increase of 17 and 49% compared to the control, but the effects were not significant for uninoculated soil (**Figure 8**).

There were no changes in soil protease activity under drought stress after the biochar amendments. However, an increase of 30% was observed in soil with rhizobial inoculation after HTC char amendment compared to soil without biochar (**Figure 6**).

The interactions of biochar × moisture level, biochar × inoculation, and biochar × inoculation showed

FIGURE 5 | Effect of HTC and MBC on soil FDA activity under well-watered (A) and drought (B) conditions. Error bars (standard error) followed by a different letter within each column are significantly different at P < 0.05 based on Tukey's HSD test.

FIGURE 6 | Effect of HTC and MBC on soil protease activity under well-watered (A) and drought (B) conditions. Error bars (standard error) followed by a different letter within each column are significantly different at P < 0.05 based on Tukey's HSD test.

significant impacts on soil acid phosphomonoesterase activity while biochar × moisture level × inoculation showed no interaction effect on soil acid phosphomonoesterase activity (**Table 5**). No interactions were found for FDA, protease, and alkaline phosphomonoesterase activities.

#### DISCUSSION