Functional Genomics Insights Into the Pathogenicity, Habitat Fitness, and Mechanisms Modifying Plant Development of Rhodococcus sp. PBTS1 and PBTS2

Pistachio Bushy Top Syndrome (PBTS) is a recently emerged disease that has strongly impacted the pistachio industry in California, Arizona, and New Mexico. The disease is caused by two bacteria, designated PBTS1 that is related to Rhodococcus corynebacterioides and PBTS2 that belongs to the species R. fascians. Here, we assessed the pathogenic character of the causative agents and examined their chromosomal sequences to predict the presence of particular functions that might contribute to the observed co-occurrence and their effect on plant hosts. In diverse assays, we confirmed the pathogenicity of the strains on “UCB-1” pistachio rootstock and showed that they can also impact the development of tobacco species, but concurrently inconsistencies in the ability to induce symptoms were revealed. We additionally evidence that fas genes are present only in a subpopulation of pure PBTS1 and PBTS2 cultures after growth on synthetic media, that these genes are easily lost upon cultivation in rich media, and that they are enriched for in an in planta environment. Analysis of the chromosomal sequences indicated that PBTS1 and PBTS2 might have complementary activities that would support niche partitioning. Growth experiments showed that the nutrient utilization pattern of both PBTS bacteria was not identical, thus avoiding co-inhabitant competition. PBTS2 appeared to have the potential to positively affect the habitat fitness of PBTS1 by improving its resistance against increased concentrations of copper and penicillins. Finally, mining the chromosomes of PBTS1 and PBTS2 suggested that the bacteria could produce cytokinins, auxins, and plant growth-stimulating volatiles and that PBTS2 might interfere with ethylene levels, in support of their impact on plant development. Subsequent experimentation supported these in silico predictions. Altogether, our data provide an explanation for the observed pathogenic behavior and unveil part of the strategies used by PBTS1 and PBTS2 to interact with plants.


