Monitoring of Rice Transcriptional Responses to Contrasted Colonizing Patterns of Phytobeneficial Burkholderia s.l. Reveals a Temporal Shift in JA Systemic Response

In the context of plant–pathogen and plant–mutualist interactions, the underlying molecular bases associated with host colonization have been extensively studied. However, it is not the case for non-mutualistic beneficial interactions or associative symbiosis with plants. Particularly, little is known about the transcriptional regulations associated with the immune tolerance of plants towards beneficial microbes. In this context, the study of the Burkholderia rice model is very promising to describe the molecular mechanisms involved in associative symbiosis. Indeed, several species of the Burkholderia sensu lato (s.l.) genus can colonize rice tissues and have beneficial effects; particularly, two species have been thoroughly studied: Burkholderia vietnamiensis and Paraburkholderia kururiensis. This study aims to compare the interaction of these species with rice and especially to identify common or specific plant responses. Therefore, we analyzed root colonization of the rice cultivar Nipponbare using DsRed-tagged bacterial strains and produced the transcriptomes of both roots and leaves 7 days after root inoculation. This led us to the identification of a co-expression jasmonic acid (JA)-related network exhibiting opposite regulation in response to the two strains in the leaves of inoculated plants. We then monitored by quantitative polymerase chain reaction (qPCR) the expression of JA-related genes during time course colonization by each strain. Our results reveal a temporal shift in this JA systemic response, which can be related to different colonization strategies of both strains.


INTRODUCTION
Plant microbiome is nowadays extensively studied, as it represents a huge potential for agriculture. Numerous studies describe the importance of microbes for plant's nutrient supply and resistance to diseases and pests (Finkel et al., 2017). Therefore, microbes are potential solutions to do the transition to a sustainable agriculture while maintaining yield (Busby et al., 2017). Especially, some rhizobacteria have been shown to have tremendous beneficial effects on plant growth (Hayat et al., 2010) and resistance to pathogens (Beneduzi et al., 2012). These beneficial effects are induced through hormonal modulations following the colonization of plant roots (Vacheron et al., 2013) and inner tissues for endophytes (Hardoim et al., 2008) as well as systemic regulations of immunity (Pieterse et al., 2014). However, the perception of microbe-associated molecular patterns (MAMPs) generally leads to an immune response called MAMPs-triggered immunity (MTI) characterized by the synthesis of antimicrobial compounds (Jones and Dangl, 2006). In the same way pathogens suppress plant immunity, and beneficial microbes are able to escape (Trdá et al., 2015) or modulate (Zamioudis and Pieterse, 2012) the immune response of plants cells. Interestingly, the suppression of MTI by both pathogenic bacteria (Millet et al., 2010) and beneficial fungi (Jacobs et al., 2011) is commonly mediated via the jasmonic acid (JA) signaling pathway. However, the global physiological response and especially the transcriptional regulations induced by plants during the interaction with beneficial bacteria are not well described.
Within the diversity of bacteria interacting with plants, the Burkholderia sensu lato (s.l.) genus of Betaproteobacteria stands out for several reasons. It contains plant pathogenic species (B glumae, Burkholderia gladioli, and Burkholderia plantarii) (Maeda et al., 2006), N 2 -fixing nodule-forming species in association with tropical legumes (Gyaneshwar et al., 2011) as well as species-forming associative symbiosis particularly with cereals such as rice (Coutinho et al., 2013;Govindarajan et al., 2008). Also, a phylogenetic separation discriminates the Burkholderia sensu stricto (s.s.) genus, which contains animal or plants pathogens as well as human opportunists, and the Paraburkholderia genus, which contains mainly plant-associated and environmental species (Sawana et al., 2014;Estrada-de los Santos et al., 2016). Also, most recent phylogenetic refinements of Burkholderia s.l. taxa defined other genera-Caballeronia, Robbsia, Trinickia, and Mycetohabitans-which are supported by both differences in genomic and ecological features (Estrada-de los Santos et al., 2018).
Interestingly, two strains of the Burkholderia s.l. genus able to fix N 2 have been described as growth promoters in rice. Paraburkholderia kururiensis M130 (hereafter Pk) is a beneficial rice endophyte (Mattos et al., 2008) related to environmental and plant-beneficial strains (Kaur et al., 2017). Burkholderia vietnamiensis TVV75 T (hereafter Bv) is a rice-associated species having positive effect on yield (Trân Van et al., 2000;Govindarajan et al., 2008), which belongs to the Burkholderia cepacia complex, a complex of species that can cause serious risks to cystic fibrosis patients (Vial et al., 2011). In order to decipher how rice perceives these two beneficial strains belonging to genera with contrasted ecologic backgrounds, we studied the transcriptional responses of rice during the establishment of the interaction. Our aim was to identify plant physiological processes and potentially key genes involved in beneficial rice-rhizobacteria interactions and also differentially regulated by each strain. We first analyzed the colonization patterns of Pk and Bv on the Oryza sativa Nipponbare genotype. We then analyzed the root and leaf transcriptional responses to the bacterial colonization by RNAseq. This led us to the identification of a co-expression JA-related network. Therefore, we monitored, throughout the establishment of the interaction, the expression of JA-related genes by reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR).

