Philippa Kaur

The GUTS modelFor completeness, we will prove the result for the most general GUTS model, which will also prove the result for all reduced forms.Theorem 1
For any model version within the GUTS framework, let (S(t; alpha )) denote the survival probability at time t for a given non-zero exposure profile (C_w(t)) scaled by some EMF value (alpha). For any chosen (x > 0) percentage effect (exposure-induced mortality), model end time (t_{E}) and background mortality (h_b) low enough such that (S(t_E; 0) > 0) there exists a unique EMF (alpha _*) such that$$begin{aligned} S(t_{E}; alpha _*) = left( 1 – frac{x}{100}right) S(t_E; 0), end{aligned}$$
(2)
(alpha _*) is the (hbox {LP}_{x}) for the exposure profile.
Baudrot and Charles10 calculated (LC_{50}) values for GUTS-RED-SD and GUTS-RED-IT. Their results implied the result of Theorem 1 for the main regulatory models. Our work makes the result explicit and generalises it to the whole GUTS framework. Another result of Theorem 1 is that the (hbox {LP}_{x}) is monotonically increasing with respect to x. For example, the (hbox {LP}_{50}) will always larger than the (hbox {LP}_{10}) for the same exposure profile. This result comes directly from (S4) in the SI.DEB modelsDue to the additional complexity of the DEB model we split the result into multiple theorems and proofs, starting by showing continuity and monotonicity of the damage ODE.
Theorem 2

Let
(C_w)
be some external concentration over time. Assume an effects model where the effects of higher exposure on growth and/or reproduction are always adverse (or zero) at all points in time. Then, defining the scaled damage ODE as
$$begin{aligned} begin{aligned} frac{D(t; alpha )}{dt} =&k_d(x_u alpha C_w – x_eD) – (x_G + x_R)D. end{aligned} end{aligned}$$
(3)

Then, for any combination of feedbacks (varvec{X} = [varvec{X}_u, 0, varvec{X}_G, varvec{X}_R]), damage is monotonically increasing with respect to (alpha), and is continuous with respect to (alpha) as long as changes to L and R are continuous. Moreover, damage is strictly monotonically increasing with respect to (alpha) whenever (D(t;1) > 0).
The limitation that (varvec{X}_e = 0) will be discussed in greater detail later. However, depending on the pMoA of the stressor we can extend the result of Theorem 2 slightly.
Corollary 3.1
If the pMoA of a substance directly affects reproduction and does not affect growth, i.e. (varvec{S} = [0, 0, 0, s_{R}, s_{H}]) then the results of Theorem 2holds for any combination of feedbacks.
Finally, we can step from the results of Theorem 2 and Corollary 3.1 to show the existence and uniqueness of a critical multiplier ((hbox {EP}_{x}) or (hbox {LP}_{x})) for growth, reproduction and survival.
Theorem 3
Consider the DEB-TKTD model of Jager11 and a substance such that at least one of Theorem 2or Corollary 3.1hold. Further, let (C_w(t)) be a non-zero exposure profile where the time of first exposure is before (t_1) as defined in Table 1. Then, for any chosen percentage effect level (x > 0) there exists a unique EMF (alpha _* >0) such that$$begin{aligned} min left( frac{L(t_{E}; alpha _*)}{L_(t_{E}; 0)}, frac{R_c(t_{E}; alpha _*)}{R_c(t_{E};0)}, frac{S(t_{E}; alpha _*)}{S(t_{E}; 0)}right) = 1 – frac{x}{100} end{aligned}$$
(4)
this (alpha _*) is the (hbox {EP}_{x}) (or (hbox {LP}_{x})) for the exposure profile (C_w(t)).

Table 1 Table of the state variables and pMoAs (including combinations of pMoAs) in the DEB-TKTD model11.Full size table
The monotonicity of effects on all state variables in the DEB model means that, for the conditions described in Theorem 3, the (hbox {EP}_{x}) (or (hbox {LP}_{x})) is also monotonically increasing with respect to x.We should note here that one can either setup an algorithm to find the critical multiplier value for growth, reproduction and survival individually and then select the minimum or setup the algorithm to directly find the minimum critical multiplier as in (4). Both will produce the same result, but the second approach is likely to be faster.One could argue that ERA should consider the combined effects of lethal and sublethal stress on the individual’s fitness. This is possible using the continuous form of the Euler–Lotka equation24$$begin{aligned} B(t) = int _0^t B(t-a) l(a)b(a) da, end{aligned}$$
(5)
where B(t) is the number of births at time t, l(a) is the fraction of females which survive to age a and b(a) is the birth rate for mothers of age a. For the offspring of a test population which all have the same age (as is the standard in long-term toxicity experiments) this integral collapses to a single point, (B(t-a) = 1) when (t=a) and zero elsewhere. The DEB model provides exactly the values which we need to calculate B(t). Namely$$begin{aligned} l(a) = S(a), quad b(a) = frac{d}{dt}R_c(a). end{aligned}$$One can now find the births per individual per time predicted by the DEB model as$$begin{aligned} B(t) = S(t)frac{d}{dt}R_c(t). end{aligned}$$
(6)
Integrating (6) over the duration of the experiment gives the expected number of offspring produced per female alive at the start of the test.There are two clear options for how to proceed. Firstly, one could calculate (int _0^{t_E} B(t) dt) for each EMF and compare it to the control, similar to finding (hbox {EP}_{x}) values for individual endpoints. Alternatively, one can use B(t) as the basis to estimate the intrinsic population growth rate25. This quantity provides an estimation of population growth based on the survival and fecundity over time of individuals. Indeed, it is listed as a potential output value in the experimental guidelines for standard Daphnia magna reproduction tests26. For the first of these options we offer an extension to Theorem 3.
Corollary 3.2
Consider a DEB-TKTD model and exposure profile such that Theorem 3holds. The number of expected offspring per female, given by$$begin{aligned} mathrm {B}(t_E; alpha ) = int _0^{t_E} S(t; alpha )frac{d}{dt}R_c(t; alpha ) dt end{aligned}$$has a unique (hbox {EP}_{x}) (alpha _*) such that$$begin{aligned}frac{mathrm {B}(t_E; alpha _*)}{mathrm {B}(t_E; 0)} = 1 – frac{x}{100}end{aligned}$$
Our results provide a rigid boundary to the applicability domain of the EMF approach both in terms of existence and uniqueness. Existence relies on the initial time in the profile when external concentration is non-zero, as described in Table 1. While it is important to know about these conditions, they will rarely inhibit an ERA, since long initial periods with zero exposure are uncommon.Cases where uniqueness cannot be guaranteed require more caution and it is unwise to use root-finding algorithms. In the next subsection we explore what can happen outside of this domain and provide suggestions for how to still produce a single reliable (hbox {EP}_{x}) value.Surface:volume scaling of eliminationThere is a reason that in Theorem 2, (varvec{X}_e = 0) was specified. In some cases when (varvec{X}_e = 1) a higher multiplier is no guarantee of higher damage for all time. Consider a substance which acts on assimilation and has surface area:volume scaled elimination (i.e. (varvec{X} = [0, 1, 0, 0])). The damage ODE under some EMF (alpha) is then$$begin{aligned} frac{dD}{dt} = k_d left( alpha C_w – frac{L_m}{L} D right) , end{aligned}$$where (L_m) is the maximum length the organism can reach. The EMF has a positive direct effect on damage, but also an opposing indirect effect. Increasing damage decreases the size of the organism which, due to the surface area:volume elimination of damage, enables faster elimination of damage. As a result, not only does Theorem 2 no longer hold but in fact a larger multiplier value can cause lower damage at some points in an exposure profile. In other words, we observe a paradoxical result whereby more exposure translates to less effect some time after exposure.Figure 3 illustrates what we will refer to as the “more is less” scenario. The exposure consists of a single pulse early in the animal’s life, modelled for two multiplier values, (alpha _2 > alpha _1). During the exposure phase the direct effect of the higher exposure causes higher damage and greater effects on size. After the pulse, external exposure is zero, and therefore the external concentration and uptake remain zero regardless of (alpha). Regardless of the EMF, scaled damage can only decrease during this phase. However, the effects of the higher multiplier are still relevant. As Fig. 2 shows, the feedback processes still influence damage dynamics. The model organism exposed to (alpha _2C_w) is smaller and therefore able to eliminate damage more rapidly because (varvec{X}_e = 1). This eventually leads to lower damage for the model organism exposed to (alpha _2C_w) (i.e. (D(t; alpha _1) > D(t; alpha _2))). The more is less phenomenon can also impact growth and cumulative reproduction, as seen in Fig. 3b,c. Sometime after exposure (L(t;alpha _2) > L(t; alpha _1)) and (R(t;alpha _2) > R(t; alpha _1)). For survival, and any additional endpoints without recovery, this “crossover” is unlikely, mortality during the exposure phase (where (D(t; alpha _2) > D(t; alpha _1))) will almost certainly dominate any mortality during the recovery phase. Figure 3d shows that for certain (x%) effect levels (vertical axis) multiple (hbox {EP}_{x}) values exist.Figure 3An illustration of the issues which can occur using the EMF approach for substances with surface area:volume scaled elimination (i.e. (varvec{X} = [0, 1, 0, 0])). The (non-multiplied) exposure is a constant (1 mu g/L) for the first 14 days and zero thereafter and effects assimilation only ((varvec{S} = [1, 0,0, 0, 0])). (a) Scaled damage, (b) length over time, (c) cumulative reproduction. (d) Endpoint value as a proportion of control after 40 days. The shape of these curves show that certain effect levels can be caused by two distinct multiplier values. Parameter values are (L_0 = 0.1), (f = 1), (r_B = 0.1), (L_p = 0.6), (L_m = 1), (R_m = 15), (kappa = 0.8), (y_P = 0.64) (z_b = 0.1), (b_b = 1), (k_d = 0.05), (varvec{X} = [0, 1, 0, 0]). See the SI for the definitions of these parameter values.Full size imageIn practice, instances of non-uniqueness such as Fig. 3 will be rare since they rely on a sudden and significant decrease in external exposure. Moreover, EMF methods for DEB-TKTD models will include a moving time window method18 consisting of many exposures constructed sequentially and assessed. Each window will produce an (hbox {EP}_{x}) value, but only the lowest will be relevant for the ERA. A time window which starts slightly earlier in the broader exposure profile would feature the same pulse later in the model organism’s lifespan and thus not allow organism recovery. Depending on the exact endpoint used, one would expect those windows to have a lower (and unique) (hbox {EP}_{x}). However, the potential for multiple (hbox {EP}_{x}) values raises concerns across all areas which impose a multiplicative margin of safety. We cannot guarantee that a multiplier resulting in (x%) effects exists nor that any value found by the algorithm is unique.Although not pictured here, maintenance and growth pMoAs and combinations of feedbacks which include (varvec{X}_e = 1) can also produce the “crossover” in the damage values and the “more is less” phenomenon seen in Fig. 3. It can also arise for scenarios which do not feature a deviation from the standard rules for growth (e.g. a starvation phase) and for other DEB based models. The SI features a similar plot to Fig. 3 showing damage crossover for a standard DEB model.Knowing this, the obvious question is how to proceed? Certainly with caution when (varvec{X}_e = 1) is necessary in model calibration and validation. Under such circumstances algorithms must ensure that the (hbox {EP}_{x}) value found is the lowest multiplier which gives (x%) effects when there is a risk of non-uniqueness. The brute force approach, incrementing from zero until the desired effect level is met or exceeded, is one example. Whether it is realistic for higher EMF values to cause reduced effects in vivo then does not alter the conservatism of the approach for ERA.Table 2 summarises the domain where the margin of safety approach can be used in conjunction with a root-finding algorithm without concern in the DEB-TKTD model of Jager11. For model configurations where non-uniqueness could emerge using another method to find the (hbox {EP}_{x}) is advisable. For example, a brute-force approach starting from an EMF of 0 in small increments (e.g. by 0.1). Without good reason, calibration should first be attempted with no feedbacks. Under this guiding philosophy of pursuing model simplicity we expect that the problem cases will be rare.Table 2 A table to mark under which scenarios the EMF approach is and is not guaranteed to produce a unique (hbox {EP}_{x}).Full size tableOther issuesThe damage crossover illustrated in the previous subsection occurs more commonly, and to a greater extent, when the pMoA is assimilation effects. This is because, at least in this standard implementation, stress can cause (100%) effect and completely cease assimilation when (s_A ge 1) (see SI for details). When this is the case, higher exposure (even from an increased multiplier) does not translate to higher stress. This differs from other pMoAs, whose stress values are unbounded. Indeed, replacing (1 – s_A) with (1/(1 + s_A)) in the model ((S5) in the SI) reduces the occurrence and scale of “crossovers” such as Fig. 3. However, the formulation of the pMoA should not be based on how it might affect the algorithm or the EMF.Certain species require further deviations from the standard model. For instance, different life-stages, growth and/or reproduction rules might be introduced to explain observed phenomena. Before models featuring these deviations are used in an EMF approach one should consider the potential issues as we have done in this section. While a proof of existence and uniqueness of the (hbox {EP}_{x}) for each model variant is ideal it is also infeasible. However, modellers should ensure that their approach is robust enough to deal with issues around existence and uniqueness. Checking that the model endpoint is reduced by (x%) when the (hbox {EP}_{x}) is applied to the exposure profile is an easy way to check accuracy and existence. An argument (if not a full, formal proof) for uniqueness should also be considered. In cases where that is not possible, the algorithm must be set up to identify the lowest (hbox {EP}_{x}), or check that no lower values exist.One common addition is to DEB-TKTD models which feature starvation is to assume that there is some maximum amount of starvation/shrinking which an animal can survive. Once that point is met or exceeded death is instantaneous27. Such death mechanisms cause problems. They can introduce a discontinuity in the response versus multiplier value for a given time window (i.e. a “jump” in plots such as Fig. 3). For instance, if in the example given in Fig. 3 the animal was not allowed to shrink, and instead died, then the multiplier of 6 would result in (100%) effects on survival (and significant growth effects). In contrast, the exposure when the multiplier is 2 is survivable and the animal can recover. Presumably, for some critical (alpha _c in (2, 6)) the exact threshold for death is reached. This (alpha _c) is a discontinuity between partial and (100%) effects relative to control. Under some circumstances this will prohibit finding a multiplier which results in exactly (x%) effects, regardless of the method used.There are two readily apparent solutions to this at the individual level. One is to set (alpha _c) as the multiplier for the window, the second is to replace such discrete responses with graded responses. In this example for instance, shrinking could add to the lethal hazard h. It is not possible to universally recommend one approach over the other as it will depend on the species’ behaviour. Once that decision has been made these issues must be recognised and reported by the modellers. More

The risks of open forest data exploitation are magnified by features of how forests are measured and who does the measuring. Generating long-term data on forest health and change involves physically measuring and identifying millions of trees. This means establishing, maintaining and revisiting plots, and curating records indefinitely. Trees are long-lived organisms so forests require decades of monitoring to properly infer change. Sustaining local observations for decades needs deep, long-term commitment to the unique but shifting combinations of people, institutions, regulations, interests and relationships that characterize each forest site. The challenge is enhanced by the great biodiversity of tropical forests. Measuring a single hectare of Amazon forest involves collecting and identifying up to ten times the number of tree species that are present in the UK’s entire 24 million hectares. There are very few people with the skills to do this.Long-term tropical-forest data measurements not only require effort and skill but also often carry risk and depend on some of the most disadvantaged actors in the global science community. Many forest workers (researchers, technicians, students, field assistants and local communities) lack basic job security, much less a career path, despite the long-term dedication that monitoring forests requires. In addition, many tropical forest workers may endure dangerous field conditions, with threats including kidnapping, armed insurgents, narcotraffickers, land-grabbers, infectious disease, snakebite, floods, fire, dangerous transport and gender-based violence. Besides these personal dangers, tropical scientists often lack the basic resources to measure and maintain their forest plots, let alone develop their research groups8.In contrast to the experiences of those monitoring forests on the ground, consider the context for satellite and aircraft-based measurements, which require ground-based data for validation. Space-based forest missions are expensive but are funded by public or private capital. Once in orbit, they stream data to analysts ‘for free’. This requires relatively few people to sustain, and although the analysts’ work is highly skilled, it carries little professional and physical risk and lacks commitment to place. Forest fieldwork is less capital-intensive, but needs sustained investment, is intensely human and carries substantial costs and risks. There are no automated collecting stations to help to identify and measure trees, so without the long-term dedication of many forest workers data collection simply stops.The risks and costs involved in acquiring and sustaining ground forest data are persistently overlooked, ignored or regarded as externalities to be picked up by the forest workers themselves. This is especially problematic because countries that hold the most tropical forests are among those least able to invest in science and development (Fig. 1, Supplementary Fig. 1). For example, monitoring the carbon balance of intact tropical moist forests has been estimated to cost US $7 million a year12, easily exceeding present support. By contrast, the USA alone spends over $90 million annually on its national forest inventory13. So, many tropical forest data are collected by skilled people working with minimal funding, in challenging conditions and facing other constraints, including complex layers of rules, agreements and research permits. Given such huge disparities, it is hardly reasonable to expect this output to be served on an open plate to the world.Fig. 1: Global distributions of per capita gross domestic product and tropical forest area.a,b, The 2008–2018 national average gross domestic product per capita (a) and tropical forest area per capita (b). Countries are coloured according to position from lowest (dark red) to highest (dark blue) within each global distribution.Full size imageIt is perhaps unsurprising that the most vocal proponents of making tropical and subtropical forest data open are often not those who actually measure and monitor them. Meanwhile, key beneficiaries include powerful publishers (usually with commercial interests), agencies and technology companies (often with commercial or political interests), and highly educated computer-savvy analysts wishing to integrate earth observation data with forest data (naturally with a career interest). Relatively few of these institutions and people are based in the tropics and subtropics. Fewer still are also data originators.And so, for many data originators the present meaning of making tropical forest data ‘open’ is to transfer the hard-won output of their labours to more privileged individuals and institutions, and lose more of the limited control they have over their professional lives. Power flows from the originators to public agencies, private companies and data scientists, mainly in the Global North. More

Lawrence, H. In City Trees: A Historical Geography from the Renaissance through to the Nineteenth Century (Charlottesville and London: University of Virginia Press, 2006, Lewis Mumford. The City in History: Its Origins, Its Transformations and Its Prospects (San Diego: Harvest Book Harcourt, 1961).Frawley, J. Campaigning for street trees, Sydney botanic gardens, 1890s–1920s. Environ. Hist. 15(3), 303–322. https://doi.org/10.3197/096734009X12474738199953 (2009).Article 

Google Scholar 
Seburanga, J. L., Kaplin, B. A., Zhang, Q.-X. & Gatesire, T. Amenity trees and green space structure in urban settlements of Kigali, Rwanda. Urban. For. Urban Green. 13(84–9313), 84–93. https://doi.org/10.1016/j.ufug.2013.08.001 (2014).Article 

Google Scholar 
Wilson, E. H. Northern trees in southern lands. J. Arnold Arbor. 4(2), 61–90 (1923).Article 

Google Scholar 
Gwedla, N. & Shackleton, C. M. Population size and development history determine street tree distribution and composition within and between Eastern Cape towns, South Africa. Urban. For. Urban. Gree. 25, 11–18. https://doi.org/10.1016/j.ufug.2017.04.014 (2017).Article 

Google Scholar 
Jacobs, A. B., Macdonald, E. & Rofé, Y. In The Boulevard Book: History, Evolution, Design of Multiway Boulevards (MIT Press, Cambridge, MA 2002), Robinson, W. The Parks and Gardens of Paris Considered in Relation to the Wants of Other Cities and of Private and Public Gardens (McMillan and Co., London , 1878).Akbari, A. H., Pomerantz, M. & Taha, H. Cool surfaces and shade trees to reduce energy use and improve air quality in urban. Sol. Energy. 70(3), 295–310 (2001).ADS 
Article 

Google Scholar 
Roy, S., Byrne, J. & Pickering, C. A systematic quantitative review of urban tree benefits, costs, and assessment methods across cities in different climatic zones. Urban For. Urban Green. 11, 351–363. https://doi.org/10.1016/j.ufug.2012.06.006 (2012).Article 

Google Scholar 
Schäffler, A. & Swilling, M. Valuing green infrastructure in an urban environment under pressure—The Johannesburg case. Ecol. Econ. 86, 246–257. https://doi.org/10.1016/j.ecolecon.2012.05.008 (2013).Article 

Google Scholar 
Santamour, F. S. Trees for urban planting: Diversity, uniformity and common sense. In Proceedings of the 7th Conference of the Metropolitan Tree Improvement Alliance (METRIA), vol. 7, 57–65 (1990).Shams, Z. I. Changes in diversity and composition of flora along a corridor of different land uses in Karachi over 20 years: caUses and implications. Urban. For. Urban Green. 17, 71–79. https://doi.org/10.1016/j.ufug.2016.03.002 (2016).Article 

Google Scholar 
Kambites, C. & Owen, S. Renewed prospects for green infrastructure planning in the UK. Plan. Prac. Res. 21(94), 483–496. https://doi.org/10.1080/02697450601173413 (2006).Article 

Google Scholar 
Cho, M. A., Malahlelac, O. & Ramoeloa, A. Assessing the utility WorldView-2 imagery for tree species mapping in South African subtropical humid forest and the conservation implications: Dukuduku forest patch as case study. Int. J. Appl. Earth. Obs. 38, 349–357. https://doi.org/10.1016/j.jag.2015.01.015 (2015).Article 

Google Scholar 
Niculescu, S., Lardeux, C., Grigoras, I., Hanganu, J. & David, L. Synergy between LiDAR, RADARSAT-2, and spot-5 images for the detection and mapping of wetland vegetation in the Danube Delta. IEEE J Sel. Top. Appl. Earth. Obs. Remote Sens. 9, 3651–3666 (2016).ADS 
Article 

Google Scholar 
Lefebvre, A., Picand, P.-A. & Sannier, C. Mapping tree cover in European cities: Comparison of classification algorithms for an operational production framework. In 2015 Joint Urban Remote Sensing Event (JURSE), IEEE, 1–4 (2015) https://doi.org/10.1109/JURSE.2015.7120511.Wyndham, C. H., Strydom, N. B., Van Rensburg, A. J. & Rogers, G. G. Effects on maximal oxygen intake of acute changes in altitude in a deep mine. J. Appl. Physiol. 29(5), 552–555 (1970).CAS 
Article 

Google Scholar 
Hegnauer, R. Chemotaxonomie der Pflanzen, vol. 3, 268–281 (Birkhäuser Verlag, Basel, 1964).Mabberley, D. J. The Plant-Book, 2nd edn. 87, 368–369 (Cambridge University Press, Cambridge, 1997).Gachet, M. S. & Schühly, W. Jacaranda—An ethnopharmacological and phytochemical review. J. Ethnopharmacol. 121, 14–27. https://doi.org/10.1016/j.jep.2008.10.015 (2009).CAS 
Article 
PubMed 

Google Scholar 
Gilman, E. F. & Watson, D. G. Jacaranda mimosifolia. Fact Sheet ST-317, Environmental Horticulture Department, Florida Cooperative Extension Service, University of Florida, Gainesville, http://www.ci.milpitas.ca.gov/_pdfs/council/2016/021616/item_04.pdf Accessed 6 June 2020 (1993).Dineva, S. B. Comparative studies of the leaf morphology and structure of white ash Fraxinus americana L. and London plane tree Platanus acerifolia Willd growing in polluted area. Dendrobiology 52, 3–8 (2004).
Google Scholar 
Liu, G., Li, Z. & Bao, M. Colchicine-induced chromosome doubling in Platanus acerifolia and its effect on plant morphology. Euphytica 157, 145–154. https://doi.org/10.1007/s10681-007-9406-6 (2007).Article 

Google Scholar 
Henry, A. & Flood, M. G. The history of the London plane, Platanus acerifolia, with notes on the Genus Platanus. Proc. R. Irish Acad Sect. B Biol. Geol. Chem. Sci. 35, 9–28 (1919).
Google Scholar 
Chavez, P. S. Image-based atmospheric corrections revisited and improved. Photogram. Eng. Rem. S. 62, 1025–1036 (1996).
Google Scholar 
Riano, D., Chuvieco, E., Salas, J. & Aguado, I. Assessment of different topographic corrections in Landsat-T. M. data for mapping vegetation types. IEEE Trans. Geosci. Remote Sens. 41, 1056–1061. https://doi.org/10.1109/TGRS.2003.811693 (2003).ADS 
Article 

Google Scholar 
Rouse J. W., Haas, R. H., Schell, J. A. & Deering, D. W. Proceedings of the Third Earth Resources Technology Satellite-1 Symposium, Greenbelt, USA: NASASP-351; 1974. Monitoring vegetation system in the great plains with ERTS, 3010–3017 (1974).R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/ (2021).Du, Y. et al. New hyperspectral discrimination measure for spectral characterization. Opt. Eng. 43(8), 1777–1786 (2004).ADS 
Article 

Google Scholar 
Bhattacharyya, A. On a measure of divergence between two statistical populations defined by their probability distributions’. Bull. Calcutta Math. Soc. 35, 99–109 (1943).MathSciNet 
MATH 

Google Scholar 
Bruzzone, L., Roli, F. & Serpico, S. B. An extension to multiclass cases of the Jefferys-Matusita distance. IEEE Trans. Pattern. Anal. Mach. Intell. 33, 1318–1321 (1995).
Google Scholar 
Kaufman, Y. & Remer, L. Detection of forests using mid-IR reflectance: An application for aerosol studies. IEEE Trans. Geosci. Remote Sens. 32(3), 672–683 (1994).ADS 
Article 

Google Scholar 
Padma, S. & Sanjeevi, S. Jeffries Matusita based mixed-measure for improved spectral matching in hyperspectral image analysis. Int. J. Appl. Earth. Obs. 32, 138–151. https://doi.org/10.1016/j.jag.2014.04.001 (2014).Article 

Google Scholar 
Kavzoglu, T. & Mather, P. M.. The use of feature selection techniques in the context of artificial neural networks. In Proceedings of the 26th Annual Conference of the Remote Sensing Society (CD-ROM), 12–14 September (Leicester, UK, 2000).Gunal, S. & Edizkan, R. Subspace based feature selection for pattern recognition. Info. Sci. 178, 3716–3726. https://doi.org/10.1016/j.ins.2008.06.001 (2008).Article 

Google Scholar 
Tolpekin, V. A. & Stein, A. Quantification of the effects of land-cover-class spectral separability on the accuracy of markov-random-field-based superresolution mapping. IEEE Trans. Geosci. Remote Sens. 47(9), 3283–3297. https://doi.org/10.1109/TGRS.2009.2019126 (2009).ADS 
Article 

Google Scholar 
Paterson, M., Lucas, R. M. & Chisholm, L. Differentiation of selected Australian woodland species using CASI data. In Proceedings IEEE International Geoscience and Remote Sensing Symposium, 643–645 (University of New South Wales, Australia, 2001).Richards, J. A. & Jai, X. Remote Sensing Digital Analysis: An Introduction, 4th edition (Springer, Berlin, 1999).Veraverbeke, S., Harris, S. & Hook, S. Evaluating spectral indices for burned area discrimination using MODIS/ASTER (MASTER) airborne simulator data. Remote Sens. Environ. 115, 2702–2709. https://doi.org/10.1016/j.rse.2011.06.010 (2011).ADS 
Article 

Google Scholar 
Breiman, L. Random forests. Mach. Learn. 45, 5–32 (2001).Article 

Google Scholar 
Georganos, S. et al. Geographical random forests: a spatial extension of the random forest algorithm to address spatial heterogeneity in remote sensing and population modelling. Geocarto Int. https://doi.org/10.1080/10106049.2019.1595177 (2019).Article 

Google Scholar 
Mellor, A., Haywood, A., Stone, C. & Jones, S. The performance of random forests in an operational setting for large area sclerophyll forest classification. Remote Sens. 5, 2838–2856. https://doi.org/10.3390/rs5062838 (2013).ADS 
Article 

Google Scholar 
Congalton, R. G. Accuracy assessment and validation of remotely sensed and other spatial information. Int. J. Wildland. Fire. 10, 321–328 (2001).Article 

Google Scholar 
Thomas, I. L., Ching, N. P., Benning, V. M. & D’aguanno, J. A. Review Article A review of multi-channel indices of class separability. Int. J. Remote Sens. 8(3), 331–350. https://doi.org/10.1080/01431168708948645 (1987).Article 

Google Scholar 
Mausel, P. W., Kramber, W. J. & Lee, J. K. Optimum band selection for supervised classification of multispectral data. Photogramm. Eng. Remote. Sens. 56(1), 55–60 (1990).
Google Scholar 
Singh, A. Some clarifications about the pairwise divergence measure in remote sensing. Int. J. Remote Sens. 5(3), 623–627. https://doi.org/10.1080/01431168408948845 (1984).Article 

Google Scholar 
Kumar, P. et al. A statistical significance of differences in classification accuracy of crop types using different classification algorithms. Geocarto Int. 32(2), 206–224. https://doi.org/10.1080/10106049.2015.1132483 (2017).Article 

Google Scholar 
McPherson, E. G., Simpson, J. R., Peper, P. J., Xiao, Q. & Wu, C. Los Angeles 1-Million Tree Canopy Cover Assessment. General Technical Report PSW-GTR-207. U.S. Department of Agriculture Forest Service Pacific Southwest Research Station. Albany, CA, 1–64 (2008).Rahimizadeh, N., Kafaky, S. B., Sahebi, M. R. & Mataji, A. Forest structure parameter extraction using SPOT-7 satellite data by object- and pixel-based classification methods. Environ. Monit. Assess. 192, 43. https://doi.org/10.1007/s10661-019-8015-x (2020).Article 

Google Scholar 
McRoberts, R. E. Satellite image-based maps: Scientific inference or pretty pictures?. Remote. Sens. Environ. 115, 715–724. https://doi.org/10.1016/j.rse.2010.10.013 (2011).ADS 
Article 

Google Scholar 
McRoberts, R. E. Probability- and model-based approaches to inference for proportion forest using satellite imagery as ancillary data. Remote. Sens. Environ. 114, 1017–1025. https://doi.org/10.1016/j.rse.2009.12.013 (2010).ADS 
Article 

Google Scholar 
Kokubu, Y., Hara, S. & Tani, A. Mapping seasonal tree canopy cover and leaf area using worldview-2/3 satellite imagery: A megacity-scale case study in Tokyo urban area. Remote. Sens. 12(9), 1505. https://doi.org/10.3390/rs12091505 (2020).Article 

Google Scholar 
Johannesburg City Parks and Zoo. 2018. The city that’s a rain forest. http://www.jhbcityparks.com/index.php/street-trees-contents-29. Accessed 14 June 2020.Tesfamichael, S. G., Newete, S. W., Adam, E. & Dubula, B. Field spectroradiometer and simulated multispectral bands for discriminating invasive species from morphologically similar cohabitant plants. GIsci. Remote Sens. 55(3), 417–436. https://doi.org/10.1080/15481603.2017.1396658 (2018).Article 

Google Scholar 
McPherson, E. G., Simpsona, J. R., Xiao, Q. & Wu, C. Million trees Los Angeles canopy cover and benefit assessment. Landsc. Urban. Plan. 99, 40–50 (2011).Article 

Google Scholar 
Baines, O., Wilkes, P. & Disney, M. Quantifying urban forest structure with open-access remote sensing data sets. Urban For. Urban Green. 50, 126653. https://doi.org/10.1016/j.ufug.2020.126653 (2020).Article 

Google Scholar 
Nowak, D. J. et al. Measuring and analyzing urban tree cover. Landsc. Urban Plan. 36, 49–57 (1996).Article 

Google Scholar 
Estoque, R. C. et al. Remotely sensed tree canopy cover-based indicators for monitoring global sustainability and environmental initiatives. Environ. Res. Lett. 16, 044047. https://doi.org/10.1088/1748-9326/abe5d9 (2021).ADS 
CAS 
Article 

Google Scholar 
Paap, T., de Beer, W., Migliorini, D., Nel, W. J. & Wingfield, M. J. The polyphagous shot hole borer (PSHB) and its fungal symbiont Fusarium euwallaceae: A new invasion in South Africa Trudy. Aust. Plant. Pathol. 47, 231–237. https://doi.org/10.1007/s13313-018-0545-0 (2018).Article 

Google Scholar  More

Cox, P. A., Banack, S. A. & Murch, S. J. Biomagnification of cyanobacterial neurotoxins and neurodegenerative disease among the Chamorro people of Guam. Proc. Natl. Acad. Sci. U.S.A. 100, 13380–13383 (2003).ADS 
CAS 
Article 

Google Scholar 
Brand, L. E., Pablo, J., Compton, A., Hammerschlag, N. & Mash, D. C. Cyanobacterial blooms and the occurrence of the neurotoxin, beta-N-methylamino-L-alanine (BMAA), in south Florida aquatic food webs. Harmful Algae 9, 620–635 (2010).CAS 
Article 

Google Scholar 
Metcalf, J. S., Banack, S. A., Richer, R. & Cox, P. A. Neurotoxic amino acids and their isomers in desert environments. J. Arid Environ. 112, 140–144 (2015).ADS 
Article 

Google Scholar 
Violi, J. P., Mitrovic, S. M., Colville, A., Main, B. J. & Rodgers, K. J. Prevalence of (beta)-methylamino-L-alanine (BMAA) and its isomers in freshwater cyanobacteria isolated from eastern Australia. Ecotoxicol. Environ. Saf. 172, 72–81 (2019).CAS 
Article 

Google Scholar 
Jonasson, S. et al. Transfer of a cyanobacterial neurotoxin within a temperate aquatic ecosystem suggests pathways for human exposure. Proc. Natl. Acad. Sci. 107, 9252–9257 (2010).ADS 
CAS 
Article 

Google Scholar 
Metcalf, J. et al. Toxin analysis of freshwater cyanobacterial and marine harmful algal blooms on the west coast of Florida and implications for estuarine environments. Neurotox. Res. 39, 27–35 (2021).CAS 
Article 

Google Scholar 
Cox, P. A. et al. Cyanobacteria and BMAA exposure from desert dust: a possible link to sporadic ALS among Gulf War veterans. Amyotroph. Lateral Scler. 10, 109–117 (2009).CAS 
Article 

Google Scholar 
Charlton, T. S., Marini, A. M., Markey, S. P., Norstog, K. & Duncan, M. W. Quantification of the neurotoxin 2-amino-3-(methylamino)-propanoic acid (BMAA) in Cycadales. Phytochemistry 31, 3429–3432 (1992).CAS 
Article 

Google Scholar 
Whiting, M. G. Toxicity of cycads. Econ. Bot. 17, 270–302 (1963).Article 

Google Scholar 
Cox, P. A., Davis, D. A., Mash, D. C., Metcalf, J. S. & Banack, S. A. Dietary exposure to an environmental toxin triggers neurofibrillary tangles and amyloid deposits in the brain. Proc. R. Soc. B: Biol. Sci. 283, 20152397 (2016).Article 

Google Scholar 
Scott, L. L. & Downing, T. G. A single neonatal exposure to BMAA in a rat model produces neuropathology consistent with neurodegenerative diseases. Toxins 10, 22 (2018).Article 

Google Scholar 
Roy, U. et al. Metabolic profiling of zebrafish (Danio rerio) embryos by NMR spectroscopy reveals multifaceted toxicity of (beta)-methylamino-L-alanine (BMAA). Sci. Rep. 7, 1–12 (2017).ADS 
Article 

Google Scholar 
Purdie, E. L., Metcalf, J. S., Kashmiri, S. & Codd, G. A. Toxicity of the cyanobacterial neurotoxin (beta)-N-methylamino-L-alanine to three aquatic animal species. Amyotroph. Lateral Scler. 10, 67–70 (2009).CAS 
Article 

Google Scholar 
Brenner, E. D. et al. Arabidopsis mutants resistant to s (+)-(beta)-methyl-(alpha), (beta)-diaminopropionic acid, a cycad-derived glutamate receptor agonist. Plant Physiol. 124, 1615–1624 (2000).CAS 
Article 

Google Scholar 
Schneider, D., Wink, M., Sporer, F. & Lounibos, P. Cycads: Their evolution, toxins, herbivores and insect pollinators. Naturwissenschaften 89, 281–294 (2002).ADS 
CAS 
Article 

Google Scholar 
Koi, S. & Daniels, J. Life history variations and seasonal polyphenism in Eumaeus atala (Lepidoptera: Lycaenidae). Florida Entomol. 100, 219–229 (2017).Article 

Google Scholar 
Koi, S. A butterfly picks its poison: Cycads (Cycadaceae), integrated pest management and Eumaeus atala Poey (Lepidoptera: Lycaenidae). Entomol. Ornithol. Herpetol. 6 (2017).Brenner, E. D., Stevenson, D. W. & Twigg, R. W. Cycads: Evolutionary innovations and the role of plant-derived neurotoxins. Trends Plant Sci. 8, 446–452 (2003).CAS 
Article 

Google Scholar 
Prado, A. The cycad herbivores. Bull. Soc. D’entomol. Quebec 18, 3–6 (2011).
Google Scholar 
Popova, A. & Koksharova, O. Neurotoxic non-proteinogenic amino acid (beta)-N-methylamino-L-alanine and its role in biological systems. Biochem. Mosc. 81, 794–805 (2016).CAS 
Article 

Google Scholar 
Salzman, S., Whitaker, M. R. L. & Pierce, N. E. Cycad-feeding insects share a core gut microbiome. Biol. J. Lin. Soc. 123, 728–738 (2018).Article 

Google Scholar 
Whitaker, M. R. & Salzman, S. Ecology and evolution of cycad-feeding Lepidoptera. Ecol. Lett. 23, 1862–1877 (2020).Article 

Google Scholar 
Zhou, X., Escala, W., Papapetropoulos, S., Bradley, W. G. & Zhai, R. G. BMAA neurotoxicity in Drosophila. Amyotroph. Lateral Scler. 10, 61–66 (2009).CAS 
Article 

Google Scholar 
Zhou, X., Escala, W., Papapetropoulos, S. & Zhai, R. G. (beta)-N-methylamino-L-alanine induces neurological deficits and shortened life span in Drosophila. Toxins 2, 2663–2679 (2010).CAS 
Article 

Google Scholar 
Mekdara, N. T. et al. A novel lenticular arena to quantify locomotor competence in walking fruit flies. J. Exp. Zool. A Ecol. Genet. Physiol. 317, 382–394 (2012).Article 

Google Scholar 
Goto, J. J., Koenig, J. H. & Ikeda, K. The physiological effect of ingested (beta)-N-methylamino-L-alanine on a glutamatergic synapse in an in vivo preparation. Comp. Biochem. Physiol. Part C: Toxicol. Pharmacol. 156, 171–177 (2012).CAS 

Google Scholar 
Okle, O., Rath, L., Galizia, C. G. & Dietrich, D. R. The cyanobacterial neurotoxin (beta)-N-methylamino-L-alanine (BMAA) induces neuronal and behavioral changes in honeybees. Toxicol. Appl. Pharmacol. 270, 9–15 (2013).CAS 
Article 

Google Scholar 
Spencer, P. S. et al. Guam amyotrophis lateral sclerosis-parkinsonism-dementia linked to a plant excitant neurotoxin. Science 237, 517–522 (1987).ADS 
CAS 
Article 

Google Scholar 
Bernays, E. A. & Chapman, R. F. Host-plant selection by phytophagous insects. In Host-Plant Selection by Phytophagous Insects. Contemporary Topics in Entomology, vol. 2, 201–213 (Springer, Boston, MA, 1994).Zandt, P. A. V. Plant defense, growth, and habitat: A comparative assessment of constitutive and induced resistance. Ecology 88, 1984–1993 (2007).Article 

Google Scholar 
Duncan, M. W. Role of the cycad neurotoxin BMAA in the amyotrophic lateral sclerosi-parkisonism dementia complex of the Western Pacific. Adv. Neurol. 56, 301–310 (1991).CAS 
PubMed 

Google Scholar 
Banack, S. A. & Cox, P. A. Distribution of the neurotoxic nonprotein amino acid BMAA in Cycas micronesica. Bot. J. Linn. Soc. 143, 165–168 (2003).Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria (2021).Therneau, T. M. A Package for Survival Analysis in R. R package version 3.2-11 (2021).Kassambara, A., Kosinski, M. & Biecek, P. survminer: Drawing Survival Curves using ’ggplot2’. R package version 0.4.9 (2021).Pennington, Z. T. et al. eztrack: An open-source video analysis pipeline for the investigation of animal behavior. Sci. Rep. 9, 1–11 (2019).Article 

Google Scholar 
Pérez, F. & Granger, B. E. IPython: A system for interactive scientific computing. Comput. Sci. Eng. 9, 21–29 (2007).Article 

Google Scholar 
Hammer, T. J., Janzen, D. H., Hallwachs, W., Jaffe, S. P. & Fierer, N. Caterpillars lack a resident gut microbiome. Proc. Natl. Acad. Sci. 114, 9641–9646 (2017).CAS 
Article 

Google Scholar 
Karlsson, O., Roman, E. & Brittebo, E. B. Long-term cognitive impairments in adult rats treated neonatally with (beta)-N-methylamino-L-alanine. Toxicol. Sci. 112, 185–195 (2009).CAS 
Article 

Google Scholar 
Whitaker, M. R. L., Salzman, S., Gratacos, X. & Tucker Lima, J. Localized overabundance of an otherwise rare butterfly threatens endangered cycads. Florida Entomol. 103, 519–522 (2021).Article 

Google Scholar 
Backmann, P. et al. Delayed chemical defense: Timely expulsion of herbivores can reduce competition with neighboring plants. Am. Nat. 193, 125–139 (2019).Article 

Google Scholar 
Yáñez-Espinosa, L. & Sosa-Sosa, F. Population structure of Dioon purpusii rose in Oaxaca, Mexico. Neotrop. Biol. Conserv. 2, 46–54 (2007).
Google Scholar 
Robbins, R. K. et al. A switch to feeding on cycads generates parallel accelerated evolution of toxin tolerance in two clades of Eumaeus caterpillars (Lepidoptera: Lycaenidae). Proc. Natl. Acad. Sci.118 (2021).Grunseich, J. M., Thompson, M. N., Aguirre, N. M. & Helms, A. M. The role of plant-associated microbes in mediating host-plant selection by insect herbivores. Plants 9, 6 (2020).CAS 
Article 

Google Scholar 
Zhang, Y. & Whalen, J. K. Production of the neurotoxin beta-N-methylamino-L-alanine may be triggered by agricultural nutrients: An emerging public health issue. Water Res. 170, 115335 (2020).CAS 
Article 

Google Scholar  More

FennelThe main effect of fertilization (F) significantly impacted all measured parameters of fennel. Intercropping (I) pattern affected all parameters except plant height and 1000-seed weight. Significant I × F interactions occurred for umbel number, seed yield, essential oil content (EO), EO yield, oil content, and oil yield (Table 2).Table 2 Analysis of variance for the effect of cropping pattern and fertilization on evaluated traits in fennel.Full size tablePlant heightThe tallest fennel plants (125.8 cm) occurred in the sole cropping (Fs), while the shortest plants (98.6 cm) occurred in 1F:2FG. Across intercropping patterns, the Fs treatment had 28%, 14%, 11% taller fennel plants than 1F:2FG, 2F:2FG, and 2F:4FG, respectively (Fig. 1A). Compared to the unfertilized control, HA and BFS increased fennel plant height by 10% and 13%, respectively (Fig. 1B).Figure 1Means comparison for the main effects of cropping patterns [Fs (fennel sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on plant height (A), and fertilization [C (control), HA (humic acid), BFS (biofertilizers)] on plant height (B) and 1000-seed weight (C) of fennel. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size image1000-seed weightCompared to the unfertilized control (3.9 g per 1000 seeds), BFS and HA increased the 1000-seed weight of fennel by 24.1% and 14.5% (4.9 and 4.5 g per 1000 seeds), respectively (Fig. 1C).Umbel numberThe fennel sole cropping fertilized with HA produced the most umbels of fennel (51.5), while 2F:4FG without fertilization produced the least (32). Averaged across fertilizer types within each intercropping system, 1F:2FG, 2F:2FG, and 2F:4FG had 21.1%, 16.3%, and 26.7% fewer umbels than fennel sole cropping, respectively. Across intercropping systems, HA and BFS increased the umbel number by 17.8% and 16.5% compared with the unfertilized control, respectively (Fig. 2A).Figure 2Means comparison for the interaction effect of fertilization [C (control), HA (humic acid), BFS (biofertilizers)] and different cropping patterns [Fs (fennel sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on umbel number (A) and seed yield (B) of fennel. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size imageSeed yieldThe different intercropping patterns had lower fennel seed yields than fennel sole cropping. Sole cropping fertilized with BFS and HA produced the highest fennel seed yields (2233 and 2209 kg ha–1, respectively), followed by unfertilized sole cropping (1960 kg ha–1). The lowest seed yields occurred in the unfertilized controls in 1F:2FG (933 kg ha–1) and 2F:4FG (1033 kg ha–1). Averaged across fertilization treatments, fennel seed yield in 1F:2FG, 2F:2FG, and 2F:4FG decreased by 41.7, 26.8, and 36.3%, respectively, compared to fennel sole cropping (Fs). Averaged across intercropping patterns, HA and BFS increased fennel seed yield by 33.3% and 39.5% compared with the unfertilized control, respectively (Fig. 2B).Essential oil content and yieldThe different intercropping patterns produced higher EO contents of fennel than fennel sole cropping. The highest absolute EO content of fennel (4.22%) occurred in 2F:2FG fertilized with BFS, although this did not statistically differ from the 2F:2FG fertilized with HA (4.04%) or 2F:4FG fertilized with HA or BFS (3.8% and 4.00%, respectively) (Fig. 3A). The lowest EO contents occurred in the unfertilized control (2.38%), HA (2.55%), and BFS (2.57%) in the Fs system. Averaged across fertilization treatments, the EO content of fennel in 1F:2FG, 2F:2FG, and 2F:4FG increased by 36%, 52%, and 44% compared to fennel sole cropping, respectively. Within each intercropping pattern, and with the exception of Fs, the HA and BFS treatments had higher EO contents of fennel, none of which significantly differed, increasing by 25% and 29%, respectively (Fig. 3A).Figure 3Means comparison for the interaction effect of fertilization [C (control), HA (humic acid), BFS (biofertilizers)] and different cropping patterns [Fs (fennel sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on essential oil content (A), essential oil yield (B), oil content (C), and oil yield (D) of fennel. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size imageMaximum EO yields of fennel occurred with HA or BFS applied in 2F:2FG (65.2 and 66.6 kg ha–1) and 2F:4FG (60.7 and 65.5 kg ha–1), respectively, while the lowest EO yields occurred in the unfertilized control in 1F:2FG (27.2 kg ha–1), 2F:2FG (33.2 kg ha–1), and 2F:4FG (32.1 kg ha–1). Averaged across intercropping patterns, the EO yield of fennel increased by 66.1% and 74.7% with HA and BFS, respectively (Fig. 3B).Fennel essential oil compositionGC–FID and GC–MS analyses identified 14 components in the fennel EO (representing 97.4–99.9% of the total composition) (Table 3), with the main constituents being trans-anethole (78.3–84.85%), estragole (3.02–7.17%), fenchone (4.14–7.52%), and limonene (3.15–4.88%). The highest percentage of (E)-anethole, estragole, and fenchone occurred in 2F:2FG with BFS. The highest limonene content occurred in 2F:4 FG with HA. The relative contents of trans-anethole, fenchone, and limonene increased by 3.9%, 16.6%, and 8.4% compared with fennel sole cropping. Notably, the contents of most compounds increased with HA and BFS. Compared to the unfertilized control, trans-anethole, fenchone, and limonene contents increased by 2.9%, 21.5%, and 7.9% with BFS and 2.3%, 22.4%, and 11.9% with HA, respectively (Table 3).Table 3 Proportion of fennel essential oil constituents under different cropping patterns and fertilization.Full size tableFennel oil content and yieldAmong the studied treatments, the highest fennel oil content occurred with HA or BFS application in 1F:2FG (16.3% and 16.6%) and 2F:2FG (16.3% and 17.4%), respectively. The lowest fennel oil contents occurred in the unfertilized control, HA, and BFS treatments (12.5%, 12.8%, and 12.9%, respectively) under fennel sole cropping, and the unfertilized control in 2F:4FG (12.6%). Averaged across fertilizer treatments, fennel oil content in 1F:2FG, 2F:2FG, and 2F:4FG increased by 22.8%, 26.0%, and 12.6% compared with fennel sole cropping, respectively. Across intercropping patterns, HA and BFS increased fennel oil content by 13.5% and 16.5%, respectively (Fig. 3C).The maximum oil yield of fennel (318.6 kg ha–1) occurred in 2F:2FG fertilized with BFS, while the lowest oil yield (129.3 kg ha–1) occurred in 1F:2FG without fertilization. Across intercropping patterns, HA and BFS increased fennel oil yield by 50.8% and 62.6%, respectively (Fig. 3D).Oil compoundsGC–FID and GC–MS analyses identified nine constituents that represented 94.3–97.9% of the total fennel oil composition. The main oil constituents were oleic acid (39.2–48.3%), linoleic acid (17.1–24.8%), stearic acid (10.9–15.4%), lauric acid (10.1–14.00%), and arachidic acid (2.2–3.4%). The highest oleic and linoleic acid contents occurred in 2F:4FG and 2F:2FG fertilized with BFS, respectively. Across fertilizer treatments, oleic and linoleic acid contents increased by 6% and 21%, respectively, under different intercropping patterns compared with fennel sole cropping. Across systems, HA and BFS enhanced oleic acid content by 1.8% and 8% and linoleic acid by 7.9% and 8.2%, respectively, compared with the unfertilized control. The highest percentage of stearic and lauric acids occurred in the unfertilized control of fennel sole cropping. Conversely, the lowest stearic and lauric acid contents occurred in 2F:2FG and 2F:4FG fertilized with BFS, 16.1% and 14.2% higher than fennel sole cropping, respectively. Finally, HA and BFS decreased stearic acid content by an average of 5.4% and 7.2%, respectively (Table 4).Table 4 Proportion of fennel oil constituents under different cropping patterns and fertilization.Full size tablePhenolic compoundsThe main phenolic compounds of fennel were chlorogenic acid (10.4–15.3 ppm), quercetin (7.0–17.2 ppm), and cinnamic acid (4.1–8.9 ppm). The highest chlorogenic acid and quercetin contents occurred in 2F:2FG fertilized with BFS and HA, respectively, while the lowest contents occurred in the fennel sole cropping system without fertilizer. Averaged across the three intercropping patterns, the chlorogenic acid and quercetin contents were 18.5% and 80.1% higher than the fennel sole cropping system. The chlorogenic acid and quercetin contents increased by 13% and 17% with BFS and 22% and 15% with HA, respectively (Table 5).Table 5 Concentration of phenolic compounds in fennel under different cropping patterns and fertilization.Full size tableFenugreekThe main effects of intercropping (I) pattern (C) and fertilizer (F) were significant for all parameters analyzed in fenugreek. Significant I × F interactions occurred for plant height, pod number per plant, seed yield, oil content, and oil yield of fenugreek (Table 6).Table 6 Analysis of variance for the effects of cropping patterns and fertilization on evaluated traits in fenugreek.Full size tablePlant heightThe 2F:2FG intercropping system fertilized with BFS produced the tallest fenugreek plants (63 cm), followed by 1F:2FG with BFS (53.3 cm) and 2F:4FG with BFS (56 cm), and 2F:2FG with HA (55 cm). The unfertilized control produced the shortest fenugreek plants (42 cm) in the sole cropping. Most fertilizer treatments across different intercropping patterns produced taller fenugreek plants than their sole cropping counterparts. Across fertilizer treatments, 1F:2FG, 2F:2FG, and 2F:4FG produced 16.2%, 26.8%, and 14.6% taller fenugreek plants than sole cropping, respectively. Across cropping patterns, BFS and HA increased fenugreek plant height by 5.7% and 15.2% compared with the unfertilized control, respectively (Fig. 4A).Figure 4Means comparison for the interaction effect of fertilization [C (control), HA (humic acid), BFS (biofertilizers)] and different cropping patterns [FGs (fenugreek sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on plant height (A) and pod number per plant (B) of fenugreek. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size imagePod number per plantThe fenugreek sole cropping with BFS and HA and 2F:4FG with BFS produced the most pods per plant (21.3, 20.3, and 20, respectively), while the unfertilized controls in 1F:2FG, 2F:2FG, and 2F:4FG produced the least (11.6, 12, and 13.3, respectively). Across fertilization treatments, 1F:2FG, 2F:2FG, and 2F:4FG had 30.1%, 25.6%, and 14.3% fewer pods per plant, respectively, than the fenugreek sole cropping system. Across cropping systems, HA and BFS increased pod number per plant in fenugreek by 25% and 33%, respectively, relative to the corresponding sole cropping (Fig. 4B).Seed number per podAcross fertilization treatments, fenugreek sole cropping produced the most seeds per pod (7.09), followed by 2F:4FG (6.02), 2F:2FG (4.93), and 1F:2FG (4.41) (Fig. 5A). In relative terms, sole cropping produced 60.5%, 43.9%, and 17.6% more seeds per pod than 1F:2FG, 2F:2FG, and 2F:4FG (Fig. 5A). Across cropping patterns, BFS and HA increased seed number per pod by 8.1% and 17.4% compared with the unfertilized control, respectively (Fig. 5B).Figure 5Means comparison for the main effects of cropping patterns [FGs (fenugreek sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on seed number per pod (A) and 1000-seed weight (C), and fertilization [C (control), HA (humic acid), BFS (biofertilizers)] on seed number per pod (B) and 1000-seed weight (D) of fennel. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size image1000-seed weightAmong different cropping patterns, sole cropping and 1F:2FG produced the highest (10.45 g) and lowest (8.34 g) fenugreek seed weights, respectively. In relative terms, fenugreek sole cropping produced 25.3%, 21.8%, and 12.4% higher seed weights than 1F:2FG, 2F:2FG, and 2F:4FG, respectively (Fig. 5C). Across cropping patterns, BFS and HA increased fenugreek seed weight by 3.7% and 5.7% compared with the control, respectively (Fig. 5D).Seed yieldMeans comparisons showed that sole cropping produced higher fenugreek seed yields than intercropping patterns. Sole cropping with BFS (1240 kg ha–1) and HA (1217 kg ha–1) produced the highest seed yields followed by the unfertilized control (Fig. 6A). The unfertilized control in 1F:2FG (437 kg ha–1) and 2F:2FG (467 kg ha–1) produced the lowest fenugreek seed yields. In all cases, and within each cropping pattern, BFS and HS produced higher fenugreek seed yields than the unfertilized control. As a result, BFS and HA increased fenugreek seed yield by 25.2% and 31.5% compared with the unfertilized control, respectively (Fig. 6A).Figure 6Means comparison for the interaction effects of fertilization [C (control), HA (humic acid), BFS (biofertilizers)] and different cropping patterns [FGs (fenugreek sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on seed yield (A), oil content (B), and oil yield (C) of fenugreek. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size imageOil content and yieldThe 2F:2FG cropping pattern with BFS produced the highest fenugreek oil content (8.3%), while the unfertilized control in sole cropping produced the lowest (5.9%). Across fertilizer treatments, 1F:2FG, 2F:2 FG, and 2F:4 FG produced 11.7%, 18.5%, and 15.7% higher fenugreek oil contents than sole cropping, respectively. In the 2F:2FG and 2F:4FG cropping patterns, BFS produced higher oil content (%) than HA. As a result, across cropping patterns, HA and BFS increased fenugreek oil content by 12.3% and 19.4%, respectively (Fig. 6B).Sole cropping with HA and BFS and 2F:2FG with BFS produced the highest fenugreek oil yields (77.1, 80.0, and 74.4 kg ha–1, respectively), while the unfertilized controls in 1F:2FG and 2F:4FG produced the lowest (27.51 and 29.8 kg ha–1, respectively). The 1F:2FG, 2F:2FG, and 2F:4FG cropping patterns produced 45.9%, 20.7%, and 41.5% lower fenugreek oil yields than fenugreek sole cropping, respectively. Moreover, except for sole cropping, BFS produced the highest fenugreek oil yield, followed by HA and the unfertilized control (Fig. 6C).Oil compoundsGC–FID and GC–MS analyses identified seven constituents (representing 91.09–99.27% of the total composition) in fenugreek oil. The main oil constituents were linoleic acid (26.1–37.1%), linolenic acid (16.9–22.4%), oleic acid (15.1–21.2%), palmitic acid (11.2–17.1%), lauric acid (5.0–12.3%), and myristic acid (3.1–6.4%). The highest linoleic and oleic acid percentages occurred in 1F:2FG and 2F:4FG with BFS. The 1F:2FG cropping pattern with BFS also had the highest linolenic acid percentage. The fenugreek sole cropping system without fertilization (control) had the lowest content of these three compounds. The intercropping patterns had 17%, 18.2%, and 17.1% higher oleic, linoleic, and linolenic acid contents than fenugreek sole cropping. In addition, HA and BFS increased oleic acid content by 15.6% and 8.8%, linoleic acid content by 12.8% and 7%, and linolenic acid content by 7.5% and 12.9%, respectively. Fenugreek sole cropping without fertilization produced the highest lauric acid and palmitic contents, 29.33% and 22.81% higher than the intercropping patterns (Table 7).Table 7 Proportion of fenugreek oil constituents under different cropping patterns and fertilization.Full size tablePhenolic compoundsThe main phenolic compounds in fenugreek were chlorogenic acid (2.01–5.49 ppm), caffeic acid (2.42–4.93 ppm), quercetin (1.98–4.45 ppm), comaric (1.09–2.43 ppm), apigenin (1.97–2.99 ppm), and gallic acid (1.76–2.92 ppm). The 2F:2FG cropping pattern with HA produced the highest quercetin and gallic acid contents, and 2F:4FG with HA produced the highest chlorogenic and caffeic acid contents. The 2F:2FG and 2F:4FG cropping patterns with BFS produced the highest comaric and apigenin contents, respectively. In contrast, fenugreek sole cropping without fertilization produced the lowest contents of the abovementioned compounds (Table 8).Table 8 Proportion of fenugreek concentration of phenolic compounds under different cropping patterns and fertilization.Full size tableLand equivalent ratio (LER)The 2F:4FG and 2F:2FG intercropping patterns treated with BFS had the highest partial LERs for fennel (0.82) and fenugreek (0.72), respectively. In addition, 2F:2FG with BFS and 1F:2FG without fertilization produced the highest (1.42) and lowest (0.86) total LERs, respectively (Fig. 7).Figure 7Partial and total land equivalent ratio (LER) for seed yields of different fennel and fenugreek intercropping patterns [1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] and fertilization [C (Control), HA (humic acid), BFS (biofertilizers)]. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size image More

Experimental design and crop managementThe study was conducted during 2019 in a field experiment on an arable soil (classified as Luvisols) in the deep root experimental facility at the University of Copenhagen, Denmark (Supplementary Table S4). The experiment was conducted with two diverse perennial deep-rooted species: the tap-rooted forage legume lucerne (Medicago sativa L. (cv. Creno); Family: Fabaceae) with the capacity to fix N2 and the intermediate wheatgrass (Thinopyrum intermedium; Family: Poaceae) kernza developed by the Land Institute (Salina, Kansas, USA). Kernza was initially sown on April 11th, 2015 and lucerne on September 9th, 2016 with a seeding density of 20 kg seeds ha−1. Every year, kernza was fertilized with NPK fertilizer (21:7:3; NH4:NO3 = 1.28) as a single dose in early spring (before the onset of plant growth). Kernza was harvested every year in August using a combine harvester and lucerne three times per year in June, August, and October. Plants were rainfed with a subsurface drain installed at both 1 and 2 m depth running between the plots.For each species, fixed frames of 0.75 m2 were inserted in the soil (ca. 5 cm) within each field plot. Specifically, three field plots of lucerne (with observable root nodulation) and kernza were used where each of the three kernza field plots contained two subplots of N fertilized kernza at 100 kg N ha−1 (K100) (i.e., the standard fertilization within this field) and N fertilized kernza at 200 kg N ha−1 (K200) (i.e., within the range of standard fertilization practices for kernza). Before the onset of plant growth, all plots received 15N (as 15NH4Cl; 98 atom%) in trace amounts (corresponding to 1 kg N ha−1) to trace N allocation from the surface to deeper layers.
13C/14C-CO2-labelingWithin each fixed frame, the 13C/14C-CO2-labeling was conducted using an atmospheric labeling chamber41. Labeling with C-tracers was done with multiple-pulse labeling (three times per week) over two months until first harvest (May 2nd to June 20th 2019). Glass beakers containing 13C labeled bicarbonate (0.1 g mL−1 labeling solution; 99 atom%), and 14C labeled bicarbonate (11 kBq mL−1 labeling solution) within a solution of NaOH (1 M) were added within each of the labeling chambers. Once chambers were sealed, hydrochloric acid (HCl; 2 M) was added to the labeling solution (in equivalent amounts) via a syringe promoting 14CO2/13CO2 evolution. Chambers remained sealed for one to three hours (between 9 am and 12 pm) depending on weather conditions (i.e., the duration and intensity of sunshine). The amount of added labeling solution sequentially increased with increasing plant growth (i.e., 5 mL per 20 cm increase in plant height) reaching a plant height of 100–120 cm at the termination of the labeling.Shoot, root, and soil samplingThe labeling plots (0.75 m2) were harvested on June 20th, 2019 to obtain the aboveground biomass of lucerne and kernza (K100 and K200). The aboveground biomass in addition to samples obtained from unlabeled parts of the field was directly stored at − 20 °C until drying at 105 °C for two days. For each plot and unlabeled samples, the plant biomass was homogenized and ball-milled for subsequent isotopic analyses.Soil cores to 1.5 m depth were taken inside all labeling plots, and cores were subdivided into four depth intervals: 0–25, 25–50, 50–100, and 100–150 cm. The soil coring was conducted in 25 cm intervals using a soil auger (6 cm inner diameter). Specifically, per depth three soil samples were taken and stored at 4–5 °C (ca. two days) and then immediately processed and stored at -20 °C until analyses. Roots, bulk soil and rhizosphere soil (adhering to the roots), were separated by sequential sieving of the soil with finer mesh sizes to 1 mm as described by Peixoto, et al.26. A subsample of the bulk soil (ca. 150 g) from each depth in all labeling plots was washed on a 250 µm sieve to recover root fragments for subsequent isotopic determination in unrecovered root fragments. Soil samples (and associated roots) from unlabeled parts of the larger field plots were used to determine natural abundance of 13C/14C/15N with depth. The collection of plant material complied with relevant institutional guidelines and seeds were gifted by University of Copenhagen.Determination of 13C/14C/15N enrichment, and C and N quantityFor each defined depth, samples of roots and soil were homogenized, freeze-dried (except PLFA samples that were stored at − 20 °C), and ground in a ball-mill for the determination of total C and N, 13C, 15N, and 14C activity. Total C, N, 13C, and 15N were measured with a FLASH 2000 CHNS/O Elemental Analyzer (Thermo Fisher Scientific, Cambridge, UK) combined to a Delta V Advantage isotope ratio mass spectrometer via a ConFlo III interface (Thermo Fisher Scientific, Bremen, Germany) at the Centre for Stable Isotope Research and Analysis (Georg August University Göttingen, Göttingen, Germany).All δ13C values are standardized to the Vienna PeeDee Belemnite international isotope standard and δ15N values standardized to the δ15N values of atmospheric N2. 13C and 15N enrichment is expressed as atom% excess as calculated by the atom% difference between the respective labeled and unlabeled samples. The 14C activity was determined by combustion in a Hidex 600 OX Oxidizer (Hidex, Turku, Finland) and counted on a liquid scintillation counter (Tri-Carb 3180TR/SL, PerkinElmer, Waltham, MA, USA). 14C enrichment is determined by the difference in the 14C activity (Bq g−1) between the respective labeled and unlabeled samples.Calculation of root C and net rhizodepositionThe amount of root C (mg C kg−1 soil) was calculated based on the root dry matter and C concentration divided by the quantity of soil sampled38. For the determination of net rhizodeposition, 14C was used due to lower detection limits in deeper soil layers42. A modified tracer mass balance approach described by Rasmussen, et al.43 with adjusted unrecovered root fragments41 was used to determine the net rhizodeposition based on the following equations where the %ClvR is the relative proportion of rhizodeposition expressed as the percent C lost via rhizodeposition:$${text{%ClvR}} = frac{{^{{{14}}} {text{C Soil (rhizosphere + adjusted bulk)}}}}{{^{{{14}}} {text{C bulk soil }} + ,^{{{14}}} {text{C rhizosphere soil}} + ,^{{{14}}} {text{C Root}}}} times 100.$$$${text{Net rhizodeposition}} = frac{{{text{%ClvR }} times {text{ root C content}}}}{{left( {100 – % {text{ClvR}}} right)}}$$The 14C soil content was the sum of the adjusted bulk soil 14C and rhizosphere 14C content for each soil sample. The 14C rhizosphere and bulk soil content for each soil sample were determined by multiplying the total quantity of C by the 14C enrichment of the soil. The adjusted bulk soil 14C content was calculated as the difference between the bulk 14C soil content by the 14C root washed content as determined by the multiplication of 14C enrichment in root fragments recovered from a subsample of soil by the total C content within the entire soil volume sampled. The 14C root content was determined by multiplying the total quantify of C in roots by the 14C enrichment. Similar equations were used to calculate the net rhizodeposition of N based on 15N enrichment within the soil and roots.Biomarker analysesPhospholipid fatty acid (PLFA)The analysis of PLFAs was done according to a modified protocol by Frostegård, et al.44 with a detailed description of the modifications provided by Gunina, et al.45. In brief, 25 μL of 1,2-Dinonadecanoyl-sn-Glycero-3-Phosphatidylcholine (C19:0) (1 mg mL–1) were added to each of the samples and used in the quantification of recovery of the phospholipids. The lipid fraction from 5–6 g of rhizosphere soil was extracted twice using a one-phase Bligh-Dyer extractant46 of chloroform, methanol (MeOH), and citrate buffer (pH 4) (1:2:0.8, v/v/v). To isolate the phospholipid fraction, a solid-phase extraction with activated silica gel and methanol elution was conducted. The derivatization into fatty acid methyl esters occurred via a sequential hydrolyzation with 0.5 mL sodium hydroxide (NaOH) (0.5 M) in MeOH for 10 min at 100 °C and methylation with 0.75 mL of boron trifluoride (BF3) (1.3 M) in MeOH for 15 min at 80 °C. An external standard stock solution containing 28 individual fatty acids (ca. 1 mg mL–1 per fatty acid) used in the quantification of PLFA content was simultaneously derivatized with the samples. The residues were dissolved in 185 μL of toluene, and 15 μL of the internal standard 2, tridecanoic acid methyl ester (C13:0) (1 mg mL–1) were added to each sample prior to measurement using an Agilent 7820A GC coupled to an Agilent 5977 quadrupole mass spectrometer (Agilent, Waldbronn, Germany). The sum of all PLFAs was used as a proxy of the living microbial biomass based on the direct relation between PLFAs and microbial biomass.Amino sugars (AS)Amino sugars were extracted according to a modified protocol by Zhang and Amelung47 with a detailed description of the procedure by Peixoto, et al.26. In brief, 0.8–1.5 g of freeze-dried rhizosphere soil were hydrolyzed with the addition of 11 mL of 6 M HCl for 8 h at 105 °C. Following hydrolysis, soil samples were filtered and HCl was removed via rotary evaporation at 45 °C to dry the filtrate. Prior to derivatization both iron precipitates and salts were removed from the filtrate and 25 μL of the internal standard 1, methylglucamine (MeGlcN) (1 mg mL–1) was added and used for quantification of recovery. The derivatization into aldononitrile acetates was conducted as described by Zhang and Amelung47. For the quantification of AS, an external standard stock solution containing the AS: N-acetylglucosamine (GlcN) (2 mg mL–1), N-acetylgalactosamine (GalN) (2 mg mL–1), N-acetylmuramic acid (MurN) (1 mg mL–1), mannosamine (ManN) (2 mg mL–1), and MeGlcN (1 mg mL–1) was derivatized and analyzed with the samples. The residues were dissolved in 185 μL of ethyl acetate-hexane (1:1, v/v), and 15 μL of the internal standard 2, tridecanoic acid methyl ester (1 mg mL–1), were added to the samples for measurement using an Agilent 7890A GC coupled to Agilent 7000A triple quadrupole mass spectrometer (Agilent, Waldbronn, Germany). Total amino sugars content was calculated as the summation of the four detected amino sugars: GlcN, MurN, GalN, and ManN.Amino acids (AA)Amino acids were extracted from both freeze-dried rhizosphere soil and root samples according to the protocol by Enggrob, et al.48. In brief, 0.8–3 g of rhizosphere soil and 0.02 g of root were hydrolyzed with the addition of 2 mL of 6 M HCl for 20 h at 110 °C to break the peptide bonds. Samples were subsequently purified via the removal of lipophilic and solid compounds by the addition of 4 mL n-hexane/dichloromethane (6:5, v/v) to the soil and root samples. Following centrifugation, the aqueous phase was filtered through glass wool and rinsed with 2 × 0.5 mL 0.1 M HCl into new glass tubes with the addition of 300 μL of the internal standard, norleucine (2.5 mM). The samples were freeze-dried and the residues dissolved in 1 mL 0.01 M HCl prior to the separation of amino acids and amino sugars (i.e., N containing compounds) on a polypropylene column with a cation exchange resin. The amino acids were eluted with a 2.5 M ammonium hydroxide solution and freeze-dried prior to derivatization of the amino acids as described by Enggrob, et al.48. For the quantification of AA, an external standard stock solution containing 14 AA was derivatized and analyzed with the samples. The amino acids were measured using a trace GC Ultra mounted with a TriPlus autosampler (Thermo Scientific, Hvidovre, Denmark) coupled via a combustion reactor (GC IsoLink, Thermo Scientific) to an isotope ratio mass spectrometer (Delta V Plus IRMS, Thermo Scientific). The total AA content of the rhizosphere soil and roots was based on the summation of the AA: alanine, Asx (asparagine and aspartate), Glx (glutamine and glutamate), glycine, isoleucine, lysine, phenylalanine, Pro/Thr (proline and threonine), serine, tyrosine, and valine.Compound-specific stable isotope probingTo determine the 13C enrichment of biomarkers, all raw δ13C were measured individually for AS and PLFA using a Delta V Advantage isotope ratio mass spectrometer via a ConFlo III interface (Thermo Fisher Scientific, Bremen, Germany). For AA, all raw δ13C were measured using a trace GC Ultra mounted with a TriPlus autosampler (Thermo Scientific, Hvidovre, Denmark) coupled via a combustion reactor (GC IsoLink, Thermo Scientific) to an isotope ratio mass spectrometer (Delta V Plus IRMS, Thermo Scientific). For each sample, chromatogram peaks identified based on retention times specific for the measured amino sugars, PLFA, and AA were integrated using Isodat v. 3.0 (Thermo Fisher Scientific). All raw δ13C values were corrected for dilution by additional C atoms added during the derivatization, amount dependence, offset, and drift (for PLFA samples)49,50,51. To determine the 13C incorporation into each biomarker, the 13C excess for each biomarker as determined by the difference between the 13C of the labeled and unlabeled biomarker was multiplied by the C content of the specific biomarker.Relative microbial stabilization (RMS)The relative microbial stabilization is based on the relation of rhizodeposited 13C in the PLFA and amino sugar pools as described in detail by Peixoto, et al.26. The underlying assumption is that 13C incorporation into the amino sugar pool indicates the transformation of rhizodeposited C into necromass52,53, and the 13C incorporation into the PLFA pool (i.e., the living microbial biomass) represents a temporary C pool as PLFAs are immediately exposed to degradation following cell lysis54. The relative microbial stabilization (RMS) is calculated as follows:$${text{Relative microbial stabilization}} = {text{log}}frac{{{text{Average weighted atom% }},^{{{13}}} {text{C excess AS}}}}{{{text{Average weighted atom% }},^{{{13}}} {text{C excess PLFA}}}}$$where the average weighted atom% 13C excess is determined by the total 13C incorporation divided by the total C content of the respective PLFA or amino sugar pools. Accordingly RMS  0 is indicative of higher stabilization of C based on the dominant entry of C into the microbial necromass. However, the RMS indicator does not imply the absolute stability of rhizodeposited C, but rather signifies the potential for microbial stabilization among contrasting experimental variables (i.e., depth and plant species).Molecular analysisDNA extractionFrom each sample, 0.5 g of freeze-dried rhizosphere soil was used for DNA extraction using the Fast DNA Spin kit for Soil (MP Biomedicals, Solon, OH, USA) according to the manufacturer’s protocol with a single modification. Following, the addition of Binding Matrix, the suspension was washed with 5.5 M Guanidine Thiocyanate (protocol from MP Biomedicals) to remove humic acids that could inhibit preceding polymerase chain reaction (PCR) steps. The DNA was eluted in DNase free water and purified using the NucleoSpin gDNA Clean-up kit following the manufacturer’s protocol (Macherey–Nagel, Düren, Germany). The purity and concentration of DNA were checked on Nanodrop and Qubit, respectively.Amplicon sequencingExtracted DNA was sent to Novogene Europe (Cambridge, United Kingdom) for library preparation and amplicon sequencing. For 16S rRNA gene amplicon sequencing of the V3-V4 regions, the primer pair 341 F and 806 R were used (Supplementary Table S5). To identify the fungal communities, we targeted the Internal Transcribed Spacer (ITS) Region 1, using the primer pair ITS1 and ITS2 (Supplementary Table S5). The constructed libraries were sequenced using a Novaseq 6000 platform producing 2 × 250 bp paired-end reads. Raw sequences were deposited in the NCBI Sequence Read Archive (Bioproject number PRJNA736561).Quantitative PCRCopy numbers of the 16S rRNA gene were determined by quantitative PCR (qPCR) using the primers 341F and 805R (Supplementary Table S5) on an AriaMX Real-Time PCR System (Agilent Technologies, Santa Clara, CA, USA). An external plasmid standard curve was made based on the pCR 2.1 TOPO vector (Thermo Fisher Scientific, Waltham, MA, USA) with a 16S rRNA gene insert amplified from bulk soil. The PCR reaction was performed in 20 µl reactions containing: 1 × Brilliant III Ultra-Fast SYBR green low ROX qPCR Master Mix (Agilent Technologies, Santa Clara, CA, USA), 0.05 µg/µl BSA (New England Biolabs Inc., Ipswich, MA, USA), 0.4 µM of each primer and 2 μl of template DNA. The thermal cycling conditions were 3 min at 95 °C followed by 40 cycles of 20 s at 95 °C and 30 s at 58 °C, and a final extension for 1 min at 95 °C. A melting curve was included according to the default settings of the AriaMx qPCR software (Agilent Technologies). The reaction efficiencies were between 97 and 102%. Fungal quantification was done by qPCR amplification of the Internal Transcribed Spacer 1 (ITS1) using the primers ITS1-F and ITS2 (Supplementary Table S5). A plasmid standard curve was made using the pCR 2.1 TOPO vector containing an ITS1 region from Penicillium aculeatum. Reaction mixture and cycling conditions were as described above for the 16S rRNA gene (Supplementary Table S5). The reaction efficiency was 84%.Quantification of functional genes involved in N cyclingThe five bacterial genes amoA, nirK, nirS, nosZ, and nifH coding for enzymes involved in N-cycling were quantified by qPCR on an AriaMx Real-Time PCR System (Agilent Technologies). Reaction mixtures and cycling conditions were as described above for the 16S rRNA gene (Supplementary Table S5). The standard curves were prepared as described in Garcia-Lemos, et al.55. The reaction efficiencies were in the range 87%-105%.Sequence processingRaw reads were treated using DADA2 version 1.14.156. In brief, reads were quality checked and primers were removed using Cutadapt v. 1.1557. We followed the protocol DADA2 using default parameters, with a few modifications. For 16S rRNA sequences, the forward and reverse reads were trimmed to 222 and 219 bp, respectively, while the maxEE was set to 2 and 5 for forward and reverse reads, respectively. Detection of amplicon sequence variants (ASVs) was done using the pseudo-pool option and forward and reverse reads were merged with a minimum overlap of 10 bp. Merged reads in the range of 395–439 bp were kept, as reads outside this range are considered too long or too short for the sequenced region. Taxonomy was assigned using the Ribosomal Database Project (RDP) classifier58 with the Silva database v.13859. For ITS region 1, quality filtered reads shorter than 50 bp were removed prior to merging the forward and the reverse reads, with maxEE set to two for both forward and reverse reads. During merging, the minimum overlap was set to 20 (default). Taxonomy was assigned with the RDP classifier using the Unite v. 8.2 database60 after removal of chimeras.As ITS region 1 has a variable length, reads can be lost during merging. Hence, to validate our dataset we ran only the forward reads through the DADA2 pipeline and compared the overall community structure with the dataset from the merging using a Mantel test. No significant changes were observed in the community structures between the two datasets (r = 0.99; p = 0.0001). To obtain the highest taxonomic resolution, the dataset based on the merged reads was used. Further analysis was done using the phyloseq v. 1.30.0 R package61.Statistical analysisAnalyses of variance (ANOVA) were conducted to examine the effects of N fertilized kernza at 100 kg N ha−1 (K100) and kernza at 200 kg N ha−1 (K200) as well as to test the effect of the deep-rooted plant species: kernza and lucerne, and soil depth on each of the dependent variables. An average across the two subplots within each of the three kernza field plots was used when measured variables did not significantly differ between K100 and K200. Subsequent pairwise comparisons of the means was conducted using the TukeyHSD post-hoc test. Homogeneity of variance and normality were confirmed (data log-transformed when required) for all comparisons using the Fligner-Killeen test of homogeneity of variances62 and the Shapiro–Wilk test of normality63. A permutational multivariate analysis of variance (PERMANOVA) using the Bray–Curtis dissimilarity matrix with the adonis function in the vegan R package was used to test the effect of K100 and K200, lucerne across both K100 and K200, and depth on the bacterial and fungal communities. The multivariate homogeneity of group dispersions or variances were confirmed for all comparisons using the function betadisper in vegan. The bacterial and fungal communities in response to the ascribed variables were visually represented as ordination plots with a Principle Coordinates Analysis (PCoA). Unique ASVs were defined for each depth and between K100, K200, and lucerne as ASVs only present in those samples belonging to a specific depth and treatment. Significance testing was conducted at p  More

Phylogenetic analysis with multiple markersThe final alignment of 5670 bp length contained 843 variable sites (14.7%). Missing data accounted for 53.5% of the alignment cells and the relative GC content was 39.5%. Our phylogeny suggests that the investigated Holarctic taxa of the niger clade sensu Ref.34 are divided into two major clades with strong statistical support (Fig. 1). The first major clade (L. niger group) consists exclusively of Palearctic species (L. niger, L. platythorax, L. japonicus, L. emarginatus, L. balearicus, L. grandis, L. cinereus, the L. alienus-complex, L. sakagamii, L. productus and L. hayashi), with the exception of an unnamed Nearctic subclade recovered as sister to the rest of the group. This unnamed subclade we describe as a new species below (L. ponderosae sp. nov.). Lasius ponderosae sp. nov. corresponds to what was previously known as the Nearctic form of “L. niger” sensu ref.17, but includes some western Nearctic populations formerly assigned to “L. alienus”17,52 as well. Monophyly of L. ponderosae sp. nov. was fully supported by Bayesian inference (pp = 1) and moderately supported by maximum likelihood (66% bootstrap support, Fig. 1). Lasius ponderosae sp. nov. is distantly related to L. niger; and L.niger is a close relative of L. japonicus and L. platythorax, as well as other Palearctic taxa. The second major clade (L. brunneus group) within the investigated Holarctic members of the L. niger clade contains both Nearctic and Palearctic species not closely related to the taxa of interest (Fig. 1).Figure 1Molecular phylogeny of 26 Holarctic ant taxa belonging to the subgenus Lasius sensu Wilson (1955) and two outgroup taxa (L. pallitarsis and L. mixtus). The phylogeny was calculated under the coalescent model and incorporates data from 9 genes (mtDNA: COI, COII, 16S, nuDNA: Defensin, H3, LR, Wg, Top1 & 28S). Names of species native to the Nearctic are shown in red and those of species native to the Palearctic in blue. Node labels show posterior probability (Bayesian inference) followed by bootstrap support (Maximum likelihood). The scale bar indicates the length of 0.01 substitutions/site.Full size imageDNA-barcodingThe native North American species L. ponderosae sp. nov. contains at least 15 COI-mitotypes (n = 28 sequenced specimens) belonging to four distinct deep lineages, with divergences of up to 5.9%. Haplotype diversity was 0.899 and nucleotide diversity was 0.012. None of the mitotypes of this species was found to be widespread or particularly abundant. In striking contrast, low genetic diversity was found in L. niger across its entire distribution (Fig. 2). No more than 7 different COI-mitotypes were detected in samples from distant localities representing most of the known range (n = 70 specimens from 12 countries), from Spain in the West to the Siberian Baikal-region in the East (Fig. 2). Their maximum pairwise divergence was only 0.6%, with a haplotype diversity of 0.682 and a nucleotide diversity below 0.001. One mitotype of L. niger is highly dominant within the native range, occurring from Western Europe to Central Siberia (mitotype h2 in Fig. 2).Figure 2Mitotype tree and distribution maps for 98 DNA-barcodes belonging to 7 mitotypes of the ant Lasius niger (blue, n = 70) and 15 mitotypes of L. ponderosae sp. nov. (red, n = 28). The red dashed line delimits the expected natural range of L. ponderosae sp. nov.53 Maps have been created using the free R-package “ggmap” v3.0.0 (https://github.com/dkahle/ggmap) in R v4.1.1. Map tiles by Stamen Design, under CC BY 3.0.Full size imageRecent Palearctic L. niger introduction to CanadaPalearctic Lasius niger was introduced to several localities in coastal Canada in recent times, where at least 11 populations were found in two metropolitan areas (Vancouver and Halifax areas, see Table S2 for details). Those populations consist of the most dominant Palearctic mitotype of L. niger (h2). However, in 3 localities in the Vancouver area, 3 specimens with a second mitotype were found (mitotype h4, Fig. 2, Table S2) in syntopy with those carrying the most common mitotype h2. This second Canadian COI-mitotype (h4) was not found among our samples from the Old World, although it only differs by a single nucleotide substitution from mitotypes found there. A review of BOLD data revealed that the Canadian barcoded specimens of L. niger were mostly collected in anthropogenic habitats such as schoolyards (Supplementary Table S2).Description of Lasius ponderosae sp. novLasius ponderosae Schär, Talavera, Rana, Espadaler, Cover, Shattuck and Vila. ZooBank LSID: urn:lsid:zoobank.org:act:22E2743A-2F1C-4870-B318-A1F2DF2B464C Etymology: ponderosae alludes to the ponderosa pine tree (Pinus ponderosa) that is at the centre of occurrence in the ponderosa pine—gambel oak communities in the western Rocky Mountains and northern Arizona.Type material: located at the Museum of Comparative Zoology, Cambridge, USA. Two paratype workers each will be deposited at the collections of University of California Davis (UCDC), the University of Utah (JTLC) and the Natural History Museum of Los Angeles County (LACM).Holotype: worker, Fig. 3a–c. Type locality: USA, Utah: Uintah Co., Uintah Mtns., 2408 m. 18.6 mi N. Jct. Rt. 40 on Rt. 191, 40.66378°N, − 109.47918°E, leg. 15.VII.2013, S. P. Cover; J. D. Rana, collection code SPC 8571. Measurements [mm]: HL: 0.899, HW: 0.823, SL: 0.821, EL: 0.239, EW: 0.189, ProW: 0.56, ML: 1.069, HTL: 0.863, CI: 92, SI: 100.Figure 3Frontal, lateral and dorsal view of the holotype worker (a–c), a paratype gyne (d–f) and a paratype male of Lasius ponderosae sp. nov. (g–i).Full size imageParatypes: 15 workers, two gynes (Fig. 3d–f), two males (Fig. 3g–i) from the same series as the holotype, morphometric data is given in the Appendix, Table S5 and Table S6. CO1 mitotype h17: Genbank Accession no. LT977508.Description of the worker caste: A member of a complex of cryptic species resembling L. niger. Intermediate in overall body size, antennal scape length and eye size and comparable to related species (Table 1). Terminal segment of maxillary palps and torulo-clypeal distance relative to head size shorter than in related Palearctic species (Table 1). Mandibles with 8 or rarely 7 or 9 regular denticles and lacking offset teeth at their basal angle. Penultimate and terminal basal mandibular teeth of subequal size, and the gap in between with subequal area than the basal tooth. Anterior margin of clypeus evenly rounded. Dorsofrontal profile of pronotum slightly angular (Fig. 4a). Propodeal dome short and flat, usually lower than mesonotum (Fig. 4a). Body with abundant and long pilosity, especially lateral propodeum, genae, hind margin and underside of head. Pilosity of tibiae and antennal scapes variable, ranging from almost no setae (“L. alienus”-like phenotype) to very hairy (“L. niger”-like phenotype). Microscopic pubescent hairs on forehead between frontal carinae long and fine. Clypeus typically with only few scattered pubescent hairs (Figs. 3, 4c). Coloration of body dark brown, occasionally yellowish- or reddish-brown or slightly bicolored with head and thorax lighter than abdomen. Femora and antennal scapes brown. Mandibles and distal parts of legs yellowish to dark brown. Specimens of all 3 castes are shown in Fig. 3a–i and morphometric data are summarized in Table 1 and raw measurements are available in Table S5 and S6.Table 1 Morphometric data of Lasius ponderosae sp. nov. and comparison to morphologically similar Palearctic species.Full size tableFigure 4Average thorax profile of Lasius ponderosae sp. nov. (a) and members of the Palearctic L. niger-complex (b). Figures were created by image averaging (L. ponderosae sp. nov n = 35; Palearctic L. niger-complex n = 30 specimens). Frontal view of head and detail of clypeus of the Holotype worker of L. ponderosae sp. nov. (c) and a non-type worker of L. niger (d).Full size imageDiagnosis: Lasius ponderosae sp. nov. workers key out to “L. niger” using Wilson’s 1955 key to the Nearctic Lasius species. However, some populations with reduced pilosity may also be identified as “L. alienus” using this key. Lasius alienus is a Eurasian species not known from North America33. The Nearctic “L. alienus” sensu Wilson (1955) includes both, L. americanus as well as populations of L. ponderosae sp. nov. with sparse setae counts on tibia and/or scapes. Lasius ponderosae sp. nov. can be distinguished from L. americanus by the presence of abundant, long setae surpassing the sides of the head in full face view (nGen  > 5 and nOcc  > 10 vs. nGen  0.8 across models and runs). The strongest predictors were: Annual Mean Temperature (mean variable importance = 0.32), Mean Temperature of Coldest Quarter (0.23), Temperature Annual Range (0.23) and Temperature Seasonality (0.24). The contribution of land cover was low (0.02). The model predicted high probabilities of occurrence of L. niger in the eastern United States and southeastern Canada, including the island of Newfoundland, and small areas of suitable habitat in southwestern Canada and the Aleutians (Fig. 6). The area with high predicted occurrence probability of L. niger in the New World includes the two sites where populations have actually established (which were not used in the modeling): Nova Scotia and Vancouver. Further areas with high occurrence probabilities are New England, Southern Ontario, the Great Lakes-region and the Northern Appalachians. Low occurrence probabilities were found for the central North American prairies as well as arctic, boreal, arid, subtropical and tropical regions (Fig. 6). Considering the highest occurrence probability range (0.8–1 on a 0–1 probability scale), the area of suitable habitats for L. niger is 4,547,537 km2 in Europe and 1,308,920 km2 in North America. For an intermediate to high occurrence probability range (0.5–1) we estimated 5,371,055 km2 in Europe and 3,054,283 km2 in North America, and for the widest probability range (0.2–1) we estimated 6,155,643 km2 of suitable areas in Europe and 6,889,745 km2 in North America (Fig. 6).Figure 6Projected occurrence probability from ecological niche modeling for the Palearctic ant Lasius niger which has been introduced to Canada, based on 19 climatic and one land use variable. The intensity of blue colour indicates the probability of occurrence on a 0–1 scale based on 180 presences (black circles) and 182 absences (white circles) in the native range in the Old World (a). The model was then projected to North America to estimate areas of suitable habitat for this introduced species (b). These maps have been created using the free R-package “ggplot2” v3.3.5 (https://ggplot2.tidyverse.org) in R v4.1.1.Full size image More

Paez-Espino D, Eloe-Fadrosh EA, Pavlopoulos GA, Thomas AD, Huntemann M, Mikhailova N, et al. Uncovering Earth’s virome. Nature. 2016;536:425–30.CAS 
PubMed 

Google Scholar 
Anderson PK, Cunningham AA, Patel NG, Morales FJ, Epstein PR, Daszak P. Emerging infectious diseases of plants: pathogen pollution, climate change and agrotechnology drivers. Trends Ecol Evol. 2004;19:535–44.PubMed 

Google Scholar 
Taylor LH, Latham SM, Woolhouse MEJ. Risk factors for human disease emergence. Philos Trans R Soc B Biol Sci. 2001;356:983–9.CAS 

Google Scholar 
White R, Murray S, Rohweder M. Pilot analysis of global ecosystems: grassland ecosystems. 2000 World Resources Institute. Washington, DC.Zhao Y, Liu Z, Wu J. Grassland ecosystem services: a systematic review of research advances and future directions. Landsc Ecol. 2020;35:793–814.
Google Scholar 
Trubl G, Jang HBin, Roux S, Emerson JB, Solonenko N, Vik DR, et al. Soil viruses are underexplored players in ecosystem carbon processing. mSystems. 2018;3:e00076–18.CAS 
PubMed 
PubMed Central 

Google Scholar 
Emerson JB, Roux S, Brum JR, Bolduc B, Woodcroft BJ, Jang HBin, et al. Host-linked soil viral ecology along a permafrost thaw gradient. Nat Microbiol. 2018;3:870–80.CAS 
PubMed 
PubMed Central 

Google Scholar 
Zablocki O, Adriaenssens EM, Frossard A, Seely M, Ramond J-B, Cowan D. Metaviromes of extracellular soil viruses along a Namib desert aridity gradient. Genome Announc. 2017;5:e01470–16.PubMed 
PubMed Central 

Google Scholar 
Jin M, Guo X, Zhang R, Qu W, Gao B, Zeng R. Diversities and potential biogeochemical impacts of mangrove soil viruses. Microbiome. 2019;7:58.PubMed 
PubMed Central 

Google Scholar 
Adriaenssens EM, Kramer R, Van Goethem MW, Makhalanyane TP, Hogg I, Cowan DA. Environmental drivers of viral community composition in Antarctic soils identified by viromics. Microbiome. 2017;5:83.PubMed 
PubMed Central 

Google Scholar 
Williamson KE, Fuhrmann JJ, Wommack KE, Radosevich M. Viruses in soil ecosystems: an unknown quantity within an unexplored territory. Annu Rev Virol. 2017;4:201–19.CAS 
PubMed 

Google Scholar 
Starr EP, Nuccio EE, Pett-Ridge J, Banfield JF, Firestone MK. Metatranscriptomic reconstruction reveals RNA viruses with the potential to shape carbon cycling in soil. Proc Natl Acad Sci. 2019;116:25900–8.CAS 
PubMed 
PubMed Central 

Google Scholar 
Wu R, Davison MR, Gao Y, Nicora CD, Mcdermott JE, Burnum-Johnson KE, et al. Moisture modulates soil reservoirs of active DNA and RNA viruses. Commun Biol. 2021;4:1–11.
Google Scholar 
Hurwitz BL, Sullivan MB. The Pacific Ocean Virome (POV): a marine viral metagenomic dataset and associated protein clusters for quantitative viral ecology. PLoS One. 2013;8:e57355.CAS 
PubMed 
PubMed Central 

Google Scholar 
Breitbart M, Bonnain C, Malki K, Sawaya NA. Phage puppet masters of the marine microbial realm. Nat Microbiol. 2018;3:754–66.CAS 
PubMed 

Google Scholar 
Wolf YI, Kazlauskas D, Iranzo J, Lucía-Sanz A, Kuhn JH, Krupovic M, et al. Origins and evolution of the Global RNA virome. MBio. 2018;9:e02329–18.PubMed 
PubMed Central 

Google Scholar 
Shi M, Lin XD, Tian JH, Chen LJ, Chen X, Li CX, et al. Redefining the invertebrate RNA virosphere. Nature. 2016;540:539–43.CAS 

Google Scholar 
Callanan J, Stockdale SR, Shkoporov A, Draper LA, Ross RP, Hill C. Expansion of known ssRNA phage genomes: from tens to over a thousand. Sci Adv. 2020;6:eaay5981.CAS 
PubMed 
PubMed Central 

Google Scholar 
Koonin EV, Dolja VV, Krupovic M, Varsani A, Wolf YI, Yutin N, et al. Global organization and proposed megataxonomy of the virus world. Microbiol Mol Biol Rev. 2020;84:e00061-19.PubMed 
PubMed Central 

Google Scholar 
Cobbin JC, Charon J, Harvey E, Holmes EC, Mahar JE. Current challenges to virus discovery by meta-transcriptomics. Curr Opin Virol. 2021;51:48–55.CAS 
PubMed 

Google Scholar 
Trubl G, Hyman P, Roux S, Abedon ST. Coming-of-age characterization of soil viruses: a user’s guide to virus isolation, detection within metagenomes, and viromics. Soil Syst. 2020;4:1–34. MDPI AG.
Google Scholar 
Santos-Medellin C, Zinke LA, ter Horst AM, Gelardi DL, Parikh SJ, Emerson JB. Viromes outperform total metagenomes in revealing the spatiotemporal patterns of agricultural soil viral communities. ISME J. 2021;15:1–15.
Google Scholar 
Adriaenssens EM, Farkas K, Harrison C, Jones DL, Allison HE, McCarthy AJ. Viromic analysis of wastewater input to a river catchment reveals a diverse assemblage of RNA viruses. mSystems. 2018;3:e00025–18.CAS 
PubMed 
PubMed Central 

Google Scholar 
Bibby K, Peccia J. Identification of viral pathogen diversity in sewage sludge by metagenome analysis. Environ Sci Technol. 2013;47:1945–51.CAS 
PubMed 
PubMed Central 

Google Scholar 
Culley A. New insight into the RNA aquatic virosphere via viromics. Virus Res. 2018;244:84–89.CAS 
PubMed 

Google Scholar 
Withers E, Hill PW, Chadwick DR, Jones DL. Use of untargeted metabolomics for assessing soil quality and microbial function. Soil Biol Biochem. 2020;143:107758.CAS 

Google Scholar 
Trubl G, Solonenko N, Chittick L, Solonenko SA, Rich VI, Sullivan MB. Optimization of viral resuspension methods for carbon-rich soils along a permafrost thaw gradient. PeerJ. 2016;4:e1999.PubMed 
PubMed Central 

Google Scholar 
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. J. 2011;17:10.
Google Scholar 
Joshi N, Fass J. Sickle: a sliding-window, adaptive, quality-based trimming tool for FastQ files. 2011.Schmieder R, Edwards R. Quality control and preprocessing of metagenomic datasets. Bioinformatics. 2011;27:863–4.CAS 
PubMed 
PubMed Central 

Google Scholar 
Kopylova E, Noé L, Touzet H. SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics. 2012;28:3211–7.CAS 
PubMed 

Google Scholar 
Li D, Liu C-M, Luo R, Sadakane K, Lam T-W. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics. 2015;31:1674–6.CAS 
PubMed 

Google Scholar 
Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2014;12:59–60. Nature Publishing Group.PubMed 

Google Scholar 
Huson DH, Beier S, Flade I, Górska A, El-Hadidi M, Mitra S. et al.MEGAN Community Edition – interactive exploration and analysis of large-scale microbiome sequencing data.PLOS Comput Biol. 2016;12:e1004957PubMed 
PubMed Central 

Google Scholar 
Hyatt D, Chen GL, LoCascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinform. 2010;11:119.
Google Scholar 
Mistry J, Finn RD, Eddy SR, Bateman A, Punta M. Challenges in homology search: HMMER3 and convergent evolution of coiled-coil regions. Nucleic Acids Res. 2013;41:e121–e121.CAS 
PubMed 
PubMed Central 

Google Scholar 
Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22:1658–9.CAS 
PubMed 
PubMed Central 

Google Scholar 
Roux S, Adriaenssens EM, Dutilh BE, Koonin EV, Kropinski AM, Krupovic M, et al. Minimum information about an uncultivated virus genome (MIUViG). Nat Biotechnol. 2018;37:29–37.PubMed 
PubMed Central 

Google Scholar 
Germain P-L, Vitriolo A, Adamo A, Laise P, Das V, Testa G. RNAontheBENCH: computational and empirical resources for benchmarking RNAseq quantification and differential expression methods. Nucleic Acids Res. 2016;44:5054–67.CAS 
PubMed 
PubMed Central 

Google Scholar 
Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. vegan: Community Ecology Package. 2019.Wickham H. ggplot2: elegant graphics for data analysis. 2016. Springer-Verlag New York.Conway JR, Lex A, Gehlenborg N. UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics. 2017;33:2938–40.CAS 
PubMed 
PubMed Central 

Google Scholar 
Katoh K. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30:3059–66.CAS 
PubMed 
PubMed Central 

Google Scholar 
Price MN, Dehal PS, Arkin AP. FastTree 2 – approximately maximum-likelihood trees for large alignments. PLoS One. 2010;5:e9490.PubMed 
PubMed Central 

Google Scholar 
Letunic I, Bork P. Interactive Tree Of Life (iTOL) v4: recent updates and new developments. Nucleic Acids Res. 2019;47:W256–W259.CAS 
PubMed 
PubMed Central 

Google Scholar 
Roux S, Emerson JB, Eloe-Fadrosh EA, Sullivan MB. Benchmarking viromics: an in silico evaluation of metagenome-enabled estimates of viral community composition and diversity. PeerJ. 2017;5:e3817.PubMed 
PubMed Central 

Google Scholar 
Ayllón MA, Turina M, Xie J, Nerva L, Marzano SYL, Donaire L, et al. ICTV virus taxonomy profile: botourmiaviridae. J Gen Virol. 2020;101:454–5.PubMed 
PubMed Central 

Google Scholar 
Krishnamurthy SR, Janowski AB, Zhao G, Barouch D, Wang D. Hyperexpansion of RNA bacteriophage diversity. PLOS Biol. 2016;14:e1002409.PubMed 
PubMed Central 

Google Scholar 
Hillman BI, Cai G. The family Narnaviridae. Simplest of RNA viruses. Adv Virus Res. 2013;86:149–76.
Google Scholar 
Obbard DJ, Shi M, Roberts KE, Longdon B, Dennis AB. A new lineage of segmented RNA viruses infecting animals. Virus Evol. 2020;6:61.
Google Scholar 
Xu X, Bei J, Xuan Y, Chen J, Chen D, Barker SC, et al. Full-length genome sequence of segmented RNA virus from ticks was obtained using small RNA sequencing data. BMC Genom. 2020;21:1–8.
Google Scholar 
Roossinck MJ. The good viruses: viral mutualistic symbioses. Nat Rev Microbiol. 2011;9:99–108. Nature Publishing Group.CAS 
PubMed 

Google Scholar 
Milgroom MG, Cortesi P. Biological control of chestnut blight with hypovirulence: a critical analysis. Annu Rev Phytopathol. 2004;42:311–38. Annual ReviewsCAS 
PubMed 

Google Scholar 
Zell R, Delwart E, Gorbalenya AE, Hovi T, King AMQ, Knowles NJ, et al. ICTV virus taxonomy profile: Picornaviridae. J Gen Virol. 2017;98:2421–2.CAS 
PubMed 
PubMed Central 

Google Scholar 
Valles SM, Chen Y, Firth AE, Guérin DMA, Hashimoto Y, Herrero S, et al. ICTV virus taxonomy profile: Dicistroviridae. J Gen Virol. 2017;98:355–6.CAS 
PubMed 
PubMed Central 

Google Scholar 
Barrios E. Soil biota, ecosystem services and land productivity. Ecol Econ. 2007;64:269–85.
Google Scholar 
Vainio EJ, Chiba S, Ghabrial SA, Maiss E, Roossinck M, Sabanadzovic S, et al. ICTV virus taxonomy profile: Partitiviridae. J Gen Virol. 2018;99:17–18.CAS 
PubMed 

Google Scholar 
Yong CY, Yeap SK, Omar AR, Tan WS. Advances in the study of nodavirus. PeerJ. 2017;2017:e3841.
Google Scholar 
Schmitt AP, Lamb RA. Escaping from the cell: assembly and budding of negative-strand RNA viruses. In: Kawaoka Y (ed). Biology of negative-strand RNA viruses: the power of reverse genetics. 2004. (Springer Berlin Heidelberg, Berlin, Heidelberg, pp 145–96.Käfer S, Paraskevopoulou S, Zirkel F, Wieseke N, Donath A, Petersen M, et al. Re-assessing the diversity of negative-strand RNA viruses in insects. PLoS Pathog. 2019;15:e1008224.PubMed 
PubMed Central 

Google Scholar 
Bejerman N, Debat H, Dietzgen, RG. The plant negative-sense RNA virosphere: virus discovery through new eyes. Front. Microbiol. 2020;11:588427.PubMed 
PubMed Central 

Google Scholar 
Wolf YI, Silas S, Wang Y, Wu S, Bocek M, Kazlauskas D, et al. Doubling of the known set of RNA viruses by metagenomic analysis of an aquatic virome. Nat Microbiol. 2020;5:1262–70.CAS 
PubMed 
PubMed Central 

Google Scholar 
Adriaenssens EM, Kramer R, van Goethem MW, Makhalanyane TP, Hogg I, Cowan DA. Environmental drivers of viral community composition in Antarctic soils identified by viromics. Microbiome. 2017;5:1–14.
Google Scholar 
Mahmoud H, Jose L. Phage and nucleocytoplasmic large viral sequences dominate coral viromes from the Arabian Gulf. Front Microbiol. 2017;8:2063.PubMed 
PubMed Central 

Google Scholar 
Koyama A, Steinweg JM, Haddix ML, Dukes JS, Wallenstein MD. Soil bacterial community responses to altered precipitation and temperature regimes in an old field grassland are mediated by plants. FEMS Microbiol Ecol. 2018;94:fix156.
Google Scholar 
Hurwitz BL, Hallam SJ, Sullivan MB. Metabolic reprogramming by viruses in the sunlit and dark ocean. Genome Biol. 2013;14:R123.PubMed 
PubMed Central 

Google Scholar  More