INTRODUCTION
The cosmopolitan Actinobacteria, one of the largest phyla within the Bacteria, are acclaimed for their metabolic versatility and their capacity to produce biologically active agents with applications in medicine, agriculture, and biotechnology (Barka et al., 2016;Lee et al., 2018). Within the Actinobacteria, members of the genus Rhodococcus have been studied extensively because of their extraordinary ability to degrade a wide spectrum of xenobiotic and organic compounds, targeting them for bioremediation. Additionally, their often large genomes (up to 10 Mb), encode highly specialized pathways that are extremely useful for biotransformation and biocatalysis purposes (Ceniceros et al., 2017;Elsayed et al., 2017;Alvarez, 2019;Gupta et al., 2019). Nevertheless, besides these valuable characteristics, some Rhodococcus species use their metabolic capacities to infect animals, humans, and plants (von Bargen and Haas, 2009;Stes et al., 2011bStes et al., , 2013Anastasi et al., 2016).
Until recently, the capacity among rhodococci to interact with plants seemed to be restricted to the plant-pathogenic species R. fascians that is now considered to be part of the "Rhodococcus fascians assemblage" (Sangal et al., 2019). R. fascians isolates typically cause diverse developmental alterations in their hosts, including excessive shoot formation, leaf deformations, and stunted growth, collectively designated the leafy gall syndrome (Stes et al., 2011b). The host range of these bacteria encompasses a broad spectrum of mostly herbaceous dicotyledonous plants, although some woody plants and few monocotyledons are sensitive as well (Putnam and Miller, 2007;Depuydt et al., 2008b). Key to symptom development is the bacterial production of an array of cytokinins via gene products of the fasciation (fas) operon encoded on a linear virulence plasmid in most pathogenic isolates (Crespi et al., 1992;Pertry et al., 2009Pertry et al., , 2010Francis et al., 2012;Creason et al., 2014;Radhika et al., 2015;Jameson et al., 2019) and the concomitant modification of the hormone landscape of the infected plant (Depuydt et al., 2008a(Depuydt et al., , 2009aStes et al., 2011aStes et al., , 2012Stes et al., , 2015Jameson et al., 2019). Additionally, the pathogenic model strain D188 possesses a chromosomal locus vic (virulence in chromosome) implicated in symptom persistence that encodes a malate synthase believed to be involved in the catabolism of leafy gall-specific nutrients (Vereecke et al., 2002). Despite the absolute requirement of the linear plasmid for virulence in most leafy gall inducers, an important contribution of co-option of chromosomal and plasmid-encoded functions seems to guarantee a successful interaction with plants (Creason et al., 2014;Francis et al., 2016). As a result of metagenomic studies in search of biostimulants, it has now become evident that diverse Rhodococcus species are established members of the plant microbiome (Francis and Vereecke, 2019). Additionally, a novel disease on "UCB-1" pistachio rootstock trees, designated Pistachio Bushy Top Syndrome (PBTS) has been associated with the presence of two Rhodococcus species that act synergistically on their host (Stamler et al., 2015a,b). So, R. fascians is no longer the sole plant interactor nor the unique plant pathogen within the Rhodococcus genus.
PBTS is characterized by an array of symptoms, including stunted growth, formation of additional shoots, and a disturbed root development. Furthermore, affected "UCB-1" plants exhibit a strongly reduced capacity to accept a Pistacia vera graft and in the unlikely event that grafting is successful, cracked gall-like tissue carrying ectopic shoots develops at the graft junctions with a weakened stock-scion union as a consequence. Altogether, these disease aspects prevent the plants from going into production (Stamler et al., 2015b). The causative agents of PBTS were designated Rhodococcus isolate 1 and Rhodococcus isolate 2 (Stamler et al., 2015b), but here will be referred to as PBTS1 and PBTS2, respectively. 16S rRNA gene sequence analysis of several PBTS1 isolates obtained from symptomatic "UCB-1" trees, showed that they cluster together with R. corynebacterioides. Although R. corynebacterioides strains have not been studied extensively, they have been reported as endophytes from the xylem of eggplant and the leaves of Arabidopsis thaliana (Traw et al., 2007;Achari and Ramesh, 2014), inhabitants of the phyllosphere of apples (He et al., 2012), and as efficient degraders of aflatoxin (formerly Flavobacterium aurantiacum; Teniola et al., 2005;Risa et al., 2018), oil (Bayat et al., 2015), and rubber products (Pan et al., 2009). They have also been described as the cause of various infections in patients (Al Akhrass et al., 2012;Kitamura et al., 2012;Vergidis et al., 2017;Khalil et al., 2019) and as cancer killers exhibiting potent anti-malignancy activity (Zhou et al., 2017). In contrast, PBTS2 isolates grouped together with R. fascians (Stamler et al., 2015b). In agreement with the similarities between the symptoms associated with the leafy gall and PBTS syndromes, in both PBTS1 and PBTS2, the presence was demonstrated by PCR of the fasD gene, encoding an isopentenyl transferase that mediates the rate-limiting step of cytokinin biosynthesis (Crespi et al., 1992), the fasA gene coding for a P450 monooxygenase involved in zeatin production (Pertry et al., 2010; Genbank accession MF278335 for PBTS2 fasA), and vic (Stamler et al., 2015b). Although the presence of fasD is 100% correlated with virulence in leafy gall inducers (e.g., Crespi et al., 1992;Stange et al., 1996;Nikolaeva et al., 2009;Creason et al., 2014), its role in the PBTS strains remains to be determined. Moreover, the synergistic activity of both species in symptom development and the absence of genuine leafy galls in the PBTS syndrome suggest that the strategies used by these bacteria differ from those of leafy gall inducers.
Unexpectedly, the fas genes could not be detected in the assembled draft genome sequences of either of the PBTS strains (Stamler et al., 2016;Randall et al., 2018) and as a consequence the pathogenic status of PBTS1 and PBTS2 has been questioned (Savory et al., 2017). Here, we addressed this controversy and used a functional genomics approach to find clues as to how the PBTS bacteria would co-occur and modulate plant development. First, we analyzed the features of the available draft genomic sequences of PBTS1 and PBTS2 and compared them to each other and to those of R. fascians D188 (Stamler et al., 2016). We reassessed the pathogenicity of the PBTS strains on "UCB-1" plants and on two Nicotiana species and tested the hypothesis that the fas genes in the PBTS strains are easily lost during subculturing in rich medium, possibly the reason for their absence from the genome sequences (Randall et al., 2018;Vereecke, 2018). Additionally, based on the consistent co-isolation of both species from "UCB-1" trees exhibiting PBTS symptoms (Stamler et al., 2015a,b), we postulated that both isolates would not compete with each other for resources and would possibly exhibit complementary activities that would allow them to grow well in each other's presence which might be a prerequisite to synergize for optimal symptom formation (Stamler et al., 2015b). To this end, we mined the chromosomes of PBTS1 and PBTS2 to predict the occurrence of particular metabolic activities and phenotypes that might be important to thrive in their habitats, to associate with each other, and to modify plant development.
For the analysis of the genome similarities of the PBTS isolates and D188, chromosome and plasmid sequences were aligned by progressive MAUVE (Darling et al., 2004). For the comparison of the six PBTS isolates, D188, and R. corynebacterioides NBRC 14404, pair-wise Average Nucleotide Identity (ANI) values were calculated with nucmer 2 and the in silico DNA-DNA hybridization (isDDH) values by means of the Genometo-Genome Distance Calculator (GGDC). The genomic DNA sequences were uploaded into the GGDC 2.1 Web server 3 (Meier-Kolthoff et al., 2013). The point estimates from Formula 2 were utilized as the isDDH estimates.
Orthologous groups of the three reference genomes were determined with OrthoMCL v. 1.4 with the default parameters (Li et al., 2003). First, the OrthoMCL program carried out reciprocal comparisons of each protein sequence with BLASTp, where after the generated reciprocal BLASTp e-values were used to create a matrix analyzed by a Markov cluster (MCL) algorithm. Based on this analysis, OrthoMCL detected orthologous and paralogous genes of the three genomes and clustered them into orthologous groups.
For functional annotation and metabolic pathway reconstruction, a Cluster of Orthologous Genes (COG) annotation was conducted by a BLASTp search (threshold e-value 10 −5 ) against the COG myva database 4 . The Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation was done with BlastKOALA 5 (Kanehisa et al., 2016). Bacteria as target taxonomy group and the species_prokaryotes database were chosen in the scoring scheme for K number assignment. Proteins assigned with a valid K number were used to generate lists of reconstructed KEGG pathways with the KEGG Mapper 6 .
The bacterial strains were grown at 23-28 • C for 2 days on D2 (Kado and Heskett, 1970), yeast extract broth (YEB Broth, PhytoTechnology Laboratories, Overland Park, KS, United States), or nutrient broth (NB; Difco TM , ThermoFisher Scientific, Waltham, MA, United States) solidified with 1.5% agar (Difco TM Agar Technical). For infections, the bacteria were scraped from D2 plates and suspended in water to an OD 600 of 2.0, unless otherwise indicated.

Pathogenicity Assays on "UCB-1" and Nicotiana Species
Unbudded clonal "UCB-1" rootstock trees in plugs, 8 weeks acclimatized from tissue culture, were purchased from different commercial suppliers. Shoot infections were done by Agdia Inc (Elkhart, Indiana, United States). A 1/1 mixture of the PBTS bacteria [10 7 (colony-forming units) CFU/ml] resuspended in a buffer containing 10 mM MgCl 2 and 10 mM 2-(Nmorpholino)ethanesulfonic acid (MES) (pH 5.5) was applied on the leaves of 3 months old plants; control plants were mock-infected with water. The plants were placed in a growth chamber at 24 • C with a 12 h light/12-dark regime and their development was followed over several months during which the plants grew well. Images were taken 12 months after inoculation. Dip infections were done in the greenhouse facilities of the Research Institute for Agriculture, Fisheries and Food (ILVO-Plant, Melle, Belgium). The PBTS bacteria grown on D2 plates were scraped and resuspended in a buffer containing 10 mM MgCl 2 and 10 mM MES (pH 5.5) (OD 600 of 0.3), and mixed in a 1:1 ratio. The aerial part of 25 plants, 5-10 cm in height, and the plugs were dipped in the bacterial suspension and planted immediately in 20 cm deep pots with a diameter of 16 cm (volume 3 l); the substrate consisted of 60% peat moss, 35% perlite, and 5% vermiculite. At the base of the stem, 10 ml of the suspension was poured. As controls, 10 plants were treated in the same way but buffer was used as mock inoculum. The plants were placed in a greenhouse and to increase humidity, the planted material was covered by plastic for 1 week. Plant development was followed over a period of 4.5 months.
The pathogenicity assays with "UCB-1" seedlings were done at CSU Bakersfield. Uncoated "UCB-1" seeds were soaked in sterile water with aeration for 24 h, where after the water was replaced with 1% bleach and the seeds were soaked for another 24 h. Subsequently, the seeds were surface sterilized with 70% ethanol for 5 min, 20% bleach for 5 min, and rinsed three times with sterile water. The sterilized seeds were placed between sterile wet paper towels in sealed sterile containers, 5 seeds per MK-5 container (Caisson Laboratories, Smithfield, UT, United States) and incubated at 30 • C. From day 5 after the surface sterilization, the seeds were checked daily for germination, upon which they were placed right below the surface of moist sterile potting mixture. Plant growth was continued at 30 • C in a plant growth chamber with a 16 h light/8 h dark photoperiod. Three weeks old seedlings were infected by placing 20 µl drops of bacterial suspensions (prepared as described above) on each leaflet and axil; water was used as mockinfected control. Symptom development was monitored over a period of 30 days. The experiment was repeated twice with three plants for each condition.
Pathogenicity assessment on N. tabacum W38 including the seedling infection assay and infection via dipping has been described previously (Maes et al., 2001). The former assay was repeated twice with 25 seedlings for each condition; the latter assay was done only once with three replicates of four plants for each condition.
For infections of N. benthamiana, seeds were first surface sterilized (agitation in 2% bleach for 2 min, in 70% ethanol for 2 min, and three rinses of 2 min in deionized water) and planted in sterile vermiculite. After germination, the plants were transplanted when approximately 5.5 cm tall into autoclaved LM-40 High Porosity Mix (Lambert, Rivière-Ouelle, Quebec, Canada) inoculated with 50 ml aqueous suspensions of either PBTS1 or PBTS2 cultures, grown in YEB and diluted to 10 7 CFU/ml; uninoculated potting mix was used as negative control. For each condition, five replicates of single plants were used.

Assessment of the Presence of fas and vic Genes
The primer sequences used for the detection of vicA, vicA2, 16S rRNA, ITS, fasC, fasD, fasE, and fasF either by quantitative (q)PCR or endpoint PCR are given in Supplementary Table S1 (Nikolaeva et al., 2012;Dhandapani et al., 2017;DeLorenzo and Moon, 2018).
qPCR was utilized to determine the amount of fasD DNA relative to DNA of chromosomal genes in PBTS1, PBTS2, and D188. PBTS1 has two homologs of the vicA gene that are amplified with different primer sets. In PBTS1, vicA2 is the chromosomal homolog equivalent to vicA in PBTS2 and D188, whereas the vicA homolog in PBTS1 should be considered as a pathogenicity gene. Genomic DNA (gDNA) was extracted from bacterial pellets using the MoBio PowerLyser PowerSoil DNA isolation kit (QIAGEN, Carlsbad, United States). The PCR reactions were conducted in triplicate on a CFX92 real-time Thermal Cycler (Bio-Rad, Hercules, CA, United States) using iTaq TM Universal SYBR R Green Supermix (BioRad), 100 ng of gDNA as template, and 0.2 µM of each forward and reverse primer for the following target genes vicA, vicA2 (PBTS1), and fasD (fasD-JR primers). Reaction conditions included an initial denaturation step for 3 min at 95 • C, followed by 40 cycles of 10 s at 95 • C and 30 s at 60 • C. A final melting curve was done from 65 • C until 95 • C at 0.5 s increments of 1 • C.
For a second independent, preliminary qPCR experiment, gDNA was extracted from three different 2 days old PBTS1-like and PBTS2-like isolates and from D188 cultures obtained after direct inoculation of glycerol stocks in 5 ml YEB. The 12 µl qPCR reaction mixture was made with GoTaq R qPCR Master Mix (Promega, Madison, WI, United States) and consisted of 6.25 µl reaction buffer, 0.625 µl 5 µM forward (FW) primer, 0.625 µl 5 µM reversed (RV) primer, 0.208 µl dye, 3.3 µl water, and 1 µl gDNA (10 ng). The 16S rRNA gene was used as a reference and the target genes were fasC (using the fasD-PJ primers that do not target fasD as reported; Dhandapani et al., 2017), fasE and fasF. Single runs for each DNA sample were done on a CFX96 TM Real-Time PCR detection system (Bio-Rad) as follows: 10 min at 95 • C, followed by 40 cycles of 20 s at 95 • C, 20 s at 58 • C and 20 s at 72 • C. After 90 s at 72 • C, the melt curve of the amplicons was determined by increasing the temperature with 1 • C every 5 s from 72 • C until 95 • C. The amplicons were sent for sequencing to LGC (Berlin, Germany).
The presence of the fasD gene was additionally analyzed with endpoint PCR. PBTS1 and PBTS2 were inoculated from glycerol stocks in 5 ml YEB medium and after 2 days of growth at 28 • C, a 1-ml sample was taken for gDNA extraction. The remainder of the cultures was diluted in 500 ml YEB and grown for another 2 days, where after another 1 ml sample was taken for gDNA extraction. The 25 µl PCR reaction mixtures were prepared as follows: 5 µl 5 × colorless GoTaq R reaction buffer (Promega), 0.5 µl dNTP mix (10 mM stock; Promega), 1 µl FW primer (10 µM stock), 1 µl RV primer (10 µM stock), 0.125 µl GoTaq R DNA Polymerase (Promega), 16 µl water, and 1 µl (10 ng) gDNA. As a reference the ITS region was used (Nikolaeva et al., 2009) and the target gene was amplified with the fasD-DV primers. gDNA from D188 was used as a positive control and a template-free reaction as a negative control. The PCR was run on a T100 TM Thermal Cycler (Bio-Rad) with the following settings: 2 min at 95 • C, followed by 35 cycles of 20 s at 95 • C, 20 s at 50 • C and 1 min at 72 • C, and a final extension of 5 min at 72 • C. The fasD amplicons were sent for sequencing to LGC (Berlin, Germany).
Genomes of Rhodo 10, Rhodo 11, Rhodo 12, and Rhodo 13 were sequenced at Oxford Genomics Center (University of Oxford, United Kingdom) on the HiSeq 4000 platform (PE150 reads; Illumina, San Diego, CA, United States), yielding at least 72 Gb of data in total or 750 Mb on average per sample (depending on the genome size, GC content, and DNA quality). The library was prepared with an in-house adapted protocol of the New England Biolabs (Ipswich, MA, United States) prep kit.

In planta Enrichment Experiment
Surface sterilized N. benthamiana seeds were directly added to pots with autoclaved LM-40 potting mix (Lambert) to allow germination. Three pots for each treatment were placed in the same glass irrigation dish and sub-irrigated for the duration of the experiment. After germination, plants were thinned to leave three plants per pot for a total of nine plants per treatment.
When N. benthamiana seedlings reached the peak of their cotyledon stage, PBTS1 and PBTS2 cultures were generated from glycerol stocks on solid nutrient agar (NA) medium and grown for 2 days. Two-µl aliquots of aqueous suspensions obtained from the plate cultures (OD 600 set to 0.7) were pipetted onto each of the two emerged leaves of each plant and 2 ml samples were taken for gDNA extraction and qPCR analysis of vic and fasD as target genes as described above. After 3 weeks, the plants from each treatment were cut above the soil line, added to individual conical tubes containing 5 ml of Dulbecco's phosphate-buffered saline, and shaken for 1 min. Then, 2.5 ml of the buffer was plated onto D2 medium to obtain epiphytic bacteria. Additionally, the plant tissues were surface sterilized with the seed-sterilization method described above, transferred to individual conical tubes with 3 ml deionized water, and ground into a green liquid, of which 500 µl was plated onto D2 medium to obtain endophytic bacteria. Individual bacteria from both sets of plates were grown on D2 medium for gDNA extraction, 300 ng of which was used for qPCR analysis of vicA and fasD as described above.

Carbon Source Utilization, Penicillin and Copper Resistance, and Antagonism Tests
Bacterial growth assays on different carbohydrates and amino acids as carbon sources and the generation of antibiograms were done as described (Francis et al., 2016). In the latter experiment, PBTS1 titers in the halo's were quantified by resuspending the bacteria present on a 5 mm diameter agar plug harvested just aside from the antibiotic-containing disc in 300 µl sterile water. Subsequent serial dilutions were plated on YEB and bacteria were counted 4 days later.
For the Cu resistance test, a plate diffusion assay was done as described (Hassen et al., 1998) with CuSO 4 .5H 2 O concentrations of 100, 150, and 200 mM. Prior to inoculation, bacterial titers were set at an OD 600 of 0.1. Bacterial sensitivity was scored by measuring the length of the growth in the inoculation streak (in mm) after 6 days of incubation at 25 • C. The occurrence of PBTS1 at the tip of the growth closest to the highest Cu concentration was evaluated with a stereomicroscope M165FC equipped with a DFC 310 FX color camera (Leica, Wetzlar, Germany).
For the antagonism test, liquid YEB cultures of PBTS1, PBTS2, D188 and Bacillus subtilis were grown for 2 days under gentle agitation at 28 • C. Fifty-µl samples (OD 600 of 1.0) were diluted in 5 ml of YEB medium with 0.5% agar cooled to 45 • C. This mixture was poured in plates containing solid YEB medium. As soon as the overlays solidified, 5 µl aliquots of the cultures (OD 600 of 1.0) were spotted on the plates. Antagonistic effects as evidenced by reduced growth of the overlay strain around the spotted strains was evaluated daily for 5 days of incubation at 28 • C.

Modulation of Plant Development Assays
Samples from asymptomatic "UCB1" trees and plants displaying PBTS symptoms were collected for morphological analyses. Fresh petiole, stem, and root tissues were dissected into 2.5% glutaraldehyde, 0.1 M imidazole buffer (pH 7.2) and fixed overnight at room temperature. Subsequently, the fixed samples were washed with 0.1 M imidazole for 30 min, stained with 0.1 M imidazole:OsO 4 for 2 h, washed and dehydrated to 100% ethanol, followed by a final 100% acetone step; all steps were done at room temperature. The samples were embedded into Spurrs resin and polymerized overnight at 60 • C. Semithin 250 µm sections obtained with a UC6 ultramicrotome (Leica), were stained with toluidine blue and analyzed with a stereoscope M165FC equipped with a DFC 310 FX CCD color camera (Leica).
For the evaluation of indole-3-acetic acid (IAA) degradation, 5 mM IAA was tested as carbon source (Scott et al., 2013;Belimov et al., 2014). PBTS1, PBTS2 and D188 cultures, grown at 28 • C in YEB for 2 days, were centrifuged for 2 min at 14,000 rpm, washed twice, and suspended in sterile water. These suspensions were used to inoculate MinA medium (Miller, 1972) to reach a start OD 600 of 0.01. Growth was measured after 2 and 4 days and compared to similar cultures grown in MinA medium without a carbon source or with sucrose as controls.
Auxin production was tested with a colorimetric plate assay using Salkowski reagent. PBTS1, PBTS2, and D188 cells of 2 days old YEB cultures were washed with water and resuspended in JM medium (Crespi et al., 1992) (OD 600 of 0.1) supplemented with 500 µg/ml tryptophan or 10 µg/ml tryptophol (indole-3-ethanol) and incubated at 28 • C for 2 days. Then, 100 µl of the cultures were mixed with 100 µl Salkowski reagent R1 consisting of 12 g/l FeCl 3 in 7.9 M H 2 SO 4 . After 30 min the pink coloration of the supernatant was measured at 540 nm (Glickman and Dessaux, 1995), neither tryptophan nor tryptophol reacted with the reagent. The production of indolic compounds was normalized to the OD 600 of the respective cultures and the auxin concentration was calculated using a standard curve based on an IAA concentration range between 0 and 25 µg evaluated with the same protocol. The experiments were done three times.
For the assessment of interference with ethylene accumulation, 1-amino-cyclopropane-1-carboxylic acid (ACC) was utilized as a nitrogen source as described (Francis et al., 2016).
Ethylene production by PBTS1, PBTS2, and D188 was evaluated with the Arabidopsis triple response assay in split plates (Guzmán and Ecker, 1990). Half of the plates were filled with 10 ml PDA, and the other half was filled with 10 ml halfstrength Murashige and Skoog ( 1 /2MS) medium, supplemented with 0.5 g/l MES monohydrate, 0.1 g/l myo-inositol, 20 g/l sucrose, and 7.0 g/l Phyto agar (P1003; Duchefa, Haarlem, The Netherlands) (pH 5.7). The bacterial strains were inoculated on the PDA, whereas 10 Arabidopsis Col-0 seeds, stratified in the dark for 4 days at 4 • C, were transferred to the MS medium. The center partition of the bipartite I-plates (90 × 14.2 mm; Plastiques Gosselin, Borre, France) ensured a physical separation of the bacterial cultures and the seeds, allowing only gaseous exchange via the headspace. In the positive control plates, the ethylene-releasing agent ethephon (Agrichim, Pescara, Italy) was applied either in the MS medium or applied on top of the PDA medium at a final concentration of 60 µM. A control treatment consisting of blank PDA was included as well. Plates were sealed with air-permeable plastic foil, exposed to 45 µmol/m 2 /s light for a period of 2 h at 22 • C to synchronize seed germination, and subsequently grown in darkness for the remaining time of the experiment. After 7 days, seedlings were evaluated with respect to hypocotyl length and apical hook formation. The assay was conducted once with one biological replicate of 10 plants each.
The effect of emitted volatiles on plant development was tested in split plates in which two compartments were filled with 10 ml PDA and the other two with the same volume of 1 /2MS medium without sucrose. PBTS1, PBTS2, or D188 bacteria were streak-inoculated from available master plates and allowed to grow for 2 days at 28 • C. Then, 2 weeks old N. tabacum W38 seedlings and 1 week old Arabidopsis Col-0 seedlings, axenically grown on 1 /2MS with 1% sucrose, were transferred to the 1 /2MS compartments (6 and 5 seedlings, respectively). The plates were sealed with air-permeable plastic foil and placed in a growth chamber at 22 • C under a 16 h light/8 h dark photoperiod (45 µmol/m 2 /s light irradiance from cool-white fluorescent tungsten tubes). The plants were evaluated after 11 days of growth: pictures were taken, the shoots were excised and their fresh weight measured, where after a second picture was taken to document the root system.

General Features and Comparative Genomics of the PBTS1 and PBTS2 Genomic Sequences
Previously, we have published the assembled draft genome sequences of PBTS1 and PBTS2 and the full genome of D188 obtained with the PacBio technology (Stamler et al., 2016). Here we analyzed the structural and functional features of these sequences (Table 1). PBTS1 has the smallest circular chromosome of the three strains consisting of 4,251,687 bp with 3,917 predicted protein-coding genes and the highest GC content of 70%; no plasmid sequences were detected. The circular chromosomes of PBTS2 and D188 have similar characteristics, 5,179,353 bp and 5,139,988 bp in size, with 4,784 and 4,747 predicted protein coding sequences, respectively, and a GC content of 64.7%. Members of the genus Rhodococcus have three to six copies of the ribosomal RNA (rrn) operon, with an average of four 7 . The number of rrn operons of the three strains is in line with these findings. The strains have a comparable number of transfer (t)RNAs, suggesting an equal ability to respond to favorable conditions and a similar efficiency to utilize resources (Lee et al., 2009).
Classification of the protein-coding genes into clusters of orthologous groups (COG) assigned 71.9% of the PBTS1, 75.4% of the PBTS2, and 73.9% of the D188 chromosomal sequences into functional classes. Overall, in the three strains, the COG distribution of the chromosome-coding sequences was equivalent. Regarding well-characterized classes, the categories of 7 http://rrndb.umms.med.umich.edu/ amino acid (E), carbohydrate (G), and inorganic (P) transport and metabolism, and the transcription category (K), had the highest representation ( Table 2). Alignment of the chromosome sequences showed that the PBTS2 and D188 chromosomes were highly co-linear
( Figure 1A). Although PBTS2 also shared a significant amount of genetic characteristics with PBTS1, the syntenic regions in this case were smaller and more scattered due to inversions and translocations ( Figure 1B). Based on the chromosome sequences, the core genome was calculated with OrthoMCL v. 1.4 (Li et al., 2003). The three strains shared 2,879 orthologous groups representing 76% of the PBTS1 proteome, 61% of that of PBTS2, and 63% of that of D188; PBTS2 and D188 additionally shared 1,414 orthologous groups, corresponding in total to a 92-93% overlap in their proteomes ( Figure 1C). Of the three strains, PBTS1 had the highest number of strain-specific coding sequences, 812 orthologous groups (870 genes; Supplementary Table S2) representing 21% of its chromosomal genes. Both R. fascians strains only had 6-8% of unique genes (381 for PBTS2 (Supplementary Table S3); 320 for D188). Considering putative differences in host range, namely woody versus herbaceous for PBTS and D188, respectively, we were particularly interested in the 43 orthologous groups that were conserved in both PBTS isolates, but lacking in D188. However, amongst the unique PBTS genes, none encoded obvious host range-related functions (Supplementary Table S4).
Besides the circular chromosome, PBTS2 and D188 have a circular plasmid that has been designated pD188 and is not associated with pathogenicity in the latter (Desomer et al., 1988). In D188, pD188 is 164,725 bp in size with a 64.3% GC content and putatively encodes 173 proteins, of which 39.9% were assigned to a COG category ( Table 1). The pD188-like plasmid in PBTS2 has 137,595 bp with a 64.7% GC content and 140 proteinencoding genes, of which 35% were classified in a COG category ( Table 1). Although the COG distribution in both plasmids is quite different, in both cases the most represented classes are involved in post-translational modification, protein turnover, chaperones (O), and inorganic transport and metabolism (P) ( Table 2). Alignment of both plasmid sequences and calculation of their shared coding sequences showed a limited co-linearity ( Figure 1D). Only 62 orthologous groups were conserved ( Figure 1E), among which genes that confer resistance to copper and cadmium (Supplementary Table S5).
Finally, D188 possesses additionally the virulence-associated linear plasmid pFiD188 (Francis et al., 2012) that is 198,917 bp in size with a 61.8% GC content and 184 protein-encoding genes ( Table 1). COG annotation classified 38% of the coding sequences with the classes involved in amino acid transport and metabolism (E), and secondary metabolites biosynthesis, transport and catabolism (Q) most strongly represented ( Table 2).

The Pathogenicity of PBTS1 and PBTS2 Is Linked With the Presence of fas Sequences
The pathogenicity of PBTS1 and PBTS2 was re-evaluated on "UCB-1" clonal plants and seedlings and was additionally tested on two tobacco species. "UCB-1" clonal plants were coinoculated with PBTS1 and PBTS2 cultures, placed in a plant incubator, and their development followed over a period of 12 months. In agreement with our previous findings (Stamler et al., 2015b), control plants grew as expected (Figure 2A), whereas PBTS-infected trees showed severe stunting (Figure 2A) and a loss of apical dominance, resulting in the outgrowth of deformed short branches with small distorted leaves ( Figure 2B). Unexpectedly, in an independent experiment, none of the 25 clonal "UCB-1" plants dip-infected with a mixture of PBTS1 and PBTS2 developed symptoms, although the plants exhibited vigorous growth (Supplementary Figure S1). Subsequent qPCR analysis demonstrated that the PBTS1 culture did not contain any of the tested fas genes. In PBTS2, the Cq value for the 16S rRNA gene, used as a reference, was approximately 14, whereas very high Cq values ranging between 36 and 39 were obtained for fasC, fasE, and fasF signifying that few cells in the culture had these genes. Similarly, whereas inoculation of "UCB-1" seedlings with either of the PBTS isolates resulted in stunted growth of the shoots and reduced development of the root system in one experiment (Figure 2C), no symptoms developed on any of the inoculated plants in a subsequent experiment with the same setup (data not shown); the presence of the fas genes was not evaluated in these experiments.
As a first assessment of putative host range restrictions, the PBTS strains were used to infect two tobacco species, which are known to be highly responsive to leafy gall-inducing R. fascians isolates. In seedling infection assays with Nicotiana tabacum, plant responses have been categorized in 3 classes that determine the pathogenicity level of the strain under study (Maes et al., 2001). Infection with D188, used as reference, completely blocked seedling development at the cotyledon stage and fully inhibited root formation (Figure 2D), consistent with previous findings (Crespi et al., 1992;Temmerman et al., 2000;Vereecke et al., 2000Vereecke et al., , 2002; nearly all plants fell in category 1 (Figure 2E). Interestingly, both PBTS isolates also caused developmental effects in the tobacco host, but, in contrast to D188, the results were more variable. Inoculation with PBTS1 resulted in an almost equal distribution of plants exhibiting an intermediate (category 2) or no (category 3) growth inhibition, with less than 3% of the plants represented in category 1 (Figure 2E). Besides their growth retardation, none of the shoots exhibited abnormal morphologies, but the root systems were significantly reduced when compared to non-infected controls ( Figure 2D). Upon infection with PBTS2, over 40% of the plants fell in category 1 and 28% in categories 2 and 3 ( Figure 2E). Interestingly, PBTS2 infection also resulted in leaf deformations in up to 60% of the plants and inhibited root formation almost to the same extent as D188 (Figure 2D). Determination of the fresh weight of the plants revealed that D188 decreased the biomass of the plants to 10% of that of the non-infected controls, whereas plants infected with either PBTS1 or PBTS2 reached a fresh weight that was approximately 20% of that of the uninfected controls ( Figure 2E). In line with these findings, when seedlings of N. benthamiana, which has also been used as a host for leafy gall inducers (Creason et al., 2014), were planted in substrate co-inoculated with PBTS1 and PBTS2, their growth was severely stunted, flowering was delayed, and the number of flowers was reduced ( Figure 2F).
Altogether these results confirm the pathogenic identity of PBTS1 and PBTS2 and extend their host range past pistachio. However, not all infections consistently provoked symptom development and the maintenance of fas genes in the inoculum was seemingly an issue that might be at the basis of the observed variability. This limited occurrence of fas genes in the PBTS cultures concurs with the absence of fasA and fasD in the genomic sequences of both PBTS1 and PBTS2 and the absence of the vicA homolog identified by PCR in the PBTS1 sequence (Stamler et al., 2015b). From what follows, inoculum preparation appears to have an important impact on the presence or absence of these putative virulence genes in the PBTS cultures and thus on the outcome of the infections.

Evidence for the Presence and Loss of pFi-Related Sequences in PBTS1 and PBTS2 Cultures
From the above and in contrast to the leafy gall-inducing model strain D188, maintenance of the fas genes does not seem to be a stable characteristic in the PBTS isolates. To get insight into this matter, we examined by qPCR the number of cells within pure cultures of PBTS1, PBTS2, and D188 that contain the chromosomal vicA gene (vicA2 in PBTS1) and the fasD gene that is located on the conjugative linear plasmid pFiD188 in strain D188 (Crespi et al., 1992). PBTS1 has two homologs of the malate synthase gene: vicA2 which is the chromosomal counterpart of vicA in PBTS2 and D188, and vicA which should be considered a pathogenicity gene in this isolate (Stamler et al., 2015b;Savory et al., 2017). The results indicated that the number of cells that contained vicA was comparable for the three strains (vicA2 in PBTS1), but that a significant discrepancy occurred amongst strains for the presence of cells positive for fasD (and vicA for PBTS1). Whereas the ratio for fasD-and vicA-positive cells was 1:1 in D188 cultures, for PBTS1 only 1 in 10,000 vicA2-positive cells contained the fasD and the vicA gene, and for PBTS2 only 1 in 10,000 cells contained both fasD and vicA. In an independent, indicative qPCR experiment, the 16S rRNA gene was used as a reference and the presence of three genes of the fas operon, fasC, fasE and fasF, was assayed in D188 cultures and in three different isolates of each PBTS species. The number of D188 cells that carried the three fas genes and the 16S rRNA gene was 100%, corroborating the presence of pFiD188 in each of the cells. For the PBTS1 isolates, however, on average only 8% of the cells had fasC (4.4, 8.8, and 12%), 4.9% fasE (1.5, 5.0, and 8.3%), and 2.8% fasF (1.6, 2.2, and 4.7%). For the PBTS2 isolates the number of cells with the fas genes was even lower ranging between 1.5 and 2.3 cells in 10,000 (fasC: 1.6, 1.8, and 3.3 in 10,000; fasE: 1.7, 1.5, and 4.6 in 10,000; fasF: 1.1, 1.2, and 2.1 in 10,000). Sequence determination of the vic and fas amplicons of both experiments confirmed their identity (Supplementary Figure S2).
The finding that only a subset of PBTS1 and PBTS2 cells retain fas genes when cultured was further supported by sequencing two additional PBTS1-like (Rhodo 11 and Rhodo 13) and PBTS2-like isolates (Rhodo 10 and Rhodo 12) with the Illumina technology. Calculation of pairwise ANI and isDDH values, revealed that the genomic sequences of the two PBTS1-like isolates were nearly identical to the published PBTS1 sequence, but very different from the genome sequence of R. corynebacterioides NBRC 14404 (Supplementary Figure S3). In fact, despite the 99.73% identity of their 16S rRNA sequences, the obtained ANI and isDDH values were well below the respective 94% and 70% threshold values determined for species delineation (Sangal et al., 2016), indicating that PBTS1 should not be considered a R. corynebacterioides strain. In contrast, whereas the genomes of the two PBTS2-like isolates were identical based on the ANI and isDDH values, they diverged more from the published PBTS2 and the D188 sequence, but were still well within the range to be considered the same species (Supplementary Figure S3); hence, PBTS2 should be considered a R. fascians strain. When the reads were assembled, no pFi-like replicon was identified in any of the four isolates, concurring with the published data (Stamler et al., 2016). However, by assuming that only a fraction of the cells in the cultures used for sequencing would contain virulence genes, the putative virulence replicon would probably not be assembled. In line with this reasoning, we mapped all the reads to pFiD188. For Rhodo 12 and Rhodo 13, a single hit was obtained, corresponding to pFi_052 and attG, respectively (Supplementary Table S6), and none for Rhodo 11. However, with Rhodo 10, 15 hits were obtained yielding sequences of 2,092 bp identical and 378 bp similar to pFiD188 (83-98% identical) (Supplementary Table S6). The hits were spread throughout the pFiD188 sequence and comprised three loci of the conserved R regions and eight loci of the unique U regions (Francis et al., 2012), including fasA, fasF, nrp2, stk4, and stk5 (Supplementary Table S6). Whereas, based on these data, no hypothesis can be put forward for PBTS1, we could speculate that PBTS2 has a virulence carrier resembling pFiD188, but this hypothesis requires further experimental confirmation. We then used the reads to calculate the number of cells with such a putative replicon in the cultures used for genomic sequencing: the obtained numbers were 1:200,000 cells (15 out of 3,324,630 pairs) for Rhodo 10, 1:4,000,000 cells (1 out 3,859,253 pairs) for Rhodo 13 and 1:3,000,000 for Rhodo 12 (1 out 31677999 pairs), which is even lower than the results obtained in the qPCR experiments.
From the above, we postulated that the fas genes in the PBTS isolates were lost upon repeated sub-culturing in rich medium. To test this hypothesis, 5 ml YEB cultures of PBTS1 and PBTS2 were initiated from glycerol stocks. After 2 days of growth at 28 • C, a 1-ml sample was taken for gDNA extraction and the remainder of the pre-culture was diluted in 500 ml YEB and grown for another 2 days, where after gDNA was extracted from a 1 ml sample. Ten ng of each gDNA was used for endpoint PCR analysis of fasD and the ITS sequence as a reference; gDNA from D188 was used as a control. A fragment of comparable intensity was detected for the ITS amplicon in all 5 cultures. In accordance with the previous experiments, a fasD amplicon was detected in PBTS1 and PBTS2 pre-cultures, but the intensity of the fragment was weaker than in D188, indicating that it is present only in a sub-fraction of the PBTS cells (Figure 3). The identity of the fasD amplicon was confirmed by sequencing (Supplementary Figure S2). Importantly, after a single additional passage in rich medium the presence of fasD could no longer be detected in both PBTS cultures (Figure 3).
We then hypothesized that the sub-population of fasDpositive cells in the PBTS1 and PBTS2 cultures would increase when the bacteria were in contact with a plant host. In support of this postulation, fasD and vic could be amplified by qPCR in gDNA extracted from leaf and stem samples of field grown "UCB-1" rootstock trees exhibiting PBTS symptoms (data not shown). Next, liquid PBTS cultures were used to inoculate N. benthamiana plants and a sample of these start cultures was taken for gDNA extraction. Three weeks after infection epiphytic colonies were isolated from plant tissue washes and endophytic colonies from surface sterilized plant tissue; gDNA was extracted from these samples. The gDNA from the starting cultures and from the re-isolated bacteria was subjected to PCR for the amplification of vic and fasD. For both PBTS1 and PBTS2, hardly any amplification of fasD was detected in the starting cultures, although the vic gene was clearly present in the populations. However, a single passage on the host resulted in a clear enrichment of fasD in epi-and endophytic cells (Figures 4A,B).
Altogether these data provide strong evidence that the maintenance of the putative virulence genes in the PBTS strains is unstable in an ex planta environment. Loss of virulence is actually a recurrent observation for all kinds of pathogens (e.g., Krokene FIGURE 3 | Presence and loss of fasD in PBTS1, PBTS2, and D188 cultures. PCR analysis with ITS amplification to validate DNA quality, equal concentration, and stability. Amplification of fasD shows its presence in the PBTS1 and PBTS2 pre-cultures (albeit in fewer cells than in D188), and its loss in the PBTS isolates after a single additional passage in rich liquid medium. λ PstI, DNA ladder; NT, no template control. The images of the original gels have been cut and rearranged for a more logic order of the samples. The original images have been provided to the reviewers for evaluation. and Solheim, 2001; Wright et al., 2003;Soto et al., 2006;Mitra, 2015) and although the presence of the linear virulence plasmid is a very solid characteristic in many leafy gall inducers, virulence instability has been recorded for R. fascians as well. Already in 1948, Margaret Lacey, a "leafy gall" research pioneer, reported that anti-auxin activity in some R. fascians isolates was lost by cultivation on artificial media (Lacey, 1948). Jacobs and Mohanty (1951) observed that several virulent R. fascians strains that had retained their pathogenicity at least 5-6 years after isolation, would become non-pathogenic after longer laboratory cultivation periods. Growth of R. fascians isolates on rich media at 37 • C resulted in a complete loss of virulence (Sabart et al., 1986). Even strain D188 could be cured from its linear plasmid by growth at 37 • C followed by other harsh treatments (Desomer et al., 1988). Additionally, plasmid instability has also been observed in R. fascians isolates in the absence of non-permissive growth conditions. For instance, Nikolaeva et al. (2009) showed that different isolates from symptomatic ornamental plants exhibited a varying degree of pathogenicity (±35% pathogenic, ±30% non-pathogenic, and ±35% with a variable phenotype). When single colonies from individual pure isolates exhibiting variable pathogenicity were analyzed, they turned out to consist of two sub-populations, one with fasD and one without fasD; bacteria without fasD lacked all tested linear plasmid markers (Nikolaeva et al., 2009). It would be interesting to compare the chromosomal and plasmid sequences of some of the unstable isolates to those of D188 to identify putative stability/instability determinants, but as far as we know, their genomes have not been sequenced. Intriguingly, as opposed to the loss of virulence upon growth outside the host, both Lacey (1948) and Sabart et al. (1986) noticed that re-isolated bacteria obtained from plants inoculated with weakly active (auxin degradation) or "regular" virulent R. fascians isolates, exhibited a strong increase in auxin degradation potential or a hypervirulent character, respectively. These findings are in line with our observation of in planta enrichment of fasD-positive PBTS isolates.
The instability of the fas genes in the PBTS isolates will make the study of their pathogenicity strategies very challenging. Nevertheless, pending the identification of their virulence carrier, we explored their chromosome sequences for genes that might be implicated in the adaptation to their habitats and their apparent co-occurrence or that may play a role in the interaction with their host. Most of these features were functionally assessed to validate the in silico predictions.

The Metabolic Features of PBTS1 and PBTS2 Suggest That Nutritional Resource Partitioning Probably Contributes to Their Apparent Co-occurrence
More than 20% of the coding sequences in the two PBTS isolates belong to the COG annotation class T "Signal transduction mechanisms" (  Table S9) indicating that they might benefit from the availability of various potential N and C sources in their habitat, including mannitol, ribose, glucose, betaine, choline, and diverse amino acids that are all released by plants as root exudates (Faure et al., 2009). KEGG analysis of the chromosome sequences of PBTS1 and PBTS2 also predicted the presence of several central metabolism pathways, including glycolysis, pyruvate metabolism, the Entner-Doudoroff pathway, the Krebs cycle, and the pentose phosphate pathway, implying that the PBTS bacteria are capable of using the wide variety of carbon sources internalized by the transporters. Interestingly, based on these analyses, PBTS1 and PBTS2 would be unable to utilize galactose and methionine and inositol, xylose, and citrate, respectively. Such a differential utilization of nutrient sources that are putatively available in the habitat might represent one mechanism to support the co-occurrence of the PBTS strains (Schlechter et al., 2019). This prediction was tested by evaluating the nutrient utilization patterns of PBTS1 and PBTS2. To this end, the PBTS bacteria were grown in minimal medium supplemented with a 0.5% final concentration of either one of 22 different carbohydrates or either one of 12 amino acids as carbon sources; strain D188 was included as reference. In agreement with the sequence data, PBTS1 was unable to grow on methionine or galactose and PBTS2 did not proliferate on myo-inositol. However, PBTS2 did grow on citrate ( Table 3). Since there were quite some differences in the carbon source utilization profiles of the three strains (Table 3), we calculated niche overlap indices (NOIs) to estimate niche similarity of the PBTS strains and D188. As described by Wilson and Lindow (1994), for each pair of strains NOIs were calculated as the number of carbon sources utilized by both strains as a proportion of the total number of carbon sources utilized by the strain in question. The smallest niche overlap was between PBTS2 and References The data are the average of two replicates measured two days after inoculation; bold characters indicate that the bacteria had not reached stationary phase at that time point and continued to grow until 4 days after inoculation when the experiment was ended. Growth is relative to the negative control without addition of any carbon source (−) and growth in rich YEB medium (++++). *ACC was given as nitrogen source, with sucrose as carbon source; growth on (NH 4 ) 2 SO 4 was used as a reference.
PBTS1 (Figure 5), suggesting that the nutritional competition between both strains might not interfere with their co-habitation of the same niche.
Altogether, these data indicate that the metabolic features of PBTS1 and PBTS2 would support their co-occurrence, in agreement with the consistent co-isolation of both PBTS strains from "UCB-1" trees exhibiting PBTS symptoms (Stamler et al., 2015a,b).

Complementary Activities in PBTS2 Might Support Habitat Fitness of PBTS1
Besides niche partitioning, complementary activities involved in habitat competitiveness might also contribute to the cooccurrence of both PBTS species. Therefore, the PBTS1 and PBTS2 chromosomes were mined for the presence of genes putatively involved in bactericide resistance and antimicrobial production. To mitigate copper micronutrient deficiencies in pistachio orchards, trees are fertigated or sprayed with copper sulfate (Kallsen et al., 2000;Soliemanzadeh et al., 2014). Additionally, copper fungicides are used for the management of Septoria leaf spot (Michailides, 2005). As a consequence of these management practices, epiphytic PBTS populations could be subjected to toxic levels of copper. The chromosome of PBTS2 contains nine genes putatively involved in copper detoxification and active copper efflux, including two copper oxidases, a copper translocating P-type ATPase, a CopC-like chaperone, a CopDlike copper resistance protein, two copper binding proteins, a copper transporter, and a CopY-like repressor. Even more, PBTS2 has an additional copper resistance system on the pD188-like circular plasmid that consists of a CopG-like copper translocating P-type ATPase, a CopC-like chaperone, a CopD-like copper resistance protein, and a CopY-like repressor (Supplementary Table S10). Although PBTS1 also contains six genes encoding proteins related to copper resistance, including three copper chaperones, one copper translocating P-type ATPase, a CopDlike copper resistance protein, and a copper transporter, no copper oxidase was identified suggesting that efflux is the main copper homeostasis mechanism in this strain (Supplementary Table S10). Copper resistance assayed with the plate diffusion method (Hassen et al., 1998), revealed that PBTS1 was indeed more sensitive to copper than PBTS2 as the extent of the growth in the inoculated streaks was consistently significantly less for PBTS1 than for PBTS2 ( Figure 6A). However, by inoculating PBTS1 together with PBTS2, PBTS1 colonies occurred at the streak tips at all copper concentrations tested, indicating that the presence of PBTS2 increases the PBTS1 tolerance to higher Cu levels ( Figure 6A). Furthermore, PBTS2 possesses all the genes required for 4-hydroxybenzoate degradation, a natural antimicrobial phenolic compound produced by many plants and highly abundant in ripe pistachio hulls (Barreca et al., 2016). The pathway is organized in a more complicated superoperonic gene cluster (pobAB, pcaJI, pcaHGBDR, pcaC, and pcaF) that differs from that described in Xanthomonas campestris pv. campestris (Wang et al., 2015), but more importantly, in PBTS1 this detoxification mechanism is seemingly incomplete (Supplementary Table S10).
AntiSMASH analysis (Blin et al., 2019) suggested that the spectrum of bioactive compounds produced by the secondary metabolite clusters of both PBTS isolates is probably quite FIGURE 5 | Calculated niche overlap indices based on the utilization of diverse carbohydrates, amino acids, and hormones as carbon sources and ACC as nitrogen source. SUI, strain under investigation; CIS, co-inhabiting strain. Table S11), implying putative complementarity in their responses toward surrounding microbes. Following an experiment in which 3 weeks old axenically grown N. tabacum plants were dip-infected with a mixture of PBTS1 and PBTS2, the excessive growth of both bacteria observed on the plant medium clearly indicated that they did not exhibit any antagonistic effect against each other ( Figure 6B). Nevertheless, when PBTS1 was co-inoculated with D188, hardly any orange-colored bacteria could be detected, suggesting that D188 produced a bioactive compound that is not synthesized by the highly similar PBTS2 (Figure 6B). Accordingly, no antagonistic activity of PBTS2 could be demonstrated against PBTS1 nor against D188 with overlay assays. When PBTS1 was used as overlay, however, the cells grew less densely around the D188 spot, but no halo developed around PBTS2 (Figure 6C). Comparison of the antiSMASH results on the PBTS2 and D188 chromosomes revealed an almost complete overlap in secondary metabolite clusters (Supplementary Table S11), with two exceptions. Cluster 12 in D188 contains a single putative type I PKS-encoding gene with an upstream gene encoding a hypothetical protein and a downstream gene encoding an aminotransferase class I/II-fold pyridoxal phosphate-dependent enzyme (A3L23_RS14450, A3L23_RS14445, and A3L23_RS 14440); these three genes are absent in the PBTS2 chromosome. Additionally, although cluster 3 in PBTS2 and cluster 9 in D188 are similar and putatively code for the production of a bacteriocin, PBTS2 has two extra genes and a small hypothetical gene that differs from that in D188, possibly affecting functionality (Supplementary Table S11). Further experimentation is required to determine whether these gene clusters are responsible for the observed antagonistic effect of D188 against PBTS1.

different (Supplementary
Soil dwelling microbes typically use antibiotic production, especially ß-lactams, to protect their niche in this nutrient-poor environment. Therefore, successful soil inhabitants often exhibit antibiotic resistances (Butterworth et al., 1979;Bucher et al., 2019). Both PBTS strains appeared to be characterized by a differential set of ß-lactamases (Supplementary Table S10) and, in agreement with the higher number of putative ß-lactamase genes in PBTS2 compared to PBTS1 (5 versus 2), antibiograms showed that the former was indeed much more resistant to different penicillins (Table 4). Quantification by serial dilutions of PBTS1 titers just aside antibiotic discs in antibiograms with four penicillins in overlays with either PBTS1 alone or both PBTS strains together, showed that the presence of PBTS2 enhanced the tolerance of PBTS1 for increased penicillin levels, with PBTS1 titers that were 10-fold higher for oxacillin, 14-fold for mezlocillin, and 20-fold for amoxycillin ( Figure 6D). Taken together these data indicate that by associating with PBTS2, PBTS1 would strongly improve its ability to survive in soil and on pistachio plants.

The PBTS1 and PBTS2 Chromosomes Encode Different Iron-Capturing Mechanisms
Plant surfaces as well as soil are environments that are highly challenging and successful survival requires adaptive skills to acquire essential nutrients from the surroundings. As iron is a growth-limiting factor, plant pathogens generally have diverse mechanisms to get iron from either their host or their habitat (Lemanceau et al., 2009). As for other bacteria (Sheldon and Heinrichs, 2015), the presence of two iron-dependent repressors indicates that the iron regulon in both PBTS strains is under tight transcriptional control (Supplementary Table S10). Although the bioavailability of inorganic iron is generally low, the PBTS isolates have an inorganic iron transport system that allows them to access ferrous iron under permissive environmental conditions (Supplementary Table S10; Sheldon and Heinrichs, 2015).
AntiSMASH analysis of the chromosomes of PBTS1 and PBTS2 identified 13 and 19 regions possibly involved in secondary metabolism, respectively (Supplementary Table S11). In PBTS2, cluster 9 was highly similar to the siderophore operon previously identified in the D188 chromosome (Francis et al., 2016). This operon consists of 17 genes encompassing three non-ribosomal peptide synthases (NRPSs), a salicylate synthase, three siderophore-interacting proteins, and two composite ABC transporters (Supplementary Table S10) that are related to enzymes involved in mycobactin production, a mycobacterial salicylic acid-based siderophore essential for pathogenicity (De Voss et al., 2000;Reddy et al., 2013). In agreement with this finding, siderophore activity could indeed be demonstrated on Chrome Azurol S (CAS)/YEB plates ( Figure 6E).
Regarding siderophore production by PBTS1, the situation is less clear. Although siderophore activity was detected (Figure 6E), none of the six NRPS-containing loci nor the 2 polyketide synthase (PKS) regions predicted by the antiSMASH analysis contained genes with homologies to siderophorerelated functions (Supplementary Table S11). The antiSMASH analysis indicated that cluster 10 was 18% similar to a gene cluster of Streptomyces sp. AA4 involved in the production of the mixed-ligand siderophore amychelin (Seyedsayamdost et al., 2011), but the overlap was too low to predict this type of siderophore in PBTS1. Besides NRPS-based siderophore synthesis, a second mechanism termed NRPS-independent siderophore (NIS) biosynthesis is often deployed (Sheldon and Heinrichs, 2015), but BLAST searches with a variety of NIS-related genes from different organisms did not yield any significant hit with the PBTS1 chromosome. Therefore, the locus FIGURE 6 | Functional evidence for competitiveness features of PBTS1 and PBTS2. (A) Copper resistance assayed with the plate diffusion method. The top panels represent an overview of the bacterial growth on the different CuSO 4 concentrations and a detailed view of the tip of the growth streak near the copper source (note the presence of the orange-colored PBTS1 colonies). The table represents the average length of the growth streaks. At each copper concentration, the growth of PBTS1 in the PBTS1 and PBTS1+PBTS2 inoculation streaks was statistically different as calculated by Student's t-tests (P < 0.01; n = 6). (B,C) Antagonistic effect of D188, but not of PBTS2, against PBTS1. (B) Differential growth of PBTS1 on 1 /2MS medium after co-inoculation of N. tabacum with PBTS2 (top) or D188 (bottom). (C) Lack of antagonistic activity by any of the strains when an overlay of PBTS2 is applied (top). In case of an overlay with PBTS1, an inhibition zone is formed around D188, but not around PBTS2 (bottom); Bacillus subtilis was used as positive control. (D) Antibiograms with four penicillins (1, mezlocillin; 2, oxacillin; 3, amoxicillin; and 4, ampicillin) and overlays of PBTS1, PBTS2 or both PBTS strains together. After 3 days of growth, the PBTS1 titers on an agar plug adjacent to the antibiotic discs were determined (see inset for sampling zone). The average PBTS1 titer is indicated and was statistically different between PBTS1 and PBTS1+PBTS2 for each antibiotic (calculated using Student's t-tests; P < 0.01; n = 6). (E) Siderophore production by PBTS1 and PBTS2 on CAS/YEB plates.
involved in the observed siderophore activity of this strain remains elusive.
Besides the capacity to produce their own siderophore, both PBTS1 and PBTS2 seem to be well-equipped to obtain siderophores from their environment ( Supplementary  Table S10). Indeed, PBTS1 has three siderophore-interacting proteins and three solitary iron-siderophore ABC transporter substrate-binding proteins, of which some might interact with the permeases of the ferric-chelate uptake system encoded in the PBTS1 chromosome. Additionally, this strain has a composite ABC transporter putatively involved in heme-iron acquisition, a putative heme binding protein, and a heme monooxygenase. Altogether, upon heme internalization, PBTS1 might be able to degrade it to release iron (Sheldon and Heinrichs, 2015). Similarly, PBTS2 contains two siderophore-interacting proteins, seven solitary iron-siderophore ABC transporter substratebinding proteins, two ferric-chelate and one ferric-citrate uptake systems, and four iron acquisition composite ABC transporters. Although PBTS2 has one putative heme-iron uptake system, the heme-binding protein and the heme monooxygenase are absent. It thus appears that both PBTS isolates invest more energy in parasitizing on xeno-siderophores produced by other organisms, than in producing and secreting their own siderophores. This mode of action has been suggested to be more energy efficient and to provide a selective advantage over other bacteria in a heterogeneous population (Lemanceau et al., 2009;Sheldon and Heinrichs, 2015).

PBTS1 and PBTS2 Have the Potential to Modulate the Hormone Homeostasis and the Development of Their Host
When the developmental changes that occur in "UCB-1" plants infected with the PBTS bacteria are analyzed, symptoms, such as the apical dominance loss and stunted growth of both shoots and roots, imply that the hormone status of the plant is strongly impacted upon colonization. In D188, the main virulence factor that affects plant development is the production of cytokinins through the linear plasmid-encoded fas pathway (Stes et al., 2011b). Additionally, basal levels of cytokinins are also produced by the chromosomally encoded tRNA degradation pathway (Francis et al., 2016), but the role of these cytokinins in the leafy gall pathology remains to be determined. Mining of the PBTS1 and PBTS2 chromosomes revealed the presence of an intact methylerythritol phosphate (MEP) pathway that provides the required cytokinin precursor molecule dimethylallyl pyrophosphate (DMAPP), a miaA and miaB gene encoding a tRNA-isopentenyltransferase and a tRNA methylthiolase, respectively, and a LONELEY GUY homolog encoding a phosphoribohydrolase activating cytokinin nucleotides (Supplementary Table S12). Thus, just as for the leafy gall inducer D188, the PBTS isolates are likely capable to produce cytokinins via the tRNA degradation pathway (Supplementary Figure S4). Although the relevance of this pathway for plant-interacting microbes has been questioned, given its omnipresence in microbes normally not associated with plants (Gray et al., 1996), more recent research has evidenced its importance. Indeed, tRNA-derived cytokinins have been shown to contribute positively to nodule development in Aeschynomene plants upon infection with Bradyrhizobium (Podlešáková et al., 2013) and to play a pivotal role in virulence of diverse fungi, including Ustilago maydis (Morrison et al., 2017), Magnaporthe oryzae (Chanclud et al., 2016), and Claviceps purpurea (Hinsch et al., 2016). In Agrobacterium tumefaciens, cytokinins, albeit not tRNA-derived, have been implicated in efficient attachment to plant tissues contributing to transformation efficiency (Sardesai et al., 2013), whereas in symbiosis, they are believed to have unexpectedly diverse and context-dependent functions, including mediating root susceptibility to rhizobial infection (Frugier et al., 2008). So, the different cytokinin levels produced by both pathways in D188 and probably also in both PBTS isolates (Figures 2-4 and Supplementary Figure S2), might have many more functionalities than the mere induction of shoots. As the typical symptoms associated with an increased cytokinin response have been attributed to the degradation of auxin rather than the production of cytokinins in some leafy gall inducers (Lacey, 1948;Kemp, 1978), the chromosomes of PBTS1, PBTS2, and D188 were examined for the presence of homologs of the IAA catabolism (iac) genes implicated in auxin degradation in diverse bacteria (Leveau and Gerards, 2008), but none were detected. Despite this finding, we assessed whether the three strains could grow on IAA as sole carbon source. In agreement with the in silico prediction, PBTS1 and D188 did not proliferate on IAA, but, surprisingly, PBTS2 could ( Table 4). Although the underlying catabolic pathway remains unknown, auxin degradation might contribute to symptom development in the PBTS disease.
It was shown that D188 possibly modulates ethylene levels in its host by degrading the ethylene precursor 1aminocyclopropane-1-carboxylic acid (ACC) via ACC deaminase activity (Francis et al., 2016). In the PBTS2 chromosome, but not that of PBTS1, an acdS gene flanked by its regulator acdR were identified (Supplementary Table S12), which is an accordance with the ability of PBTS2, but not PBTS1, to FIGURE 7 | Effect of PBTS1 and PBTS2 on plant development. (A) Triple response assay used to evaluate ethylene production. PBTS1, PBTS2 and D188, as a reference, were inoculated on PDA medium in I-split plates and grown for 2 days at 28 • C, then 10 Arabidopsis Col-0 seeds were placed in the adjacent sector on 1 /2MS medium with 2% sucrose. PDA medium without any additions was used as negative control and 60 µM ethephon was either mixed in the 1 /2MS medium or spotted on the PDA medium as positive controls. The length of the seedling hypocotyls (indicated below the images) was measured after 7 days and statistical differences compared to the untreated control were calculated with Student's t-tests ( * significant difference at P < 0.05; n = 5-10). (B) Secondary differentiation of the vasculature of petioles of deformed leaves of "UCB-1" plants exhibiting PBTS symptoms compared to non-symptomatic leaves, suggesting an increased auxin response. (C) Production of plant growth promoting volatiles assessed in split plates on Arabidopsis Col-0 (right top quadrant) and N. tabacum W38 (left lower quadrant) grown on 1 /2MS medium without sucrose. The bacteria were grown on PDA. After 11 days of growth the shoots were excised to determine their fresh weight (indicated below the images) and the root system was photographed. Statistical differences in the average fresh weight compared to the control were calculated with Student's t-tests ( * significant difference at P < 0.01; n = 6 for N. tabacum (Nt); n = 5 for Arabidopsis (At)). use ACC as a nitrogen source (Table 4). To further explore a possible role of ethylene in the interaction between the PBTS bacteria and plants, we carried out an Arabidopsis triple response assay in split plates (Guzmán and Ecker, 1990). When ethephon, used as the positive control, was mixed in the MS medium in which the Arabidopsis seeds germinated, a full triple response was obtained, i.e., shorter and thicker hypocotyls with an exaggerated apical hook ( Figure 7A). However, when ethephon was applied on the PDA medium in the adjacent plate section, only an intermediate response was triggered, as illustrated by the significantly shorter hypocotyl length compared to the untreated control. Interestingly, both PBTS isolates provoked a similar intermediate response hinting at the production of low levels of ethylene, whereas D188 did not ( Figure 7A). Bacterial ethylene production has been reported in several genera and it is generally accomplished by one of two pathways. Ethylene can either be spontaneously produced at trace amounts via oxidation of 2-keto-4-methylthiobutyric acid or it is generated from α-ketoglutarate and arginine by the ethylene-forming enzyme (EFE) (Eckert et al., 2014). Importantly, the non-enzymatic production of ethylene by Pseudomonas solanacearum is eliminated in the dark (Swanson et al., 1979). Assuming that this is also the case in the PBTS isolates, this would suggest that the observed ethylene production is enzyme-mediated. Nevertheless, because no significant hits were found with different EFE sequences from other microbes nor with ACC oxidase 1 enzyme (ACO1) sequences of plants, the pathway responsible for the putative ethylene production in the PBTS strains remains unresolved. In the symbiotic interaction between the legume Sesbania rostrata and Azorhizobium caulinodans, which is mediated via crack entry rather than via root hair colonization, plantderived ethylene is required for rhizobial colonization. The combination of the bacterial Nod Factors and ethylene induces programmed cell death that allows the intercellular progression of the bacteria through the cortex toward the developing nodule primordium (D'Haeze et al., 2003;Guinel, 2015). Additionally, in the relatively recently evolved symbiosis between the non-legume Parasponia andersonii and rhizobia such as Bradyrhizobium elkanii, ethylene is also required for endophytic infection (Seifi Kalhor et al., 2016). Even more, in Pseudomonas syringae pv. glycinea, efe mutants were significantly impaired in their ability to grow in planta (Weingart et al., 2001). Thus, we speculate that the putative bacterium-derived ethylene plays a role in the endophytic colonization process. Alternatively, ethylene production by soil microorganisms has been shown to affect the viability of other microbes (Smith and Cook, 1974), so a role in niche protection could also be envisioned.
ACC deaminase activity is often detected in bacteria that produce auxins (Glick, 2014). Strain D188 produces auxin via a chromosome-encoded indole-pyruvate pathway in which tryptophan serves as a precursor (Vandeputte et al., 2005), and accordingly, a putative indole-pyruvate decarboxylase ipdC gene has been identified in its chromosome (A3L23_RS23370) (Francis et al., 2016). The much more efficient IAA production when tryptophol (indole-3-ethanol) is fed as a precursor to D188, signifies that the IpdC activity represents the rate-limiting step in auxin production (Vandeputte et al., 2005). In D188, auxin is believed to play a crucial role in the increased auxin signaling that is at the basis of the neovascularization of the developing leafy gall (Dolzblasz et al., 2018). Interestingly, light microscopy analysis of petioles of control and symptomatic "UCB-1" leaves, clearly showed thickening of the tissue and aberrations in the vasculature (Figure 7B), that are highly reminiscent of those observed upon infection of Arabidopsis with D188 (Dolzblasz et al., 2018), implying the putative involvement of auxin in the PBTS symptomatology. In this respect, in the PBTS1 and PBTS2 chromosomes, all the genes of a tryptophan biosynthesis pathway similar to that of D188 were identified as well as a gene homologous to A3L23_RS23370 (Supplementary Table S12 and Supplementary Figure S5).
Despite the weak homology of A3L23_RS23370 to the IpdC of Paenibacillus polymyxa (31% identity, 47% similarity) (Phi et al., 2008) and its ability to increase IAA production in Escherichia coli (Francis et al., 2016), the protein is annotated as acetolactate synthase. In the genome of P. polymyxa, two homologs of this gene are present: PPSC2_RS39450, which is the IpdC, and PPSC2_RS35795, annotated as a biosynthetic type of acetolactate synthase (large subunit); both proteins are 30% identical and 48% similar. Downstream of the latter, PPSC2_RS35800 encodes an acetolactate synthase (small subunit) followed by two genes, PPSC2_RS35810 and PPSC2_RS35815, involved in leucine biosynthesis (Ma et al., 2011). In the D188 and PBTS chromosomes, only one acetolactate synthase gene is present and the encoded protein has a slightly higher homology to PPSC2_RS35795 (50% identity, 66% similarity) than to the IpdC. Additionally, in all three chromosomes, an acetolactate synthase (small subunit) gene and two genes involved in branched-chain amino acid biosynthesis occur immediately downstream of this gene (Supplementary Table S12), an organization that is very similar to that of the second gene cluster in P. polymyxa. Thus, the ipdClike gene in the Rhodococcus chromosomes might function in branched-chain amino acid biosynthesis rather than in auxin biosynthesis (Supplementary Figure S5). Nevertheless, because the only copy of this gene present in the PBTS bacteria and in D188 has a low homology level to both of the proteins in P. polymyxa and since no other genes putatively related to auxin biosynthesis appear to be present in the The data are the average of three independent experiments.
Frontiers in Microbiology | www.frontiersin.org Rhodococcus chromosomes, we speculate that this IpdC-like protein might function in the two pathways by utilizing both pyruvate and indole-pyruvate as substrates. Recently, several indole-3-acetaldehyde dehydrogenases have been identified in Pseudomonas syringae that convert indole-3-acetaldehyde into IAA (McClerklin et al., 2018). In the three Rhodococcus strains, homologs are found to PSPTO_0092, the enzyme with the highest activity in P. syringae (Supplementary Table S12), that possibly encode the last step in IAA production via the indole-pyruvic acid pathway (Supplementary Figure S5). To validate the prediction that PBTS1 and PBTS2 can produce auxin, the bacteria were grown in the presence of tryptophan and tryptophol, the two auxin precursors that are used by strain D188 with different efficiencies (Vandeputte et al., 2005), and the occurrence of indolic compounds was evaluated using Salkowski's reagent (Glickman and Dessaux, 1995). When PBTS1, PBTS2, and D188 were grown in YEB medium supplemented with either of the precursors, no auxin could be detected with the colorimetric assay in any of the supernatants ( Table 5). In contrast, detectable auxin production by the three Rhodococcus strains could be demonstrated when tryptophan, but especially tryptophol, was provided, in line with previous findings on auxin production by D188 (Vandeputte et al., 2005). In the presence of tryptophol, PBTS2 and D188 produced a comparable level of auxin, which was about 5 times higher than that measured for PBTS1 under the experimental conditions used (Table 5). Altogether, in analogy with strain D188, these data suggest that PBTS1 and PBTS2 might produce auxin via the indole-pyruvate pathway with the IpdC activity representing the rate limiting step. Besides its role as an intermediate in the biosynthesis of branched-chain amino acids, acetolactate can spontaneously convert to diacetyl that can be used by an acetoin reductase to form acetoin (3-hydroxy-2-butanone), that is further transformed into 2,3-butanediol by a 2,3-butanediol dehydrogenase, a common bacterial pathway known as butanediol fermentation (Xiao and Xu, 2007). Interestingly, the latter two compounds have been identified as important volatiles produced by beneficial bacteria that exhibit significant plant growth-promoting effects and trigger induced systemic resistance (Ryu et al., 2005). Because of the presence of the two genes required for butanediol fermentation in the three Rhodococcus chromosomes (Supplementary Table S12 and Supplementary Figure S5), we evaluated whether the production of plant growth-promoting volatiles could be demonstrated in a plate assay with both Arabidopsis and N. tabacum as test plants. A strong biostimulating effect on the shoots as well as the roots of both plants was observed upon exposure to the volatiles of the three bacteria ( Figure 7C). The PBTS bacteria have been postulated to be beneficial rather than pathogenic (Savory et al., 2017). However, in contrast to the effect of the plasmid-free D188 derivative (Francis et al., 2016), no positive impact on plant development was ever observed under our experimental conditions when the PBTS bacteria were in contact with their host (Figure 2). Furthermore, the clear plant growth-promoting effect accomplished by the volatiles of the PBTS isolates, should not be considered as a proof of their alleged beneficial character because many notorious phytopathogens have been shown to produce volatiles that positively affect plants (Sánchez-López et al., 2016). Eventually, a function in niche occupation can be envisioned, since 2,3-butanediol has recently been shown to be important for rhizosphere robustness in the rhizobacterium Bacillus subtilis (Yi et al., 2016).

CONCLUSION
By combining previous results (Stamler et al., 2015b) with the data presented here, we conclude that subpopulations of the PBTS bacteria carry a fas locus that is highly similar to that present in leafy gall inducers. Additionally, for PBTS2, the fas operon is probably present on a linear plasmid, but for PBTS1, not enough data are currently available to support the same hypothesis. In contrast to many leafy gall inducers, the elusive virulence carrier in the two PBTS strains is highly unstable outside the plant host, affecting the outcome of pathogenicity assays. Nevertheless, the chromosomes of the PBTS strains encode multiple functions that could contribute to the interaction with and developmental modification of the plant host, including those implicated in the modulation of auxin and ethylene levels and in the production of cytokinins and growth-affecting volatiles. These findings suggest that, similar to leafy gall inducers, co-option between the virulence carrier and the chromosome is likely an important factor in the virulence strategy of these bacteria.
Our data also revealed that the co-occurrence of PBTS1 and PBTS2 is supported by a diversification in nutrient utilization patterns and in complementary activities that could stimulate niche partitioning. PBTS1 would strongly benefit from the presence of PBTS2 to occupy habitats that could otherwise be too hostile for this strain. In this context, it is tempting to speculate on the origin of PBTS: by associating with the pathogenic PBTS2, PBTS1, as an initially harmless inhabitant of pistachio, found itself in a situation permissive for virulence gene exchange which resulted in synergistic activities that strongly impacted the development of their host. Until the virulence carriers of PBTS1 and PBTS2 are identified and thoroughly analyzed, this scenario cannot be validated, but the chromosome analysis presented here allows to get a first glance at the mode of action of these recently emerged pathogens.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
DV, IF, RS, JK, and JR designed the experimental work. DV, IF, PL, JV, RS, JK, and JR performed the experiments and analyzed the data. YZ performed the bioinformatics analysis. DV wrote the paper. All authors critically read the manuscript and approved its final version.

FUNDING
Part of this research was supported by the California Pistachio Research Board and Redox Chemicals LLC. IF was awarded a Spring 2019 Research Development Grant program from the CSU Program for Education and Research in Biotechnology (CSUPERB). The funder, Board and Redox Chemicals LLC was not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication.