text:
Imaging of plant current pathways for non-invasive root Phenotyping using a newly developed electrical current source density approach

The flow of electric current in the root-soil system relates to the pathways of water and solutes, its characterization provides information on the root architecture and functioning. We developed a current source density approach with the goal of non-invasively image the current pathways in the root-soil system. A current flow is applied from the plant stem to the soil, the proposed geoelectrical approach images the resulting distribution and intensity of the electric current in the root-soil system. The numerical inversion procedure underlying the approach was tested in numerical simulations and laboratory experiments with artificial metallic roots. We validated the method using rhizotron laboratory experiments on maize and cotton plants. Results from numerical and laboratory tests showed that our inversion approach was capable of imaging root-like distributions of the current source. In maize and cotton, roots acted as “leaky conductors”, resulting in successful imaging of the root crowns and negligible contribution of distal roots to the current flow. In contrast, the electrical insulating behavior of the cotton stems in dry soil supports the hypothesis that suberin layers can affect the mobility of ions and water. The proposed approach with rhizotrons studies provides the first direct and concurrent characterization of the root-soil current pathways and their relationship with root functioning and architecture. This approach fills a major gap toward non-destructive imaging of roots in their natural soil environment.

Introduction
New root phenotyping technological developments are needed to overcome the limitations of traditional destructive root investigation methods, such as soil coring or "shovelomics" (Trachsel et al. 2011). Mancuso (2012) and Atkinson et al. (2019) provide extensive reviews on the methodological advances on non-destructive root phenotyping, including Bioelectrical Impedance Analysis (BIA), planar optodes, geophysical methods, and vibrating probe techniques. These techniques aim to mitigate key limitations of traditional root phenotyping, especially addressing the need for a better and more convenient characterization of the finer roots and root functioning. Advances in non-invasive and in-situ approaches for monitoring of root growth and function over time are needed to gain insight into the mechanisms underlying root development and response to environmental stressors.
Geophysical methods have been tested to nondestructively image roots in the field. Ground Penetrating Radar (GPR) approaches have been used to detect coarse roots (Guo et al. 2013;Liu et al. 2018). Electrical Resistivity Tomography (ERT) and ElectroMagnetic Induction (EMI) approaches have been used to image and monitor soil resistivity changes associated with the Root Water Uptake (RWU) (Cassiani et al. 2015;Shanahan et al. 2015;Whalley et al. 2017). Recent studies explored the use of multi-frequency Electrical Impedance Tomography (EIT) to take advantage of the root polarizable nature Kemna 2017, 2019).
Despite these advantages, geophysical methods to date share common limitations regarding root characterization. Geophysical methods developed to investigate geological media: in the case of roots they measure the root response as part of the soil response, see Fig. 1a for the ERT acquisition. Because of the natural soil heterogeneity and variability the resolution and signal characteristics of geophysical methods strongly depend on soil type and conditions. As such, interpretation of the rootsoil system response is non-unique, hindering the differentiation between roots of close plants and the extraction of specific information about root physiology from the electrical signals.
Unlike geophysical methods, the BIA for root investigation developed to specifically target the impedance of plant tissues, limiting the influence of the growing medium. A practical consequence is that BIA involves the application of sensors (i.e., electrodes) into the plant to enhance the method sensitivity. BIA measures the electrical impedance response of roots at a single frequency (resistance and capacitance methods) or over a range of frequencies (bioelectrical impedance spectroscopy). The measured BIA responses have been used to estimate root characteristics, such as root absorbing area and root mass (Dietrich et al. 2012;Mancuso 2012;Ozier-Lafontaine and Bajazet 2005). Estimation of these root traits is based on assumptions on the electrical properties of roots (see section Current Pathways in Roots). A key assumption is that current travels and distributes throughout the root system before exiting to the soil (Fig. 1b), with no leakage of current into the soil in the proximal root position (Fig.  1c). It is only in the former case, that the BIA signal would be sensitive to root physiology. Despite the physiological relevance of the BIA assumptions and the number of BIA studies, a suitable solution for the characterization of the current pathways in roots is missing. Thus far, only indirect information obtained from invasive and time-consuming experiments have been available to address this issue (e.g., excavation and trimming). Mary et al. (2018Mary et al. ( , 2019, and Mary et al. 2020) tested the combined use of ERT and Mise A La Masse (MALM) methods for imaging grapevine and citrus roots in the field. An approach hereafter called inversion of Current Source Density (iCSD) was used to invert the acquired data. The objective of this inversion approach is to image the density and position of current passing from the plant to the soil. The current source introduced via the stem distributes into "excited" roots that act as a distributed network of current sources (Fig. 1d). Consequently, a spatial numerical inversion of these distributed electric sources provides direct information about the root current pathways and the position of the roots involved in the uptake of water and solutes. The numerical approach used to invert for the current source density is a key component required for such an approach. Mary et al. (2018) used a nonlinear minimization algorithm for the inversion of the current source density. The algorithm consisted of gradient-based sequential quadratic programming iterative minimization of the objective functions (L1/L2 approach therein) described in Mary et al. (2018). The algorithm was implemented in MATLAB, R2016b, using the fmincon method. Because no information about the investigated roots was available, the authors based these inversion assumptions and the interpretation of their results on the available literature data on grapevine root architecture. Consequently, Mary et al. (2018) highlighted the need for further iCSD advances and more controlled studies on the actual relationships between current flow and root architecture.
In this study, we present the methodological formulation and evaluation of the iCSD method, and discuss , and iCSD (d) methods. The red cross and the blue dash are the current injection electrodes. The black dashes are the potential electrodes, used to image the electric potential field (ERT and iCSD). a) The ERT method characterizes the RWU through its impact on the soil moisture content and, consequently, on the soil resistivity. b) The BIA method investigates root traits based on the impedance response. In this case, since the current flows through the roots before entering the soil, the BIA is sensitive to the roots properties. c) In this BIA case, the current passes to the soil without first flowing through the roots. Under these conditions, the BIA is not sensitive to the root properties. d) The iCSD method images the positions where the current passes from the plant to the soil (small red crosses) by measuring the distribution of the electric potential its applications for in-situ characterization of current pathways in roots. We perform our studies using laboratory rhizotron experiments on crop roots. The main goals of this study were: 1) develop and test an iCSD inversion code that does not rely on prior assumptions on root architecture and function; 2) design and conduct rhizotron experiments that enable an optimal combination of root visualization and iCSD investigation of the current pathways in roots to provide direct insight on the root electrical behavior and validate the iCSD approach; and 3) perform experiments to evaluate the application of the iCSD method on different plant species (maize and cotton) and growing medium (soil and nutrient solution) that are common to BIA and other plant studies.

Current pathways in roots
The relationship between hydraulic and electrical pathways has been the object of scientific debate because of its physiological relevance and methodological implications for BIA methods (Urban et al. 2011;Dietrich et al. 2012).
A key and open question concerns the distribution of the current leakage (i.e., the root portions where the current passes from the root to the soil; Fig. 1b and c). The distribution of the current leakage is controlled by 1) the electrical radial and longitudinal conductivities (σ cr and σ cl ), and 2) by the resistivity contrast between root and soil.
With regard to σ cr and σ cl , when σ cl is significantly higher than σ cr , the current will predominantly travel through the xylems to the distal "active" roots, which are mostly root hairs. Based on the link between hydraulic and electrical pathways, this is consistent with a root water uptake (RWU) process where root hairs play a dominant role while the more insulated and suberized roots primarily function as conduits for both water and electric current (Aroca 2012, and references therein). On the contrary, if the σ cr is similar to σ cl , the electrical current does not tend to travel through the entire root system but rather starts leaking into the surrounding medium from root proximal portions. The coexistence of proximal and distal current leakage is in line with studies that suggest the presence of a more diffused zone of RWU, and a more complex and partial insulation effect of the suberization, possibly resulting from the contribution of the cell-to-cell pathways (e.g., Kramer 1946;Steudle 2000;North and Baker 2007;Ranathunge and Schreiber 2011;Schreiber 2010).
Soil resistivity can affect the distribution of the current leakage by influencing the minimum resistance pathways, i.e., whether roots or soil provide the minimum resistance to the current flow. In addition, soil resistivity strongly relates to the soil water content, which, as discussed, affects the root physiology. Therefore, information on the soil resistivity, such as the ERT resistivity imaging, has the potential for supporting the interpretation of both BIA and iCSD results. Dalton (1995) proposed a model for the interpretation of the plant root capacitance results in which the current equally distributes over the root system. Because of the elongated root geometry this model is coherent with the hypothesis of a low resistance xylem pathway (σ cl > σ cr ). Numerous studies have applied Dalton's model documenting the predicted correlation between root capacitance and mass (Dietrich et al. 2012, and references therein). In fact, recent studies with wheat, soy, and maize roots continue to support the capacitance method (Čermák et al. 2013;Preston et al. 2004;Postic and Doussan 2016;Cseresnyés et al. 2018).
Despite accumulating studies supporting the capacitance method, hydroponic laboratory results of Dietrich et al. (2012) and other studies have begun to uncover potential inconsistencies with Dalton's assumptions. In their work, Dietrich et al. (2012) explored the effect of trimming submerged roots on the BIA response and found negligible variation of the root capacitance. Cao et al. (2010) reached similar conclusions regarding the measured electrical root resistance (overall electrical resistance for the electric circuit obtained with injection into the stem and return electrode in the hydroponic solution, Fig. 1a). Urban et al. (2011) discussed the BIA hypotheses and found that the current left the roots in their proximal portion in several of their experiments. Conclusions from the latter study are consistent with the assumption that distal roots have a negligible contribution on root capacitance and resistance.
Because of the complexity of the hydraulic and electrical pathways, their link has long been the object of scientific research and debate. For recent reviews see Aroca (2012) and Mancuso (2012); for previous detailed discussions on pathways in plant cells and tissues see Fensom (1965), Knipfer and Fricke (2010), and Findlay and Hope (1976); see Johnson (2013) and Maherali et al. (2009) in regard to xylem pathways. See Jackson et al. (2000) and Hacke and Sperry (2001) for water pathways in roots. Thus, above discrepancies in the link between electrical and hydraulic root properties can be, at least to some degree, attributed to differences among plant species investigated and growing conditions. Among herbaceous plants, maize has been commonly used to investigate root electrical properties (Frensch and Steudle 1989). For instance, Ginsburg (1972) investigated the longitudinal and radial current conductivities (σ cl and σ cr ) of excited root segments and concluded that the maize roots behave as leaking conductors. Similarly, Anderson and Higinbotham (1976) found that σ cr of maize cortical sleeves was comparable to the stele σ cl . Recently, Rao et al. (2018) found that maize root conductivity decreases as the root cross-sectional area increases, and that primary roots were more conductive than brace roots. By contrast, BIA studies on woody plants have supported the hypothesis of a radial isolation effect of bark and/or suberized tissues (Čermák et al. 2006;Aubrecht et al. 2006). Plant growing conditions have been shown to affect both water uptake and solute absorption due to induced differences in root maturation and suberization (e.g., Schreiber 2010). Redjala et al. (2011) observed that the cadmium uptake of maize roots grown in hydroponic conditions was higher than in those grown aeroponically. Tavakkoli et al. (2010) demonstrated that the salt tolerance of barley grown in hydroponic conditions differed from that of soil-grown barley. Zimmermann and Steudle (1998) documented how the development of Casparian bands significantly reduced the water flow in maize roots grown in mist conditions compared to those grown hydroponically. During their investigation on the effect of hypoxia on maize, Enstone and Peterson (2005) reported differences in oxygen flow between plants grown hydroponically and plants grown in vermiculite.
The results reported above and in other investigations (Aroca 2012, and references therein) are conducive to the hypothesis that root current pathways are affected by the growing conditions, as suggested in Urban et al. ( 20 11 ) . F o r ex a m p l e, t h e ob s e r v a t i on s b y Zimmermann and Steudle (1998) and Enstone and Peterson (2005) may explain the negligible contributions to the BIA signals from distal roots under hydroponic conditions (Cao et al. 2010;Dietrich et al. 2012). At the same time, the more extensive suberization in natural soil and weather conditions could explain the good agreement between the rooting depth reported by Mary et al. (2018) based on the iCSD and the available literature data for grapevines in the field.
To minimize these ambiguities and to develop a more robust approach for non-invasive in-situ root imaging, we aim to develop iCSD inversion code that does not rely on prior assumptions on root architecture and function and use rhizotron experiments to validate the iCSD approach.

Materials and methods
The phrase "inversion of Current Source Density" (iCSD) was introduced by Łęski et al. (2011) to describe the 2D imaging of current sources associated with the brain neural activation. Similar inversion methodologies have been developed for the interpretation of the selfpotential data, where the distribution of naturally occurring currents is investigated (Minsley et al. 2007, and references therein). With regard to active methodologies, Binley et al. (1997) developed an analogous approach for detecting pollutant leakage from environmental confinement barriers. Although there are physical and numerical intrinsic differences between application of the iCSD to detect brain neuronal activity and current pathways in roots, we decided to adopt the term iCSD as the general physical imaging of current source density remains valid. With iCSD, we indicate the coupling of ERT and MALM through the proposed numerical inversion procedure for the imaging of the current source density, and its correlation with root architecture.
We introduce the necessary aspects regarding the ERT and MALM methods in this section. However, we direct the interested readers to more in-depth discussion about the ERT method (Binley and Kemna 2005), and to Schlumberger (1920) and Parasnis (1967) with regard to the MALM method.
In the following discussion we use ρ med to represent the 2D or 3D distribution of the electrical resistivity in the growing medium (i.e., hydroponic solution or soil, and with possible influences from RWU and ligneous root mass). CSD represents the 2D, or 3D, distribution of the Current Source Density within the same medium. In the case of roots, the CSD is controlled by the current conduction behavior of the roots, specifically by the leakage pattern of the root system (e.g., proximal or distal current leakage, Fig. 1).
Both ERT and MALM are active methods. In these methods the current is forced through the medium by applying a potential difference between two current electrodes. In ERT, both current electrodes are positioned in the investigated medium, while for MALM the positive current pole is installed in the plant stem, similar to BIA (Fig. 1). The potential field resulting from the current injection depends on CSD, resistivity of the medium (ρ med ), and boundary conditions. The boundary conditions are known a priori and their impact on the potential field can be properly modeled. In ERT, the current sources correspond to the electrodes used to inject current, allowing us to invert for ρ med . Then, the iCSD accounts for the obtained ρ med and explicitly inverts the MALM data to obtain current source distribution.

Laboratory setup and data acquisition
The rhizotrons used in this study were designed to enable the concurrent direct visualization of the roots and electrical measurements. Rhizotron dimensions were 52 cm (x) × 53 cm (y) × 2 cm (z), see Fig. 2. Figure 2a shows the rhizotron setup with 64 silver/silver chloride (Ag/AgCl) electrodes located on the back viewing surface. The viewing surfaces were covered with opaque material to stop the light from affecting the development of the roots. The back viewing surface was removable, allowing homogeneous soil packing for the plant experiments and convenient access to the electrodes. Besides the top opening, the rhizotrons were waterproof to enable hydroponic experiments and controlled evapotranspiration conditions during the soil experiments and plant growth. All the experiments were performed in a growth chamber equipped with automatic growth lights (LumiGrow Pro650e) and controlled temperature and humidity. The temperature varied with a day/night temperature regime of 25/20°C. The humidity ranged from 45 to 60%.
For both ERT and MALM methods, the electrical potential field is characterized by a set of potential differences (ΔV) measured between pairs of electrodes. It is important to properly arrange the electrodes on the rhizotron viewing surface and design a suitable acquisition sequence to obtain a good sensitivity coverage (hereinafter coverage) of the investigated system (Slater et al. 2002). This is particularly true for the iCSD, as both ERT and MALM acquisitions affect its result.
The 64 electrodes were arranged in a 8 by 8 grid on the back viewing surface of the rhizotron, leaving the front surface clear for the observation (Fig. 2). For the ERT, the designed arrangement of the electrodes offers a good compromise between a high coverage on the central part of the rhizotron, which encompasses the root zone, and a sufficient coverage on the rhizotron sides to avoid an excessive ERT inversion smoothness.
For the MALM, the arrangement of the electrodes is highly sensitive to the position of the investigated current sources. Because of their central positions, the electrodes are closer to the expected sources of current (i.e., the roots) and thus in the region of maximum potential gradient. Hence, this electrode configuration maximizes the changes in both magnitude and sign of the measured ΔV associated with a change in the CSD distribution.
The electrode diameter was 1.5 mm. The penetration of the electrodes into the rhizotron was 4 mm ± 1 mm. To evaluate the possible distorting effects of the densely populated electrodes on the potential field distribution, a test was performed with low conductivity water (20 μS cm −1 ). The test showed no resistivity anomalies, which may be caused by the presence of the electrodes (data not shown). Therefore, while rhizotron setups with electrodes only on the sides were successfully adopted (Weigand and Kemna 2017), we found that the current setup represents a better solution for iCSD experiments (Fig. 2).
Data were acquired with a MTP DAS-1 resistivity meter with 8 potential channels. For the ERT acquisition over the 2D grid of electrodes, we chose a dipole-dipole skip 2 configuration. For each skip 2-couple of injection electrodes (e.g., electrodes 1 and 4) the remaining skip-2 couples of electrodes were used as potential dipoles (e.g., electrodes 5 and 8, 6 and 9, etc.). The associated complete set of reciprocals was also acquired, the resulting acquisition sequence contained 3904 data points (Binley and Kemna 2005;Parasnis 1988;Mary et al. 2018).
Following the ERT data acquisition, the MALM data acquisition required little setup adjustments and time. As the two current electrodes are fixed, the use of a multichannel resistivity meter significantly reduced the acquisition time and, consequently, supported the acquisition of more robust data sets. Electrode 1 was used to inject the current into the plant stem, while electrode 64 was used as a return electrode in the growing medium (Fig. 1d). The remaining 62 electrodes were used to map the resulting potential field. A sequence with 204 ΔVs was used. Considering the grid in Fig. 2a, the sequence included the vertical, horizontal, and diagonal ΔVs between adjacent electrodes. While 61 ΔVs would provide all the independent differences, the 204 ΔV sequence was preferred because of its redundancy and consequent lower sensitivity to acquisition errors. The acquisition time remained relatively short (2-3 min) as the multichannel instrument was optimized with fixed current electrodes that allowed 8 ΔVs to be measured at once.

Data processing and ERT inversion
ERT and MALM data were filtered considering: 1) reciprocal error (5% threshold), 2) stacking error (3 stacks, 5% error threshold), 3) minimum measured ΔV (minimum 0.5 mV threshold), 4) apparent resistivity (the range changed among different data sets), and 5) maximum electrode contact resistance (30 kohm). The filtering was implemented to enhance the control on the ERT inversion and obtain a reliable ρ med for the successive iCSD inversion.
After the processing, the data were inverted with the BERT inversion software to obtain ρ med (Günther and Rücker 2013;Rücker et al. 2006;Günther et al. 2006). Data error was set to 5% in line with the stacking and reciprocal thresholds used for the data filtering. The regularization was adjusted using the lambda optimization algorithm provided by the BERT software. Generally speaking, a rhizotron is treated as a 2D geometry. However, this bounded and thin geometry leads to some complications into the ERT inversion. Specifically, nocurrent-flow boundary conditions were set for all the rhizotron surface and the inversion had to be adapted for the resulting pure Neumann problem (Bochev and Lehoucq 2005). For a higher quality forward calculation, the rhizotron volume was discretized in 3D by extruding a 2D mesh with 5 layers. The discretization allowed the refinement of the mesh near the electrodes while maintaining high mesh quality. At the same time, in order to limit the inversion time and force the inversion to be two-dimensional (x, y), the elements in z direction (resulting from the extruded layers) were grouped together and inverted as a unique variable. This way, a 3D forward calculation was implemented within a 2D inversion (Ronczka 2016).

iCSD inversion
The iCSD inversion that we developed was based on the physical principles of a bounded system in which linearity and charge conservation were applied to decompose the investigated CSD distribution into the sum of point current sources. This provided a discrete representation of the root system portions where the current leaks from the roots into the surrounding medium.
Because of the linearity of the problem, the collective potential field from multiple current sources is the linear combination of their individual potential fields. As such, the measured ΔV can be viewed as and decomposed into the sum of multiple ΔVs from a set of possible current sources. These possible current sources are named ViRTual electrodes (VRTe). As purely numerical (virtual) electrodes, they are simulated by mesh nodes representing possible current sources, but with no direct correlation with the real (physical) electrodes used during data acquisition. Basically, the VRTe were distributed to represent a grid over which the true CSD distribution is discretized. In order to account for any possible CSD, a 2D grid of 306 VRTe was arranged to cover the entire rhizotron (Fig. 2b).
The charge conservation law implies that the sum of the current fractions associated with the VRTe (hereinafter VRTe weights) has to be equal to the overall injected current, which is provided by the resistivity meter. If we normalize the injected current to be equal to 1, the sum of the VRTe weights has to be 1 as well. Briefly, for Ohm's law, normalizing the current to 1 is equivalent to calculating the resistance, R, from ΔV. Then, the use of R simplifies the presentation of the numerical problem.
Once the VRTe nodes are added to the ERT-based ρ med structure, the potential field associated with each of the VRTes is simulated with BERT. From these simulated potential fields, the same sequence of 204 R is extracted, each corresponding to a single VRTe. Each extracted sequence contains the resistances that would be measured in the laboratory if all the current sources were concentrated at the VRTe point (i.e., VRTe Green's function).
The VRTe weights are the unknowns that the iCSD strives to estimate. Once the VRTe weights are estimated and associated with the respective VRTe coordinates, they provide a discrete visualization of the investigated CSD. The problem of estimating the VRTe weights that decompose the measured R sequence into a sum of the simulated VRTe R sequences is analogous to a linear vector decomposition (Strang 1976), and is expressed by: Where A is a matrix, its columns are the simulated VRTe R sequences; x is a vector containing the unknown VRTe weights; b is a vector containing the measured sequence of 204 R. Each row in A corresponds to the relative R in the acquisition sequence, e.g., A 1,1 is the first resistance extracted from the potential field simulated with injection at the first VRTe.
The charge conservation is implemented by appending a row of 1's to A and a corresponding 1 to the vector b. This forces the sum of the VRTe weights to be equal to 1. A weight of 1000 was used to guarantee the charge conservation; with a lower weight observations and regularization may dominate the charge conservation.
Since the problem is undetermined, a first order spatial regularization is added (Menke 1989). Rows are added to express the differences between adjacent (vertically and horizontally) VRTe, e.g., the row [1 − 1 …] is the difference between the first two VRTe weights. The differences are added for the entire VRTe grid and set to 0 by adding corresponding 0's to b.
Lastly, the trade-off between data misfit and solution regularization is controlled by a diagonal weight matrix W. The numerical routine includes a "pareto" functionality wherein regularization and model-to-measurement fit are traded off while changing the regularization weight, i.e., running the inversion with different regularization weights. The obtained set of solutions can be used to construct the "pareto front" (L-curve), which is a widely accepted way to estimate the optimum regularization weight (Hansen and Dianne 1993). Eq. 2 shows the resulting system WAx = bW. Where S indicates the n VRTe sources, R the m resistances, and w the weights used to control the solution regularization.
Lastly, the solution is further constrained by forcing the linear solver to seek only positive VRTe weights (i.e., inequality constraint), as the negative source of current is known to correspond uniquely to the return electrode.
The linear problem formulation is conducive to the inversion optimization during the calculation of the pareto front. The calculation time of the Pareto front can be further reduced by code optimization as the calculations that do not depend on the regulation weights can separated from the inversion routine and performed only once during the initialization of the linear problem. The initialization phase includes the processing of the MALM experimental data, forward calculation of the VRTe responses for the given ρ med , inclusion of the continuity constraint, and construction of the matrices.
Continuity constraint, bounded-value constraint, and first-order spatial regularization stabilize the inversion while limiting the impact of the spatial regularization strategy. The impact of the spatial regularization was evaluated by monitoring the relative components of the misfit and the resulting distribution of the current source. In both synthetic and laboratory tests, as well as in plant experiments the iCSD results are often limited to few current sources (Figs. 3, 4, 7, and 9).

Synthetic and experimental iCSD testing
Synthetic numerical and laboratory experimental tests were performed in order to evaluate the capabilities of the setup and inversion routine to couple the ERT and MALM approaches for the iCSD. In the numerical tests both the true source response and VRTe responses were calculated with BERT. Figure 3 shows an explanatory numerical test with inversion of a point source, and the associated Pareto front that was used to select the optimum regularization strength. As this first experiment was performed to specifically test the inversion routine, a homogeneous ρ med was used in order to avoid influence from the baseline resistivity distribution complexity.
For the second experiment, the laboratory tests were conducted. Because of the ρ med heterogeneity of any experimental system, these laboratory tests need to include the ERT inversion, and the use of the obtained ρ med as input in the iCSD. The true current sources were obtained using insulated metallic wires inserted into the rhizotron (Fig. 4). The insulating plastic cover was removed at the tips (1 cm, red dots in Fig. 4a) of the metallic wires to obtain the desired current sources. Six experimental tests were performed using different numbers and positions of these current sources. The rhizotron was filled with tap water and left to equilibrate to achieve steady state conditions of water temperature and salinity, thus minimizing ρ med heterogeneity and changes during the experiment. Changes in ρ med during the ERT and MALM acquisition periods would make the ERT-based ρ med less accurate and compromise the iCSD. To make sure ρ med was stable, a second ERT was performed after the MALM acquisition and compared with the initial measurement. The conductivity of the solution was also measured in several locations of the rhizotron with a conductivity meter to validate the Fig. 3 Numerical synthetic test. a) The synthetic test evaluates the iCSD capability to locate a current point source (white dot). The white circles represent the VRTe and the white square the return electrode. The source of current was centered between four VRTe to highlight the distribution of the current density. While the fractions of current source refer to the VRTe nodes, the results are linearly interpolated for a better visualization. b) Pareto front used to select the proper regularization weight for the iCSD, the blue-circled inversion was chosen as best compromise between fitting and regularization ρ med obtained from the ERT inversion. This setup allowed the acquisition of good quality data sets since less than 5% of the data were discharged during the data processing. Because of the controlled laboratory conditions, the ρ med obtained with the ERT was stable and consistent with the direct conductivity measurements. The quality of the ERT inversion was also confirmed by comparing the model responses with the acquired data (i.e. data misfit). Similarly, the acquired iCSD data were plotted against the resistances calculated with the CSD distribution obtained from the iCSD.
The tests also allowed a more informed definition of the VRTe grid. For our setup, a spacing of 3 cm provided a good compromise between resolution, stability, and duration of the iCSD routine. The 3-cm spacing also agrees with the ERT resolution, which would not support a higher iCSD resolution.
Successive numerical tests were based on the 8source laboratory tests shown in Fig. 4. These tests aimed to 1) link laboratory and numerical tests to evaluate the influence of the numerical iCSD routine and laboratory setup on the overall iCSD stability and resolution; 2) account for a more complex CSD, given by the 8 wire-tip sources that were used to simulate distal current pathways; and 3) account for possible ρ med heterogeneity. To address goals 1 and 2, the position of the 8 sources was replicated in the numerical tests and a test with homogeneous ρ med was included to simulate the water resistivity of the laboratory tests. To address goal 3, heterogeneous ρ med were tested. Fig. 4 Laboratory tests. a) Rhizotron filled with water for the tests with metallic wires. The plastic insulation of the metallic wires avoided the current leakage, which occurred only where the insulation was removed (eight red dots). This test simulated a complex and root-like structure, with the current leaking at the tips. b) result of iCSD test with two metallic sources in water. c) result of iCSD test with eight metallic sources shown in a). In both tests, the unequal distribution of the current among the wires is due to different proximity of the respective sources to the return electrode: the shorter the water path from the source to return electrode, the lower the path resistance, and the bigger the current fraction In order to account for the heterogeneous ρ med the following modeling steps were carried out. First, a true ρ med was assigned to the mesh cells of the rhizotron ERT model. We included homogeneous, linear, and quadratic resistivity profiles in the y direction, see Fig. 5. Second, the ERT acquisition was simulated (i.e., forward calculation) with the ERT laboratory sequence and 3% of Gaussian error, in line with reciprocal and stacking errors observed in the laboratory data sets. Third, the forwarded ERT data sets were inverted following the exact laboratory procedure. A refined and different mesh was used for forward and inverse problems to, respectively, increase the simulation accuracy and avoid the inverse crime (Wirgin 2004). The ERT forward calculation was then repeated over the inverted ρ med. The obtained inverted responses were compared with the responses of the true models.
As for ERT, we compared true and inverted MALM responses. First, the true response was simulated with the 8 current sources overt the true (not inverted) ρ med . Second, a MALM response was calculated over the inverted ρ med and inverted to obtained the inverted CSD. Third, the obtained inverted CSD was used to forward calculate the inverted MALM response over the inverted ρ med . True and inverted MALM responses were then compared.

Plant experiments
We performed hydroponic and soil experiments using maize and cotton plants. In all the plant experiments, the injection electrode was positioned in the plant stem at a height of 1 cm from the surface of the growth media.
For the hydroponic experiments, the plants were first grown in columns with aerated nutrient solution (Hoagland and Arnon 1950). They were then moved to the rhizotron for the experiments. As in the metallic roots test, the rhizotron was filled 1 day before the experiment to reach stable and homogeneous temperature and salinity conditions. The plant was positioned at the center of the rhizotron with soft rubber supports. The plants were submerged at the same level as in the growing column to avoid discrepancies caused by the plant tissue adaptation to the submerged and aerated conditions, as discussed above with regard to the growing conditions. Consequently, the root crown was approximately 3 cm below the water surface.
For the soil experiments, seedlings were grown directly in the rhizotron to avoid damaging the roots and altering the root-soil interface. The soil was prepared by mixing equal volumes of sandy and clay natural soils acquired from an agricultural study site run by U.C. Davis, CA (Russell Ranch). The plants were irrigated with double strength Hoagland solution (Hoagland and Arnon 1950).
Two soil experiments were performed. In the first experiment, four cotton plants were grown for four months. For these experiments, the plants were positioned with the root crown approximately 8 cm deep (y = 0.44 m). In the second experiment, a pregerminated maize seed was planted 3 cm deep and then grown for four months (Fig. 10).
Results iCSD testing Figure 3 shows the result of a synthetic numerical test performed to evaluate the iCSD resolution, inversion stability, and influence of imposed constraints. The obtained CSD matches the true position of the simulated current source. The sum of the current sources equals 1 as expected and required by the continuity constraint. The resolution of the CSD is in line with the electrode interspace (5 cm). The first-order regularization (smoothness constraint) does not hinder the reconstruction of simulated point source. Figure 4 shows the results of a laboratory experiment where the iCSD method was tested with a known distribution of current sources obtained with metallic wires. The use of metallic wires offered a comprehensive solution to test the overall correct functioning of the laboratory setup and inversion routine. The iCSD correctly characterized both position and intensity of the test sources with no need for prior information to constrain the solution. The asymmetrical distribution CSD agrees with the principle of path of least resistance, i.e., the closer to the return electrode, the shorter the water path, the higher the current source density. The obtained distribution of current source agrees with true CSD and resolution expected on the basis of previous numerical tests. The iCSD resolution was sufficient for the imaging of a root-like structure, e.g., differentiate between distal, proximal, or equally distributed current pathways in the roots (Fig. 4c). In this case, a distal CSD was investigated (tips of the metallic wires) and the iCSD correctly found no proximal current source (near the immersion point of the root-like metallic structure). Figures 5, 6, and 7 show the results of the synthetic tests based on the 8-source laboratory tests. Figure 5 compares the true (a) and inverted (b) ρ med to highlight the overall influence of rhizotron meshing, acquisition sequence, processing, and inversion. Figure 6a offers a more quantitative comparison between true and inverted ERT responses. Similarly, Fig. 6b presents the correlation between true iCSD response (true ρ med and true CSD) and the inverted iCSD response (inverted ρ med and inverted CSD). ERT and iCSD inverted responses correlate with the associated true responses (all correlation coefficients >0.95). Figure 7 shows the inverted CSD distribution and confirms the iCSD capability to characterize both intensity and position of the current sources with a resolution of 3 cm.
Plant Experimental Results.
The plant experiments investigated the root current pathways combining the iCSD with the direct observation of the roots. Hydroponic and soil experiments were performed with maize and cotton plants. Figure 8a shows the root architecture for one of the two plant hydroponic experiments with maize plants. Although the roots were visually observed to extend deeper in the rhizotron, both the inverted CSD distributions show that the current leakage in this hydroponic system primarily occurred in the proximal root part at the top of the rhizotron. In particular, the highest current density was predominantly concentrated near the root crown with the remaining current leakage already occurring along the submerged stem section. The ρ med obtained from the ERT inversion confirmed the expected homogeneous temperature and salinity conditions.
In the soil experiments with cotton plants, the iCSD method successfully located the root position of the four plants along the soil surface and produced consistent signals corresponding to the root crowns (Fig. 9). The results also show a localized current leakage near the root crowns with no appreciable current leakage along the four stem sections in the soil or in the distal portions of the root systems. By contrast to the cotton results, the CSD distribution obtained from the maize experiment conducted within the soil presents current leakage along the stem section (Fig. 10).
The ERT-based ρ med obtained from the cotton soil experiments presents a strong heterogeneity with a highly resistive layer at the top of the rhizotron, representing a dry soil layer (Fig. 9b). The evolution of this dry layer is due to water loss through evapotranspiration. The ERT results highlight the need to account for the heterogeneity of ρ med in order to correctly calculate the responses of the VRTes during the iCSD. The ERT results also show how the ERT characterization can support the interpretation of the iCSD results by providing information on the soil characteristics and thus on possibly associated physiological responses, such as drought and suberization.

iCSD routine results
Both numerical and laboratory iCSD tests successfully image the CSD without requiring prior assumptions on the actual distribution of the roots and CSD. These assumptions were present in previous works (e.g., Mary et al. 2018) but were here replaced by application-specific constraints (smoothness, charge conservation, and positivity). In this sense, not only the explicit linear formulation allowed a direct implementation and evaluation of the these constraints, but also enabled full characterization of the resulting liner problem (e.g., condition number, uncertainty analysis; see Strang (1976) and John (2015)).
Relative to the field setup of Mary et al. (2018), the laboratory setup offered the advantages of more flexible positioning of the electrodes and consequently better ERT and MALM data coverage. The other significant advantages over previous studies were the reduced number of VRTe in a 2D distribution and the new iCSD routine resulting from the 2D distribution of the VRTe and trans-dimensional ERT inversion.
The inversion times was significantly reduced by the proposed iCSD, particularly thanks to its linear formulation, use of specific constraints and code optimization. Solving the linear system described above (Eq. 2) took less than 3 s on a standard laptop, compared to a few minutes using the generic optimization process used in Mary et al. (2018). The reduced inversion time also offers the advantage of a faster calculation of the Pareto front, which allowed a more informed inversion regularization.
Finally, the coupling between ERT and iCSD through the python geophysical library eased the optimization of the entire routine (Rücker et al. 2017). Coupling ERT Fig. 6 Correlations between ERT true and inverted responses (a), and MALM true and inverted responses (b). A homogeneous resistivity of 10 Ω m was used for the homogeneous tests; a linear increment of the resistivity from 10 to 100 Ω m (bottom to top) was used for the linear tests; a quadratic increment of the resistivity from 10 to 750 Ω m was used for the quadratic test (Fig. 5). All linear correlation coefficients >0.95 Fig. 7 CSD obtained from the 8source numerical test with quadratic resistivity vertical profile (Fig. 5). The test is based on the laboratory experiment with metallic roots, see Fig. 4 and iCSD is a relevant aspect of our procedure because the ERT data processing/analysis routine has a strong influence on the successive iCSD inversion, which may require the entire procedure to be repeated to test different ERT inversion parameters (e.g., Fig. 9).

Discussion
The iCSD results showed that the current leakage occurred in the very proximal regions of the root systems in both soil and hydroponic conditions. The proximal leakage was observed despite the return electrode being placed at the bottom of the rhizotron to allow deep current pathways. Nonetheless, the expected influence of the return electrode position was observed in the laboratory test with metallic roots and motivates the use of this electrode configuration in future laboratory experiments.
Our results are consistent with the early studies on maize root electrical properties (Ginsburg 1972), and corroborate the recent works that questioned the assumptions of the BIA methods (Dietrich et al. 2012;Urban et al. 2011). The high resistivity of the top layer of soil (Fig. 9) is expected to induce root suberization (see Introduction): our results would support the physiological hypothesis of water and nutrient absorption through older and possibly suberized roots (Aroca 2012;Kramer 1946).
It is worth noticing that the absence of current leakage along the section of the cotton stems in the top dry soil supports the assumption that the electrical structure of roots controls their current conduction behavior. Suberized epidermal cells can affect the movement of ions and, consequently, the current conduction in roots (Cseresnyés et al. 2013;Čermák et al. 2006). Therefore, it is feasible for current to be conducted along deeper portions of the more woody roots with minor leakage. The BIA experimental results that have observed positive correlations between electrical signals and root area are likely a result of physiological correlations between the root regions that contribute to the current flow and hair roots, which contributes most to functional root surface area (Fig. 1c). While correlations between BIA electric signals and investigated root traits appear to be indirect, the correlations observed across experimental platforms and species continue to validate its value for in-situ root phenotyping (Dietrich et al. 2012).
In their study on field grapevines, Mary et al. (2018Mary et al. ( , 2019Mary et al. ( , 2020 concluded that the CSD could be used to infer the root depth distribution. Because of the significant suberization of major roots in grapevines and orange trees, the electric current could penetrate deep into the root system before significant leakage occured. The deeper current penetration allowed the iCSD method to access and phenotype the root system. On the contrary, limited current results in lower sensitivity of the iCSD to distal and younger parts of the root system that are likely dominated by less suberized, finer roots. We attribute Fig. 8 Hydroponic experiments with a maize plant. a) Position of the plant and roots, with the root crown approximately 3 cm below the water surface. b-c) CSD obtained from the iCSD: the current leakage in proximal root part is dominant, while no appreciable current fraction reaches the distal part of the root system differences in current penetration among root systems of maize, cotton, grapevine, and orange tree to the differences in physiological traits such as the extent of suberization and lignification. If on one hand the sensitivity to the root physiological traits is a promising opportunity, on the other hand it has to be accounted for when phenotyping more herbaceous roots.

Conclusions
In this study we present a novel numerical inversion routine for inversion of current source density in a bounded and non-homogeneous medium. The new inversion procedure offers an explicit matrix formulation that takes advantage of the problem Fig. 9 Experiment with four cotton plants in soil. a) Rhizotron at the moment of the experiment with the root architectures and the dry soil layer at the top. The root crowns were approximately 8 cm deep (white dots). The gray line highlights the extension of the roots. b) ρ med obtained from the ERT inversion highlighting the strong variability of the soil resistivity that has to be taken into account for a proper iCSD. ρ med also confirms the presence of the top dry layer visible in a), which is caused by the evapotranspiration. c-g) iCSD results, in all the four experiments the current leakage concentrates near the root crown with no leakage along the stem above linearity and constraints. The improvements provided by our numerical inversion technique enable a faster and more controlled imaging of the current source density. The numerical inversion routine was also coupled with an existing ERT inversion code to ease the overall inversion procedure. The iCSD inversion code was based on top of an open source geophysical Python library (Rücker et al. 2017).
After testing the method using synthetic numerical and laboratory experiments, we designed a specific rhizotron setup for laboratory plant experiments with maize and cotton. We found that with the obtained electrical coverage, the iCSD correctly characterized the different CSD distributions used in the tests, including the imaging of a root-like structure (Figs. 4, 6, and 7). By contrast to the results reported by Mary et al. (2018), the experimental setup and the new numerical inversion procedure did not require prior assumptions on the root distribution to stabilize and regularize the inversion. By eliminating the need for imputing prior information, our inversion protocol excludes the possibility of obtaining biased solutions.
In all seven plant experiments, the current left the root pathway in its proximal portion, radially leaking into the surrounding growth medium. These results corroborate recent studies that questioned the BIA assumptions (Urban et al. 2011;Dietrich et al. 2012). While these studies provided a thoughtful discussion on the topic, they were based on time-consuming, invasive and less direct observations that restricted the ability to elucidate the complexity of root types and adaptation (e.g., excavation and trimming). This study describes the first direct and concurrent characterization of the root-soil current pathway and root architecture. Intrinsic differences in the root electrical properties of maize, cotton, grapevines, and orange trees emerge from the combined interpretation of our results with those presented by Mary et al. (2018Mary et al. ( , 2019, and Mary et al. 2020).
The laboratory rhizotron setup offered significant advantages for the development of the iCSD approach. Such advantages include the attainment of optimal electric coverage, direct root observation, and prior information on the growing medium. However, the development of the iCSD approach is well-timed and in line with the development of the geoelectrical technology for field-scale agricultural applications. The technological development driven by the use of ERT and SIP can readily be transferred to the iCSD approach. For example, the increased use of borehole electrodes and resistivity meter capabilities allow high electric coverage at the field-scale. The implementation of long-term monitoring has also become more common and provides precious information to characterize the dynamics of soil and vegetation. In this sense, the iCSD is a noninvasive solution that could be used to monitor the changes in root current pathways in response to water show the plant and root architecture at the moment of the experiment: the roots reach the bottom and sides of the rhizotron. c) iCSD result, the current passes from roots to soil in the most proximal part of the root system, with significant contribution from the stem section in the soil stress and other environmental factors. In conclusion, our inversion approach fills a major gap toward the full implementation of methods that allow non-destructive imaging of roots in their natural soil environment.


extfieldsofstudy: Medicine Chemistry

s2fieldsofstudy:
Agricultural And Food Sciences
Chemistry

provenance: AgriculturalAndFoodSciences,2019-2020-Agricultural And Food Sciences-2020.jsonl.gz.gz

text:
Contamination of Animal Feed with Undeclared Tetracyclines—Confirmatory Analysis by Liquid Chromatography–Mass Spectrometry after Microbiological Plate Test

The presence of tetracycline (TC) antibiotics was determined in animal feed that had been previously screened with a microbiological plate test. Feed samples were screened by a microbiological plate test on a pH 6.0 culture medium seeded with Bacillus cereus ATCC 11778 able to pre-reveal the presence of tetracyclines. Subsequently, confirmation and quantification were performed using a validated HPLC method with mass spectrometric detection. In 2013–2018, 353 feed samples were analysed to detect antibacterial substances, of which 186 (52.7%) were suspected to contain tetracyclines. Forty-two out of 186 (22.6%) samples analysed by the chromatographic method contained undeclared tetracyclines, which were determined at concentrations from 0.3 to 49 mg kg−1. The most frequently identified contaminating tetracyclines were doxycycline and chlortetracycline.

Introduction
Tetracyclines (TCs) are antibiotics that exhibit activity against infections caused by both Gram-positive and Gram-negative bacteria as well as chlamydia, mycoplasmas, protozoan parasites and rickettsiae [1][2][3]. They are widely administered to intensive farming animals because of their broad spectrum of activity and cost-effectiveness [1,4]. Since their discovery in the late 1940s, TCs have been used widely for the treatment and prophylaxis of bacterial diseases and as feed additives to promote growth in farm animals [5][6][7]. The emergence of bacterial resistance against tetracyclines, and the emergence of bacterial resistance against antibacterial agents in general, has caused severe drawbacks in their use. Concerns were raised on the potential problems due to the use of human antibiotics for growth promotion in animal husbandry and the transference of resistance to humans [8][9][10]. The discussion concluded in 1974 with the elimination of TCs as feed additives along with other antibiotics used in humans [11].
Several properties make tetracyclines the most commonly used of all antibiotics: they are active against most common pathogens, show good oral absorption, exhibit low toxicity and cause only few allergic reactions, and they are relatively inexpensive [12]. This has led to an intensive use of TCs for the treatment and prophylactic control of bacterial infections not only in humans, but also in animals. Of the eight forms of commercially available tetracyclines, four are frequently used in animal husbandry: chlortetracycline (CTC), doxycycline (DC), oxytetracycline (OTC) and tetracycline (TC). According to data described by Grave et al. [13] on the consumption of veterinary antimicrobial agents intended for food-producing animals, in all European countries tetracyclines, penicillins and sulfonamides accounted for more than half (range 53-88%) of the total amount of antimicrobial agents sold by country, expressed as mg/PCU. This trend continues and, according to the ESVAC report [14] in 2017, the highest sales concerned tetracyclines, followed by penicillins and sulfonamides, accounting for 30.5%, 26.9% and 9.8%, respectively of the total sales of veterinary antimicrobial agents.
The usage of antimicrobials in farm animals is strictly regulated in the European Union (EU) to protect consumers, as the excessive and inappropriate use of antibiotics can lead to the emergence of bacteria resistant to antibiotics or the presence of antibacterial residues on food products of animal origin [15][16][17]. In the EU, antibacterial substances permitted as feed additives to help fatten livestock were prohibited for non-medicinal purposes in 2006. The control and analysis of antibacterial substances in feed for farm animals has become an important issue, as only authorized feeds can be medicated under specific conditions, as stated in Regulation (EU) 2019/4 [18]. The antibacterial substances, in addition to medicated feeds, can enter feed via different paths, for instance due to unintentional cross-contamination in the feed mills or by unauthorized use for prophylaxis or to fraudulently promote growth [19]. Cross-contamination between medicated and nonmedicated feed can occur with any type of drug added to the feed, particularly when the cleaning process between batches is inefficient [20]. The fact that antibiotics could be present as contaminants in feed without farmers' knowledge implies that a withdrawal period will not be considered and, in addition to the development of antibiotic-resistant bacteria, antibiotic residues could remain in animal products (e.g., meat, eggs, milk) [21].
The effective and efficient control of the possible illegal use of undeclared antibiotics requires the availability of multi-screening and confirmation methods which can be implemented in an overall control strategy [22][23][24]. For the rapid detection of antibiotics in feed or food, microbiological assays are routinely used because they are easy to perform and inexpensive. However, the main disadvantage of these tests is their lack of specificity; that is, they often detect growth inhibition of a sensitive test strain, whether a substance with antimicrobial activity is present or not. The technique of choice for the confirmatory analysis of antimicrobial agents in feed is liquid chromatography coupled with mass spectrometry (MS). However, the analysis of animal feed has proved to be quite a challenge because of the high complexity and the variability of the matrix composition [25,26].
This study presents an antibiotics detection system based on a microbiological screening followed by a confirmation and quantification step consisting of a suitable chemical technique, such as HPLC with MS detection. This paper focuses on the analysis of the group of antibiotics approved by the EU for therapeutic use in pigs, poultry, cattle and sheep. The tetracyclines may be added to animal feed for treating the animals as a medicated feed. However, production practices in feed mills have been identified as the main source of feed contamination [27].
Because the use of antibiotics in feeds has been prohibited in the EU for non-medicinal purposes, monitoring for undeclared or illegal use of these substances is conducted within national programmes of feed control. This paper reports the investigation on the occurrence of tetracyclines in feed used in farm animals.

Results
Out of 353 feed samples, 186 were suspected for the presence of tetracyclines, because of the presence of an inhibition zone on the B. cereus ATCC 11778 plate. The LC−MS method was adopted for the confirmatory analysis of 186 feed samples found positive by the microbial assay. After LC−MS analysis, the presence of TCs was confirmed in 42 (22.6%) out of 186 suspect samples. Overall, the percentage of positive samples for TCs among all the analysed samples was 11.9%. The results are presented in Table 1. The most frequently identified antibiotic was doxycycline, followed by chlortetracycline ( Figure 1). Doxycycline was found in 20 samples, and chlortetracycline in 11 samples. Moreover, oxytetracycline was detected in 9 samples. Furthermore, in 2 samples the presence of doxycycline in addition to tylosin was revealed.  Table 2. Representative chromatograms of blank feed sample, fortified feed sample and real feed samples with chlortetracycline and doxycycline, analysed by the described method, are shown in

Discussion
Animal feeding has become an important economic activity and has a clear impact on food safety. Animal feeds, which are the main input into livestock production, not only require adequate quality from a nutritional point of view, but must also comply with legal limits regarding contaminants and be free from the unauthorized or unintentional presence of veterinary drugs. For this reason, the analysis and control of antibacterial agents in animal feed is a very important issue. Screening tests are often used to detect the presence of antimicrobials used in farm animals. The application of microbial assays for the screening of antibacterial substances in food or feed has been widely reported in the literature [15,[22][23][24][28][29][30]. Microbiological assays are screening methods that are among the most commonly used techniques for the detection of the majority of antibiotic classes, especially if a large number of samples are to be analysed. On the other hand, screening methods can have acceptable false-positive result rates. Having false-positive results is tolerable for a screening test that must be very sensitive, but not very selective. Hence, due to the risk of false-positive samples, results from microbiological assays require confirmation.
The eight-plate screening assay adopted in this study is based on the combination of pH conditions to favour or inhibit the antibiotic activity. Moreover, it is based on the sensitivity or resistance of the test microorganisms to various antibiotics. The combination of these factors enables a rough identification of antibiotics or antibiotic groups. In our method, eight plates used for a presumptive identification of six antibacterial families, including tetracyclines, proved to be a useful tool for the screening procedure. When all the plates are used, the pattern of inhibition can reduce unnecessary confirmatory analyses by indicating the antibacterial group most likely to be present. One of the eight plates is designed for the detection of tetracyclines owing to its low pH and a sensitive bacterial strain. Originally, the tetracycline family was the only one detected on the plate inoculated with B. cereus ATCC 11778.
To confirm the results obtained with the microbiological assay and to supply quantitative data, an HPLC analysis with mass-spectrometric detection was performed on the suspected and positive samples. Confirmation with HPLC or a similar technique is always necessary to identify and quantify analytes after a screening test. Due to the risk of false-positive samples, results from microbiological assays require confirmation by a confirmatory method, allowing for the selective, sensitive and accurate detection and quantification of antibiotics for effective surveillance.
Animal feed is a very complex matrix. Different animal species require different nutritive doses to get an equilibrated diet and be able to develop their specific productive functions. Poultry feeds in particular are characterized by a complex composition. In poultry, mainly commercial feeds are used; this lies in contrast to pigs, which are fed mostly production feeds. In general, feed composition may include, for example, amino acids, minerals, vitamins and salts, but also additives like organic acids, herbal mixtures or coccidiostats, or any other non-specific antibacterial substances (inhibitors) which induce matrix interference that may cause false-positive results in microbial assays. The lack of confirmed samples among poultry feeds as well as in concentrates/premixes shows that the presence of these interfering components can influence the occurrence of false-positive results. Nevertheless, microbiological methods play a crucial role in the effective official feed control and analyses.
The results showed that tetracycline antibiotics were identified in 42 (11.9%) of the analysed samples. Doxycycline and chlortetracycline were found in a total of 31 (73.8%) samples. The levels of antibiotics ranged from 0.3 to 5 mg kg −1 . Low levels of contamination suggest cross-contamination during manufacturing in the feed mill or on the farm, although in Canada pigs can receive feed supplemented with subtherapeutic doses of chlortetracycline, at 5.5 mg kg −1 [31]. However, the presence of oxytetracycline may be surprising, because in Poland there are no approved medicated premixes with this antibiotic to manufacture medicated feed; thus, it may suggest the illegal use of oxytetracycline on farms. This could be confirmed by quantified relatively high concentrations ranging from 0.9 to 49 mg kg −1 . Moreover, two feed samples contained doxycycline together with tylosin at concentrations of 7.9 and 10.9 mg kg −1 and 3.7 and 0.4 mg kg −1 , respectively, indicating either the illegal use of the antibiotics or cross-contamination. Cross-contamination during the production process in feed mills can be an important reason for the contamination of feed with undeclared drugs [27,30]. In the survey described by Lynas et al., the most frequently detected contaminant was chlortetracycline which was found in 15.2% of all analysed feedingstuffs. Furthermore, a study conducted in the Netherlands confirmed that nonmedicated feed batches were contaminated with antibiotic residues such as tetracyclines, penicillins and sulfonamides [20]. Out of 140 analysed samples, 87% of the samples contained detectable amounts of antibiotics, with the concentrations ranging from 0.1 to 154 mg kg −1 . However, our own research revealed the unconscious use of chlortetracycline-contaminated feed on few pig farms. The feed was manufactured by a medicated-feed-producing mill. The concentrations of contaminating chlortetracycline were between <0.3 and 2.7 mg kg −1 .
Usually, cross-contamination results in subtherapeutic concentrations in feeds which are unlikely to equal those prescribed for medication. In spite of this, they can sometimes be sufficient to cause violative residues in meat or meat products. Thus, the antibacterial substances that enter the feed via different paths can pass through the food chain to the consumers, causing possible allergic reactions, toxic effects and bacterial resistance problems. Feed contaminated with pharmacologically active substances can lead to the sale of food of animal origin containing drug residues, for which the zero level applies, or at concentrations greater than the authorized maximum residue limits MRLs. For example, Segato et al. [32] conducted studies in order to obtain information on the transfer of doxycycline to chicken tissues when feed contaminated by this antibiotic at the level of 4 mg kg −1 was administered to poultry. Under the adopted experimental conditions, the doxycycline concentrations in each tissue (muscle, liver, kidney) reached average levels lower than the corresponding MRL, but the analyte was present. This study showed the role played by feed; in the absence of the withdrawal period, contaminating drugs can pose a risk to the consumers due to the occurrence of residues.
Animal feeds must have the required quality, especially in the context of food chain safety. This is very significant in the event that antimicrobials might provoke allergies and contribute to the development of resistant bacterial strains if they reach the food chain. An exposure to low concentrations of antimicrobials may induce resistance in the normal gut bacteria of food animals, which could then be transferred to pathogenic organisms. Tetracyclines are one of the most commonly used groups of antibiotics in animal production. At the same time, they are categorized as highly important antimicrobials in human medicine [33]. The evidence for a link between non-human sources (animals, food, or the environment) and the potential to cause human disease is greatest for certain bacteria (e.g., non-typhoidal Salmonella, Campylobacter spp., Escherichia coli). Monitoring of antimicrobial resistance in commensal Escherichia coli isolated from slaughtered broilers, laying hens, turkeys, swine and cattle carried out in Poland between 2009-2012 revealed that antimicrobial resistance reached the highest values for tetracycline (43.3%), followed by ampicillin (42.3%) and ciprofloxacin (39.0%) [34]. Researchers consider that antibiotic use in the life cycle of food animals makes it more likely that zoonotic foodborne bacteria such as Salmonella and Campylobacter harboured by the animals will be resistant to common antibiotics-especially if they are used at subtherapeutic concentrations. An increase in the prevalence of tetracycline-resistant Salmonella or Campylobacter strains may have implications for the antibiotic treatment of salmonellosis or campylobacteriosis cases in humans. However, studies carried out by Holman [31] indicated that chlortetracycline appears to have no effect on the measurable chlortetracycline resistance when given at a concentration of 5.5 mg kg −1 feed.  Table 1.

Microbiological Plate Test
In the first stage, the tetracycline antibiotics were detected by a microbiological 8-plate screening method for detection of antibacterial substances in animal feedingstuffs, which was developed and validated in our laboratory [16,35]. Combination of agar medium, pH, and test microorganisms was optimised to yield a suitable sensitivity and optimal conditions for the detection of a single group of antibiotics. For tetracycline antibiotics, agar medium pH 6.0 and test microorganism Bacillus cereus ATCC 11778 were used. The limit of detection of this method for tetracyclines was 0.3 mg kg −1 .
To perform a sample assay, sterile petri dishes (diameter, 120 mm) were filled with 20 mL of inoculated medium and left to harden at room temperature on a horizontal surface. When the agar was solidified, two wells were punched with a cork borer of 11 mm diameter on the plate. The plates were subjected to a quality control when unknown samples were analysed. Two antimicrobial disks (tetracycline 30 µg, Oxoid) were placed on the control petri dish to confirm standardization.
A ten-gram portion of the sample was then weighed into 200 mL Erlenmeyer flask. The samples were extracted with 50 mL of methanol/phosphate buffer (pH 8.0) mixture (1:1, v/v). The extraction time was 30 min. After the extraction, the samples were centrifuged for 10 min at 3000 × g and 100 µL of the supernatants was pipetted into wells on the agar surface on the plates. Test plates were incubated overnight at 30 • C, and then inspected for inhibition zones around the wells. Results of the microbiological assay were obtained by analysing the effect of diffusion on the agar medium inoculated with bacterial strain Bacillus cereus ATCC 11778. The presence of antibacterial substances was shown by the formation of growth inhibition zones around the punch hole after overnight incubation, as shown on Figure 6. The diameter of the zones was measured and a positive result was defined by an inhibition zone greater than 14 mm. The samples showing growth inhibition zones on the test plate designed for TCs were suspicious of TCs and analysed and confirmed by liquid chromatography mass spectrometry.

Liquid Chromatography-Mass Spectrometry Analysis
An analytical method for simultaneous determination of oxytetracycline, tetracycline, chlortetracycline and doxycycline using liquid chromatography-mass spectrometry was developed and validated as described by Patyra and Kwiatek [36].
The HPLC-MS consisted of HPLC Agilent 1200 and MS Agilent 6140 (Agilent Technologies, Santa Clara, CA, USA) controlled by the software ChemStation (Agilent Technologies, Santa Clara, CA, USA). The chromatographic analyses were performer by injecting 15 µL of extract into a Kinetex C18 column (100 × 4.6 mm i.d., 2.6 µm) (Phenomenex, Torrance, CA, USA). Mobile phases A and B were mixed on a gradient mode with a flow rate of 0.5 ml min −1 . The column temperature was set at 20 • C. Analytes were detected in positive electrospray mode (ESI+). The monitored precursor ions for OTC, TC, CTC and DC were 461, 445, 479 and 445 m/z, respectively.

Chemicals, Reagents and Stock Solutions
Oxytetracycline, tetracycline, chlortetracycline, doxycycline, formic acid and ethylenediaminetetraacetic acid disodium salt (Na 2 EDTA) were purchased from Sigma-Aldrich (St. Louis, MO, USA). Acetonitrile and methanol (HPLC grade) were purchased from J.T. Baker (Deventer, the Netherlands). Sodium hydroxide and sodium hydrogen phosphate were obtained from Avantor Performance Materials (Gliwice, Poland) and citric acid from Acros Organics (Geel, Belgium). All reagents used were of analytical grade and analytically pure. Water was deionised (>18 MΩ cm −1 ) by the Milli-Q water purification system (Bedford, MA, USA).
The McIlvaine-Na 2 EDTA buffer was prepared by dissolving 11.406 g Na 2 EDTA in 115.65 mL 0.2 M phosphate buffer and 184.65 mL 0.1 M citric acid. The pH was adjusted to 4.0 by adding 0.1 M citric acid or 0.2 M phosphate buffer. Mobile phase A consisted of Milli-Q water acidified to 0.1% with formic acid, and mobile phase B consisted of acetonitrile acidified to 0.1% with formic acid.
To prepare the individual stock solutions of tetracyclines, 10 mg of tetracycline was dissolved in 10 mL of methanol and stored at −20 • C for up to six months. A working standard solution of tetracyclines was prepared by diluting the stock solution of each tetracycline to yield a final concentration of 100 µg mL −1 (the purity was consider when calculating appropriate concentrations) and stored at −20 • C for up to six months.

Sample Preparation and Clean-Up
Firstly, matrix-matched calibration curves were prepared. Different volumes of the tetracyclines working solutions were added to 5 g blank feed samples. Concentrations of tetracyclines in the matrix-matched feed samples were 0, 0.3, 0.5, 1.0, 2.5 and 5.0 mg kg −1 .
To extract the tetracyclines from the feed samples, 5 g of grounded feed and 25 mL of McIlvaine-Na 2 EDTA buffer were added to a 50 ml polypropylene tube. After shaking for 45 min on a horizontal shaker, the sample was centrifuged at 4000× g. Then, 5 mL of the supernatant was loaded onto a Strata-X-CW SPE cartridge (300 mg, 3 mL) that had been preconditioned with 3 mL of methanol and 3 mL of water. The cartridges were then rinsed with 6 mL of water, 6 mL of methanol and 3 mL of acetonitrile with 2% formic acid. The analytes were eluted with 3 mL of methanol with 2% formic acid. The eluate was evaporated to dryness under a nitrogen stream and the final residue was dissolved in 1 mL of 0.1% formic acid in water. Then, the extract was put into 2-mL autosampler vials and 15 µL were injected into the HPLC−MS system.

Conclusions
In conclusion, animal feed can be screened for the presence of antibacterial substances, including tetracyclines with a microbiological inhibition test, using a solid medium at pH 6.0 and B. cereus as a test organism. Confirmation with the HPLC technique is always necessary to gain information about the substance identification and quantification. Our study shows that tetracycline antibiotics may occur as contaminating agents in animal feed. A major implication of contamination of feed with undeclared antimicrobial agents is, firstly, the contribution to the development of bacteria resistant to drugs used to treat infections, and secondly, the risk concerning residues of antibiotics in edible tissues and products that can produce allergic or toxic reactions in consumers, particularly if it concerns antibiotics that are absorbed in the digestive tract [37,38].

Conflicts of Interest:
The authors declare that there is no conflict of interest regarding the publication of this article.