Plants and Bacterial Cultivation
O. sativa L. ssp. japonica cv Nipponbare was used in this study. For all experiments, seeds were dehusked and sterilized as follows: 70% ethanol for 10 min and 9.6% NaClO supplemented with 1% Tween 20 for 30 min. Treated seeds were rinsed twice with sterile distilled water, twice with 2% thiosulfate solution, and finally four times with sterile distilled water. To confirm surface sterilization, 100 µl of the last rinsing solution was plated in tryptic soy agar (TSA) medium (Sigma-Aldrich). Seeds were then put in sterile distilled water at 28°C for 24 h and transferred on 8% H 2 O agar plate for 30 h. Homogeneously germinated seeds were transferred to sterile magenta boxes (SPL Lifesciences Co. Ltd) containing 150 ml of autoclaved perlite and 200 ml of sterile hydroponic medium (recipe in Supplementary Table 1). Plants were grown in a growth chamber (16 h light; 8 h dark; 28°C; 70% humidity).
All bacterial strains (listed in Supplementary Table 2) were cultured as follows: Glycerol stocks (20%) of bacterial cells conserved at −80°C were plated in low-salt lysogeny broth (LB) (Sigma-Aldrich) agar plates and incubated for 72 h at 28°C. Liquid low salt LB medium of 20 ml was then inoculated in 50-ml Falcon tubes and incubated for 16 h under agitation (180 rpm) at 28°C. Overnight culture of 500 µl was inoculated in fresh liquid medium for 2 h. Bacterial cells were then centrifuged for 5 min at 4,000 rpm and resuspended in sterile distilled water. Each plantlet was inoculated with 10 7 bacterial cells 4 days after sowing in hydroponic system.

Bacterial Transformation
P. kururiensis M130 and B. vietnamiensis TVV75 cells were transformed by electroporation with the pIN29 plasmid (Vergunst et al., 2010). The plasmid chosen to transform the strains, pIN29, comprises a chloramphenicol resistance gene as well as the DsRed gene under the control of a constitutive TAC promoter. After 24 h of incubation of selective medium low salt LB Cm (200 μg·ml −1 ) at 28°C, the most fluorescent colonies were selected.

Rice Root Colonization Assays
The roots of plants were harvested at 1, 7, and 14 days postinoculation (dpi), weighted, and grinded in sterile water with a sterile ceramic bead using a FastPrep-24 ™ 5G at 6 m·s −1 for 40 s. The solution was then diluted and inoculated in low-salt LB selective medium containing 200 μg·ml −1 of chloramphenicol and incubated at 28°C for 24 h. Colony-forming units were then enumerated. The size of the root-associated bacterial population was measured during two independent experiments, each comprising nine plants. In order to measure the size of the endophytic population, the inoculated rice roots were surface disinfected for 1 min using a solution of 1% chloramine T (Sigma-Aldrich) supplemented with 0.1% Tween 20. Roots were then rinsed six times with sterile water. Controls of disinfection were performed by plating rinsing water in TSA medium (Sigma-Aldrich) overnight. Surface-disinfected roots were then treated as previously described. The size of the endophytic population was measured on five plants.

Microscopy
All microscopic observations of the bacterial colonization were restricted to the primary root in order to compare the colonization patterns on roots that have been in contact with the bacterial population for the same amount of time. Primary roots were harvested at 7 and 14 dpi and mounted between slide and slips cover and directly examined with the microscopes. Epifluorescence observations were performed using a Nikon Eclipse Ni-E microscope. Confocal Laser Scanning observations were performed using a Zeiss LSM880 confocal microscope.

RNA Extraction
For the analysis of root and leaf transcriptional profiles, both plants' tissues were harvested at 6 h postinoculation (hpi), 1 dpi, 7 dpi, or 14 dpi with live bacterial cells. Each biological replicate consisted of five pooled root system or five pooled last mature leaves harvested from a single hydroponic system. For each time point and each inoculated strain, three biological replicates were harvested. Roots and leaves of untreated plants were collected at the same time points. After harvest, samples were snap-frozen in liquid nitrogen and stored at −80°C.
Rice roots were homogenized in liquid nitrogen using cooled mortar and pestle. Rice leaves were grinded using a TissueLyser II (Retsch) set to 30 Hz for 30 s. Total RNA extraction using TRIreagent (Sigma) was performed according to manufacturer's instructions. All samples were treated with DNase I (Ambion) and purified using the RNA Clean & Concentrator kit (Zymo) according to manufacturer's instructions. The integrity and quality of the total RNA were confirmed using a NanoDrop ™ 1000 spectrophotometer (Thermo Fisher) and a 2100 BioAnalyzer (Agilent).

RNA Sequencing and Mapping of Reads
Quality of RNA was checked by determining the RNA Integrity Number (RIN) with a Fragment Analyzer (Agilent). For the library preparation, samples with a RIN value > 6 were used. Eighteen RNA libraries were prepared using an Illumina TruSeq stranded mRNA sample preparation kit by MGX-Montpellier GenomiX core facility (MGX) France (https://www.mgx.cnrs.fr/). Library construction and sequencing were performed as described in Karmakar et al. (2019) on an Illumina HiSeq 2500. The quantitative and qualitative validation of the library was performed by qPCR, Roche LightCycler 480, and a Fragment Analyzer (Agilent) using a Standard Sensitivity NGS kit. Quality control and assessment of raw Illumina reads in FASTQ format were done by FastQC software (version 0.11.5) to obtain per base quality, Guanine-Cytosine (GC) content, and sequence length distribution. Clean reads were obtained by removing the lowquality reads, adapters, and poly-N-containing reads by using Trimmomatic v0.36 software (Bolger et al., 2014). RNAseq reads were aligned to the IRGSP 1.0 version of the rice genome using HISAT2 v2.0.5.1 (Kim et al., 2015). The number of reads mapped to each gene locus was counted using HTSEq-count v0.6.0 (Anders et al., 2015).

Differential Gene Expression and Gene Ontology (GO) Term Enrichment Analysis
DESeq2 v3.7 (Love et al., 2014) was used to calculate differential gene expression between non-inoculated and inoculated conditions. All genes having an adjusted p-value inferior to 0.01 were considered as significantly differentially expressed. All functional enrichment analyses were performed using g:Profiler (version e95_eg42_p13_f6e58b9) with g:SCS multiple testing correction method applying significance threshold of 0.05 (Reimand et al., 2007;Raudvere et al., 2019).
Quantification of mRNA Levels Using RT-qCPR cDNA was produced from 350 ng of DNase-treated total RNA using the SuperScript III Reverse Transcriptase (Thermo Fisher Scientific). The cDNA reaction was diluted five times before qPCR using MESA BLUE qPCR Master Mix for SYBR ® assay (Eurogentec) on an Mx3005P qPCR system (Agilent Technologies). The relative expression level was calculated according to Pfaffl (2001). Three independent samples were analyzed for each condition, and each sample was assayed in triplicate. Primers used are listed in Supplementary Table 3.

Analysis of Root Colonization
We used hydroponic culture of rice plants, grown in axenic condition and inoculated with DsRed-tagged strains, to monitor the bacterial colonization at 1, 7, and 14 dpi. First, the root's colonization was measured by counting the bacterial populations on the rhizoplan and endosphere (see Materials and Methods). The roots of rice plants were rapidly colonized by both bacterial strains ( Figure 1A). The populations of Pk and Bv reach a median value of 4.16 × 10 6 and 1.73 × 10 6 cfu·g −1 , respectively, at 1 dpi. At 7 dpi, the maximum measured size of the bacterial population is reached for both strains: 7.49 × 10 8 and 3.18 × 10 8 cfu·g −1 for Pk and Bv, respectively. Then, between 7 and 14 dpi, the size of the total root-associated population decreases for each strain, reaching a median value of 4.6 × 10 8 and 1.46 × 10 8 cfu·g −1 for Pk and Bv, respectively. Between the same time points, the variation of the endophytic population size differs between the two strains. Indeed, the endophytic population size of Pk decreases from 1.25 × 10 6 to 1.17 × 10 4 cfu·g −1 between 7 and 14 dpi, while the median number of endophytic Bv cells remain stable with 7.34 × 10 4 and 7.5 × 10 4 cfu·g −1 at 7 and 14 dpi, respectively.
Microscopic observations of the primary roots of inoculated rice plants demonstrated that both bacterial strains tagged with the DsRed gene colonized the root surface after inoculation (Figures 1B, C). Moreover specific zones were more densely colonized by bacteria such as the surface of root hairs and the emergence of lateral roots (Supplementary Figure 1). By comparing the colonization of the two strains, differences can be observed in the way they colonized the surface of the primary root. Several epidermal plant cells seemed colonized intracellularly by the tagged Bv cells, while the phenomenon was observed less frequently in Pk-colonized roots ( Figure 1C; Supplementary  Figure 1). Observations by confocal microscopy revealed that for Pk, the highly colonized epidermal cells appear to be only colonized on their surface as the whole outline of epidermal cells as well as in intercellular spaces (Figures 2A, C), while for Bv, most of the observations show that bacterial cells were able to cross the cell wall and were observed in the cytoplasm of the cell (Figures 2B, D). Thus, both strains colonized the roots of the Nipponbare cultivar both externally and endophytically but through apparently different intensities and entry roads. We then wondered if the host plant induces contrasted transcriptional regulations associated with the differential colonization pattern of each strain.

Transcriptional Response of Rice to Bacterial Inoculation
In order to identify changes in O. sativa transcriptome in response to both strains, we performed RNAseq on leaves and roots of non-inoculated controls, Pk-inoculated and Bv-inoculated plants at 7 dpi. We chose this time point to avoid the initial plant defense burst due to a bacterial inoculation in hydroponic system (hours to 1 dpi) and allow an advanced colonization stage of the roots but without any visible developmental effect such as increased growth. To confirm that the inoculated and noninoculated control plants were at the same developmental stage, we measured the dry weight of the plants and could not detect any significant impact of the inoculation on biomass production (Supplementary Figure 2) nor on plant height. Therefore, we assume that the differences between the transcriptomes of inoculated plants and non-inoculated controls should be related to bacterial colonization of the roots rather than a developmental impact of inoculation.
A total of 843 million reads were sequenced with an average of 47 million reads per sample ( Table 1). An average of 68% of reads was uniquely mapped to the O. sativa genome per sample. A principal component analysis also discriminates the transcriptome of non-inoculated plants compared with the inoculated ones as well as the plant responses to each bacterial strain (Figures 3A, B). Differential expression analysis yields a total of 4,951 and 5,275 significantly (p < 0.01) differentially expressed genes (DEGs) in response to Bv and Pk, respectively. Comparing leaf and root transcriptomes reveals that there are five and eight times more DEGs detected in leaves in response to Bv and Pk, respectively, than in roots (Figures 3C, D). When comparing the response to Pk and Bv, a large proportion of DEGs are commonly regulated in leaves (around 50% in response to both strains) contrarily to roots in which the transcriptional regulation appears to be more specific to each inoculated strain. Indeed, only 33% of the DEGs in response to Bv are also differentially expressed in response to Pk. The proportions of commonly up-regulated DEGs in roots represent 20% and 27% of DEGs in response to Bv and Pk, respectively.  Table 4). Three main biological processes are commonly regulated during the interaction with both strains in leaves and roots: response to stimuli, and metabolic and also developmental processes. First, stress-related genes are enriched in the commonly DEGs. Indeed, in leaves, the "response to abiotic stimulus" as well as the "response to oxidative stress" terms are enriched in the up-regulated DEGs, whereas defense-related GO terms are enriched in the down-regulated DEGs in both leaves and roots. Also, several hormone-related GO terms are enriched in the DEGs in both leaves and roots. In leaves, GO terms related to auxins and abscisic acid (ABA) response are enriched in the up-regulated genes, while in roots, the up-regulated DEGs are enriched in cytokinines (CK), brassinosteroids (BR), and ethylene response terms. Finally, in roots, down-regulated DEGs are enriched in gibberellic acid (GA) and salicylic acid (SA) response-related terms. The analysis also revealed that metabolic processes are transcriptionally regulated in response to the interaction with both strains: In leaves, the up-regulated DEGs are enriched in the "photosynthesis" and "translation" terms, while the down-regulated DEGs are enriched in "starch metabolic process" and "membrane lipid catabolic process" terms. Furthermore, in leaves, the "iron ion homeostasis" term is enriched in up-regulated DEGs, while the "metal ion transport" term is enriched in the down-regulated DEGs of roots. Although the inoculation of both strains did not significantly impact the biomass of rice plants (Supplementary Figure 2), several development-related GO terms are enriched in the DEGs in both leaves and roots. First, in leaves, among others, the "anatomical structure development" term is enriched in the up-regulated DEGs, and the "glucan biosynthetic process" term is enriched in the down-regulated DEGs, which correspond to genes implicated in cellulose and callose synthesis. Finally, root transcriptional response is also enriched in development-related GO terms such as the "xylem development" term enriched in the up-regulated DEGs and the "lignin biosynthetic process" term enriched in the down-regulated DEGs.
Additionally, in order to identify the biological processes specifically induced during the interaction with each strain, we performed a GO term enrichment analysis on the DEGs specifically regulated by each strain (Supplementary Figures 4 and 5; Supplementary Tables 5 and 6). This analysis revealed that leaves of plants inoculated with Pk up-regulate genes related to biosynthetic process and translation, while the interaction with Bv induces the up-regulation of genes coding for components of the photosystem II and also the downregulation of protein folding genes (Supplementary Figure 4, genes listed in Supplementary Table 7). Interestingly, the up-regulated DEGs in response to each strain are enriched in one hormonal signaling pathway. Indeed, the interaction with Pk induced the up-regulation of JA-related genes while cytokininrelated genes are enriched in the up-regulated DEGs in the leaves of Bv-inoculated plants. The analysis of the specific root transcriptome also revealed processes specifically induced by each strain. Particularly, the response to Pk in roots encompasses the down-regulation of oxidative stress response-related genes and also chitin catabolic process, while the interaction with Bv induced the down-regulation of genes involved in defense, JA signaling, and response to stimuli (Supplementary Figure 5, corresponding genes in Supplementary Table 8).
In order to deepen the transcriptome analysis, a functional categorization of the top 200 up-regulated and down-regulated DEGs identified in response to each strain was carried out using all available databases for rice annotation (UniProt, Kyoto Encyclopedia of Genes and Genomes (KEGG), RAP-DB, and Oryzabase). The proportion of genes related to each category is presented in Figure 4 and exemplifies rice-specific responses to the two bacterial strains. In leaves, more stress-related genes as well as genes involved in secondary metabolism are up-regulated in response to Pk compared with the response to Bv. Also, in response to Pk, 16 genes related to chromatin remodeling are down-regulated compared with two in response to Bv. Inversely, the inoculation of Bv induced the down-regulation in leaves of 22 stress-related genes, whereas only three are down-regulated in response to Pk. In roots, twice as much transcription regulation and protein degradation-related genes are up-regulated in response to Bv than in response to Pk. Conversely, approximately twice as much "nutrition/transport" and signaling-related genes are down-regulated in response to Pk than in response to Bv. In order to identify plant key genes that could be related to the observed patterns of root colonization, we focused our analysis on defense (leaves in Table 2, roots in Table 3) and hormone-related genes (leaves in Table 4, roots in Table 5). Indeed, these processes are enriched among the strain-specific DEGs (Supplementary Figures 4 and 5) and could pinpoint differences in the physiological response of rice following the perception of each strain.
In roots, 13 defense-related genes are commonly regulated, while 33 and 18 are specifically Pk and Bv regulated (Table 3). Additionally, 12 hormone-related genes (implicated in JA, CK, ethylene, GA, and ABA syntheses or signaling) are commonly regulated in response to both strains, while 14 and 19 hormonerelated genes are Pk and Bv regulated, respectively ( Table 5). The latest encompasses genes implicated in auxins, BR, SA, and strigolactone synthesis, or signaling. In leaves, 41 defenserelated genes are regulated, six being commonly regulated, while 11 and 24 are specifically induced by Pk and Bv, respectively ( Table 2). Also hormone-related DEGs are detected in leaves; six genes implicated in GA, JA, auxins, ethylene and ABA synthesis, transport, or signaling are commonly regulated. Finally, in response to each strain, 20 hormone-related specific DEGs are detected. Interestingly, only the interaction with Pk induced the regulation of BR-and SL-related genes, while only Pv induced the regulation of a SA-related gene. Other hormonal pathways, such as CK, ABA, ethylene, JA, auxins, and GA, are modulated by each strain but induced the regulation of different genes ( Table 4).

Validation of RNAseq Data by qCPR
To validate the RNAseq data, we selected genes implicated in defense, hormone signaling, or development regulated at 7 dpi by one or both strains (Supplementary Table 9). We measured the level of expression of the selected genes by qRT-PCR in an independently conducted experiment. According to RNAseq analysis, in leaves, two defense-related genes are specifically regulated by each strain: ALD1 is specifically down-regulated in response to Bv, and WRKY71 is specifically up-regulated in response to Pk. Also, two hormone-related genes were detected as specifically regulated by each strain; namely, bHLH148, a JA signaling component, is down-regulated in roots inoculated with Bv; while RR9, a CK signaling component, is up-regulated in the leaves of Pk-inoculated plants. Finally, in roots, RSL9 is up-regulated and SHR5 is down-regulated in response to both strains. Gene expression changes obtained through qPCR analysis demonstrated a pattern similar to RNAseq (Figure 5) except for the over-expression of RR9 in response to Pk, which was not detected as significant by DESeq2.

Temporal Analysis of Strain-Specific Marker Genes
Following the confirmation of strain-specific transcriptional regulations, we wanted to identify genes that are differentially expressed in both conditions but with opposite regulations. Among all DEGs detected in roots and leaves, only three DEGs are potential differential markers: RERJ1 (Os04g0301500), a JA-responsive gene; ATL15 (Os01g0597600), a putative amino acid transporter; and DREB1B (Os09g0522000), a droughtresponsive gene. Noteworthy, these three genes are only detected as DEGs in leaves and follow the same pattern of expression, as they are up-regulated in response to Pk and down-regulated in response to Bv.
Among these three genes, two of them, ATL15 and RERJ1, are part of a co-expression network recovered from RiceFREND database , which contains four JA-related genes (JAZ6, 10, 12, and AOS1) ( Figure 6A). As JA is one of the main phytohormones implicated in defense (Pozo et al., 2004), we further wanted to know if the whole co-expression network could act as a differential marker of the response to each bacterial strain. In order to describe the transcriptional regulation of this co-expression network throughout the establishment of the  interaction, we produced a transcriptional kinetic of rice leaves tissues at 6 hpi, 1 dpi, 7 dpi, and 14 dpi and analyzed all six genes by qRT-PCR. First, at 7 dpi, we can confirm that depending on the inoculated strains, the regulation of the co-expression network is opposite (Figure 6; Supplementary Table 10). On the one hand, JAZ6, JAZ10, JAZ12, ATL15, and AOS1 are up-regulated following the inoculation with Pk; and on the other hand, JAZ6 and JAZ10 are down-regulated in response to Bv. Interestingly, and in sharp contrast to the 7 dpi response, for the short-term response (6 hpi), we observed the up-regulation of five out of the six genes of the co-expression network in rice leaves after inoculation with Bv, whereas only JAZ12 is up-regulated in response to Pk. Then, the whole network appears to decline for the rest of the kinetics in response to Bv. Eventually, at 14 dpi, the expression levels of the six genes are quite comparable between the two conditions being all down-regulated compared with non-inoculated controls.

DISCUSSION
The aim of this study was to describe the transcriptional regulations induced by rice following the perception of naturally associated beneficial bacteria. Additionally, we took advantage of the particular phylogenetic organization of the Burkholderia s.l. genus to compare the responses of plants with two closely related beneficial species with different phylogenetic backgrounds in terms of ecology.

P. kururiensis M130 and B. vietnamiensis TVV75 Differentially Colonize Roots of O. sativa Cultivar Nipponbare
We first observed that both bacteria were able to efficiently colonize the rice Nipponbare roots. The amount of culturable bacterial cells associated with rice roots appears to be coherent with the literature, as Compant et al. (2010) described that generally a range between 10 7 and 10 9 cfu·g −1 of root fresh weight is found colonizing roots both externally and internally. The same goes for the endophytic population, which generally ranges between 10 5 and 10 7 cfu·g −1 . The decrease of the bacterial population size between 7 and 14 dpi may be due to the increase of root biomass by the formation of newly emerged roots, which are not importantly colonized at least in the time of our experiments. However, by comparing the colonization of the two strains, two main differences can be observed in the dynamic of root colonization. First, Pk forms a significantly bigger population while colonizing rice roots surface than does Bv at every time postinoculation. Second, the dynamic of the endophytic population of both strains is very different, as the number of endophytic Pk cells declines between 7 and 14 dpi which is not the case for Bv. This could be due to the fact that the endophytic colonization by Pk is restricted by the plant throughout time in contrast to Bv, which maintains its population size. From this observation, two hypotheses can be proposed: Either the plant is not able to control the colonization by Bv, or this strain is more efficiently colonizing the newly emerged roots.
Both the maturation zone and the area of lateral root emergence were identified as hotspots for Pk and Bv colonization. A similar area of colonization was identified for Paraburkholderia phytofirmans strain PsJN in grape plants (Vitis vinifera L.) (Compant et al., 2005) as well as for B. vietnamiensis MGK3, which also intensively colonizes the same areas of the roots of rice plants (Govindarajan et al., 2008). Root exudates, which contain essential nutrients for microbes, are released in the lateral root emergence zones (Badri and Vivanco, 2009). This may aid colonization and allow the possible entry of bacteria via mechanisms such as "crack entry" into the internal tissues (Hardoim et al., 2008).
Both strains were observed massively colonizing the maturation zone, on both the outside and even the inside of some root hair cells (Supplementary Figure 1). The accumulation of bacterial cells in the root hair zone has already been described for Pk (Mattos et al., 2008), and it was proposed to be a common hotspot for rhizobacterial colonization (Compant et al., 2010), as it is correlated to a higher local exudate concentration (Gamalero et al., 2004). Apparent intracellular colonization of epidermal cells was observed for both species but more frequently for Bv. In the case of Bv, the  fact that the vacuole is still intact when observing through the colonized epidermic cell offers evidence that the cell is still living (Figure 2D). Also, as the bacterial cells seem to have passed through the cell wall, the fluorescence signal appears to be cytoplasmic. Intracellular colonization of rice root epidermal cells by bacteria was also observed during the interaction with Azoarcus BH72  as well as during the colonization of ryegrass roots by Paraburkholderia bryophila Ha185 (Hsu et al., 2018). Furthermore, intracellular colonization of root hair cells was also observed by Prieto et al. (2011). They showed by confocal microscopy that Pseudomonas putida PICP2 and Pseudomonas fluorescens PICF7 are able to colonize root hair cells of olive trees (Olea euroapea L.) and subsequently move into epidermal cells. This observation led them to propose a new route of entry in intern plant tissues for endophytic bacteria, which starts with the colonization of single root hairs.

Root Colonization Induces More Transcriptional Regulations in Aerial Parts Than in Roots
The analysis of root and leaf transcriptomes revealed that colonization of rice roots by both strains induced more transcriptional regulations in leaves than in roots ( Figure 3B). This could be due to the fact that the inoculated plants were harvested at 7 dpi and therefore at an established state of the interaction between bacteria and rice. In contrast, previous studies revealed that the inoculation of beneficial rhizobacteria induced the regulation of at least 1,000 genes in rice roots at earlier time points after inoculation (Drogue et al., 2014;Brusamarello-Santos et al., 2019;Rekha et al., 2018). Moreover, compared with the only study that analyzed the leaf transcriptional response of rice to the inoculation by beneficial rhizobacteria that retrieved only 2,414 DEGs at early stages of the interaction (Wu et al., 2018), when our study identified data of at least 4,000 DEGs in leaves in response to both strains.

Roots Trigger Contrasted Transcriptional Response Depending on the Inoculated Strain
In response to the colonization by both strains, transcriptional reprogramming of defense-related genes occurs in rice roots ( Table 3). Only 13 of them appeared commonly regulated in the same way. First, three WRKY (62, 76, and 86),genes are commonly down-regulated; interestingly, those three genes are negative regulator of defense (Peng et al., 2008;Yokotani et al., 2013)-putatively for WRKY86 (Choi et al., 2017). Also, two negative regulators of NH1, the rice NPR1 homolog of the major SA response regulator (Yuan et al., 2007;Chern et al., 2012), namely, RH2 and RH3, are also down-regulated (Chern et al., 2014). Taken together, as these genes are described as negative regulators of defense, their down-regulation may reflect a common defense response towards root colonizing bacteria. The remaining 51 genes were differentially regulated in response to each strain (33 responding to Pk and 18 to Bv, respectively). Indeed, roots colonized with Pk induce the down-regulation of 11 pathogenesis-related (PR) genes as well as five chitinases and two xylanase inhibitors. One of the down-regulated PR genes is PBZ1 (Os12g0555200), and it is described as a defense marker in leaves, but it is also up-regulated in response to root invasion by Magnaporthe oryzae (Marcel et al., 2010). Moreover, bZIP79, which is described as a phytoalexin synthesis suppressor (Miyamoto et al., 2015), is up-regulated. On the other hand, roots colonization by Bv induce the down-regulation of only one chitinase, whereas two PR genes are up-regulated, of which one gene, BetvI, is targeted by parasitic nematode to suppress root defense (Chen et al., 2018). Furthermore, four WRKY genes of which two are thought to encode negative regulators (WRKY28  and WRKY79 (Choi et al., 2017)) and WRKY40, which is up-regulated in striga-resistant rice roots (Swarbrick et al., 2008), are here down-regulated. Taken together, these results show that the colonization of rice roots by Pk is associated with a downregulation of PR genes, while the colonization by Bv is associated with the down-regulation of defense-suppressing WRKY genes. This transcriptional regulation of defense-related genes could be the consequence of the more invasive colonization pattern of Bv cells described above.
As previously described (Vacheron et al., 2013;Drogue et al., 2014;Rekha et al., 2018), the inoculation of beneficial rhizobacteria induces important hormone-related transcriptional regulation in rice roots (Table 5). Similarly, the inoculation of Pk and Bv induced the regulation of genes encoding for signaling components of CK and ethylene. Particularly, both bacterial strains induced the up-regulation of two genes encoding for ACO, which are enzymes  implicated in the synthesis of ethylene (Ravanbakhsh et al., 2018). Interestingly, in response to the inoculation of other beneficial bacteria, genes coding for ACOs were down-regulated in rice roots (Drogue et al., 2014;Brusamarello-Santos et al., 2019). Also, important transcriptional regulations of ethylene responsive factors (ERFs) occur in response to both strains. In the same way, cytokinine signaling is impacted as both conditions induce the up-regulation of at least two RR genes and one cytokinine oxidase encoding gene (Tsai et al., 2012); however, it is not the same genes that are regulated by each bacterium.

Systemic Responses Associated With Root Inoculation Suggest a Time Shift in Defense Response Between the Two Species
Interestingly, the RNAseq analysis revealed that the inoculation of rhizobacteria induced major transcriptional regulations in leaves (Tables 2 and 4; Supplementary Table 5). As previously described (Persello-Cartieaux et al., 2003;Verhagen et al., 2004;Campos-Soriano et al., 2012), the colonization of beneficial microbes induces transcriptional regulations of defense-related genes in the leaves of the inoculated plants.
Both strains induced the up-regulation of WRKY7 as well as the down-regulation of three putative R genes and two PR genes. Pk induced the up-regulation of WRKY71, which is known to confer enhanced disease resistance to Xoo (Liu et al., 2007) and the down-regulation of TFX1, which is known as a susceptibility gene for Xoo (Sugio et al., 2007). On the other hand, Bv induces the up-regulation of five R genes: three (NB-ARC, MLO4, and STV11) (Wang et al., 2014), NPR4, as well as the down-regulation of three WRKY genes (9, 25, and 69), of which two are SA responsive (Liu et al., 2007;Choi et al., 2017), and also the down-regulation of ALD1, a basal immunity regulator needed for the accumulation of SA (Jung et al., 2016). . For RNAseq, the data presented are the mean count of normalized reads. For qRT-PCR, transcript levels were normalized to that of the reference gene EF 1α (Os03g0177400). Stars represent significant differences of mean expression compared with the control condition according to post-hoc test on a generalized linear model. * corresponds to p < 0.05 and ** to p < 0.0001. Error bars represent the standard deviation (n = 3). Bv, Burkholderia vietnamiensis TVV75 T ; Pk, Paraburkholderia kururiensis M130; DEGs, differentially expressed genes.
It is also in this organ only that differential markers have been detected, some of which designated the JA signaling pathway as a putative marker of the interaction with the two bacterial species. Also, interestingly, in contrast with the regulation of JA-related genes, several GA-related genes seem to be down-regulated in response to Pk and up-regulated in response to Bv. This antagonism between JA and GA has been described (Yang et al., 2012;Yimer et al., 2018) and proposed as one of the ways for plants to fine-tune the balance between growth and defense (Huot et al., 2014). Nonetheless, this apparent contrast in terms of JA-related genes transcriptional regulation at 7 dpi may only be a consequence of a delayed JA systemic signaling in response to Pk compared with the response to Bv, which induced the up-regulation of the co-expression network as soon as 6 hpi. As JA-induced plant defenses have been proposed to contribute to the restriction of endophytic colonization in grasses (Miché et al., 2006), this temporal shift in the induction of JA-related genes between the responses to the two bacterial strains could be associated with a delayed JA-induced defense signal. We propose the following interpretation: the perception of Bv induces very early (6 hpi) the up-regulation of JA-related genes in leaves to restrict the colonization, whereas in response to Pk, this JA signal happens at 7 dpi and results in the decrease of the bacterial population ( Figure 1A).
The comparison of the interactions between rice and Pk and between rice and Bv revealed important differences in the process of root colonization and rice transcriptional regulations induced by each strain, which we summarize in Figure 7. First, the numerous intracellular colonization of root epidermic cells by Bv resembles a pathogen infection compared with the apoplastic colonization observed for Pk and other beneficial endophytes (McCully, 2001;Reinhold-Hurek et al., 2007). However, such intracellular colonization by beneficial endophytes has already been observed (Reinhold-Hurek and Hurek, 1998). Also, the specific root response to Pk is characterized by the downregulation of gene coding for chitinase and PR proteins, while Bv colonization specifically induced the down-regulation of several WRKY genes (Table 3). Moreover, Bv specifically induced the down-regulation of JA signaling genes (Supplementary Figure 5; Supplementary Table 6) in addition to the common SA-related down-regulation in roots. Therefore, the strategies to circumvent the immune system of roots appear to be different between the two strains. Finally, the fact that the root inoculation of Bv induced the up-regulation of JA-related genes in leaves in only 6 hpi compared with the delayed similar signal in response to Pk colonization at 7 dpi also supports the fact that Bv induces a response similar to that of pathogens. Indeed, the rice root colonization by M. oryzae induced the up-regulation of JA-related genes in leaves at 3 and 4 dpi (Marcel et al., 2010). All these elements support the fact that Bv appears to have a much more invasive colonization strategy in terms of both patterns and modulation of the plant immune system. This statement is in accordance with the opportunistic and pathogenic background of the Burkholderia s.s. genus (Eberl and Vandamme, 2016). Consequently, it would be of interest to analyze the response of the cultivar Nipponbare to a larger diversity of plant-associated Burkholderia and Paraburkholderia species to investigate if the FIGURE 6 | JA-related co-expression network transcriptional regulations induced by Pk and Bv root colonization. (A) Co-expression network of differential markers (circled in blue); edges' width is proportional to the co-expression ratio between each gene according to RiceFREND database. (B) Gene expression dynamics of the co-expression network in response to Pk and Bv quantified by qRT-PCR. Transcript levels were normalized to that of reference gene EF 1α (Os03g0177400). Values between 0.5 and −0.5 are colored in gray. Data presented are the mean of log2 fold change (n = 3). JA, jasmonic acid; Bv, Burkholderia vietnamiensis TVV75 T ; Pk, Paraburkholderia kururiensis M130. differences observed between the two strains can be extrapolated to other species of the respective clades.
Another interesting conclusion from this study is the importance of the JA signaling in the interaction with beneficial rhizobacteria. As previously discussed, this major component of plant defense classically associated with the resistance to herbivores and necrotrophic pathogens appears to be involved in a larger diversity of biotic interactions (Thaler, 2004). We have demonstrated that there is a temporal delay in the induction of JA-related genes in leaves following the colonization by Pk comparatively to the response to Bv. It would be of interest to investigate the role of JA in the establishment of these interactions in terms of colonization level and plant defense status using JA-deficient mutants.

DATA AVAILABILITY
These sequence data for this study can be found in the EMBL database under accession number PRJEB31936.

AUTHOR CONTRIBUTIONS
EK, PC, and LM contributed to the conception and design of the study. EK, AW, IR, CB, and AK performed data collection and analysis. EK wrote the first draft. PC and LM contributed to manuscript revision. All authors read and approved the submitted version.

FUNDING
The authors acknowledge the CGIAR Research Program on Rice Agri-food Systems (RICE) and the MIC-CERES Project (FC Project ID 2013-1888; AF Project ID 1301-003, jointly supported by Agropolis Fondation through the "Investissements d'avenir" programme with reference number ANR-10-LABX-0001-01, and Fondazione Cariplo) for funding. EK and AW were supported by a fellowship from the French Ministry of Higher Education, Research and Innovation.

ACKNOWLEDGMENTS
We would like to thank Annette Vergunst for providing the pIN29 plasmid. We would like to thank Antony Champion and Trang Hieu Nguyen for providing the JAZ qPCR primers.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.01141/ full#supplementary-material FIGURE 7 | Local and systemic transcriptional regulations of rice in response to Pk and Bv root colonization. Common responses include the up-regulation of development-related genes and the down-regulation of defense-related genes in both leaves and roots. Strain-specific regulations in roots encompass the downregulation of different defense-related genes: While Pk induced the down-regulation of pathogenesis-related (PR) genes, Bv induced the down-regulation of WRKY transcription factors. Root colonization also induced a systemic up-regulation of jasmonic acid-related (JA) genes in leaves, however, not at the same time for each strain: While Bv induced this signal at early stages of the interaction, Pk induced a transient delayed systemic up-regulation of JA-related genes. Bv, Burkholderia vietnamiensis TVV75 T ; Pk, Paraburkholderia kururiensis M130.