Coupling the endophytic microbiome with the host transcriptome in olive roots.

The connection between olive genetic responses to environmental and agro-climatic conditions and the composition, structure and functioning of host-associated, belowground microbiota has never been studied under the holobiont conceptual framework. Two groups of cultivars growing under the same environmental, pedological and agronomic conditions, and showing highest (AH) and lowest (AL) Actinophytocola relative abundances, were earlier identified. We aimed now to: i) compare the root transcriptome profiles of these two groups harboring significantly different relative abundances in the above-mentioned bacterial genus; ii) examine their rhizosphere and root-endosphere microbiota co-occurrence networks; and iii) connect the root host transcriptome pattern to the composition of the root microbial communities by correlation and co-occurrence network analyses. Significant differences in olive gene expression were found between the two groups. Co-occurrence networks of the root endosphere microbiota were clearly different as well. Pearson's correlation analysis enabled a first portray of the interaction occurring between the root host transcriptome and the endophytic community. To further identify keystone operational taxonomic units (OTUs) and genes, subsequent co-occurrence network analysis showed significant interactions between 32 differentially expressed genes (DEGs) and 19 OTUs. Overall, negative correlation was detected between all upregulated genes in the AH group and all OTUs except of Actinophytocola. While two groups of olive cultivars grown under the same conditions showed significantly different microbial profiles, the most remarkable finding was to unveil a strong correlation between these profiles and the differential gene expression pattern of each group. In conclusion, this study shows a holistic view of the plant-microbiome communication.

Introduction

The olive tree (Olea europaea L.) is iconic in the Mediterranean Basin and has strong economic, social and ecological impacts in this region [1]. While the origin and domestication of this crop has not yet been undoubtedly documented, it is likely that took place in the Levant area from where olive cultivation gradually expanded along with human migrations [2], [3]. Multilocal selection and backcrosses between wild (O. europaea L. subsp. europaea var. sylvestris (Mill.) Leh.) and cultivated (O. europaea L. subsp. europaea var. europaea) olives have created the huge number of cultivars currently available [4]. Around 10.5 million ha are now dedicated to olive cultivation worldwide, creating an agro-ecosystem of utmost relevance in countries of the Mediterranean Basin where 98% of the cultivated area is located [5]. The major product of this tree crop is virgin olive oil, a source of fat providing nutritional and human health benefits that include reduction in the prevalence of cardiovascular diseases, prevention of some types of cancers and type 2 diabetes, and stimulation of higher diversity of beneficial gut bacteria [6], [7], [8].

In contrast to most of the fruit crops, the olive tree is well adapted to environmental conditions usually found in Mediterranean-type climatic regions (e.g. wide thermal amplitude, salt stress, water deprivation, etc.) [9], [10]. Tolerance to some of these conditions (i.e. drought and/or saline soils) appears to be genotype dependent although studies are only limited to few cultivars [11], [12], [13]. While the plant genotype is clearly decisive to explain adaptation to diverse environmental situations, the contribution of its associated microbiota is just starting to be unveiled. This is particularly true for the plant root/belowground microbiota pair since crucial multitrophic interactions for the development, fitness and health of the plant holobiont take place in the rhizosphere/root endosphere compartments [14]. However, the connection between the global host genetic responses to different environmental and agro-climatic conditions and the composition, structure and functioning of its inhabiting microbiota remains mostly unknown. While pedological characteristics, plant age and phenology, variable weather conditions and root exudates [15], [16], [17] are relevant, the plant genotype often appears as the most decisive factor to shape the composition and functioning of the belowground microbiome, particularly at the root endosphere level [18], [19], [20], [21], [22], [23].

Our knowledge about the link between olive-associated microbial communities and growth, development, fitness and adaptation to a/biotic constraints of this long-living tree is still very fragmentary. Nevertheless, recent studies are providing significant information about the factors shaping the olive-associated microbiota and how their constituents contribute to host stress tolerance and adaptation to specific soil/climatic conditions. For instance, different soil management practices affected the composition of the arbuscular mycorrhizal fungi (AMF) community [24]. The olive endophytic bacteriome has been shown to originate from the soil, reaching aerial organs through the xylem sap, and the composition of these communities seemed to be largely influenced by agronomic practices implemented in the orchard [25], [26]. Our own studies have demonstrated that agricultural management produces differences in the belowground microbiota compared to that in adjacent soils left under natural conditions. Agronomic practices also removed differences in the microbial communities associated with olive cultivars [27]. Other works have shown that olive cultivars played a more decisive structuring effect than plant organs on the composition of the epiphytic and endophytic bacterial communities of the phyllosphere [28]. Most of the available studies on the olive microbiota are related to biotic constraints caused by pathogens and arthropods such as Xylella fastidiosa, causing the olive quick decline syndrome [29], [30], [31], Pseudomonas savastanoi pv. savastanoi, the etiological agent of the olive knot disease [32], or Bactrocera oleae, the olive fruit fly [33]. Our own studies have revealed that the presence of the soil-borne fungus Verticillium dahliae, causing one of the most serious disease affecting olive cultivation (Verticillium wilt) [34], did not produce relevant alterations in the structure and functionality of the belowground microbiota of two olive cultivars differing in susceptibility to the disease. Remarkably, however, notable differences were found in their co-occurrence networks in response to the inoculation with the pathogen, particularly in the root endosphere communities [35].

Concerning the host side of the olive holobiont, the most recent whole transcriptome-based studies are providing interesting information on aspects such as how and where the domestication process might have taken place [36], the basis of the alternate olive fruit bearing [37], the complex gene expression patterns in different organs [38], the genetics controlling drupe maturation and response to environmental cues [39], or the development process from germination to the juvenile stage [40]. Likewise, global genetic responses to different stressing situations have been studied using transcriptomic approaches as for the case of olive adaptation to cold temperatures [41], [42] or salinity [43]. RNA-seq studies have been very useful to unravel the differential responses of olive cultivars upon the attack of pathogens such as X. fastidiosa [44] or V. dahliae [45]. Moreover, dual or co-transcriptomics approaches offered the opportunity to unveil the interaction between olive and V. dahliae at full transcriptome level and during the interaction between both partners [46], [47]. However, the cross-talk between olive and its inhabiting microbiota has not been studied from a holistic perspective. Recently, available RNA-seq datasets from root tissues of V. dahliae-infected and non-infected plants have been re-examined in an attempt to describe the disease through a systems biology approach [48]. While data mining from RNA-seq studies offers a novel and interesting perspective, the role of important components of the plant-associated microbiota to explain the reaction to a specific stress situation can be overlooked to some degree due to the methodological approaches implemented [48], [49].

We recently described the belowground microbiota (rhizosphere and root endosphere) associated to different olive cultivars grown under the same climatic, agronomical and pedological conditions [20]. Overall, results showed that the olive genotype was the main factor shaping the root-associated microbial communities, and that the olive root endosphere microbiota was less diverse than that of the rhizosphere. Remarkably, Actinophytocola (Actinobacteria, Pseudonocardiaceae family) was found to be the most abundant bacterial genus in the olive root endosphere, suggesting it may play an important role in olive adaptation to environmental stresses. Indeed, representatives of Actinophytocola spp. have been found in highly arid environments like the Sahara, Tenger Desert (China) and Mongolia [50], [51], [52], [53]. Interestingly, Actinophytocola was not only isolated from soils but also from inside the root of Thai glutinous rice plants (Oryza sativa sp.) [54].

Considering the previous antecedents, and to further explore the meaning of the significant differences in Actinophytocola relative abundance previously found in roots of different olive cultivars, we aimed in this study to: i) compare the olive root gene expression, at the whole transcriptome level, of two groups of olive cultivars harboring significantly-different abundances of the genus Actinophytocola; ii) analyze the rhizosphere and root-endosphere co-occurrence networks of the microbial communities associated to these two groups; and iii) connect the host whole-transcriptome expression to the composition of the root residing microbial communities by correlation and co-occurrence network (combining host relative gene expression and operational taxonomic units (OTUs) relative abundance data) analyses. The hypothesis to-be-tested is that diverse olive varieties cultivated under the same (agronomic, climatic, pedological) conditions show different root transcriptome patterns, and that these global gene expression differences can affect keystone OTUs (e.g. Actinophytocola) of the root endophytome that might be relevant for the fitness of the olive holobiont. This global interaction can be demonstrated by implementing an approach involving host transcriptomics and microbiomics.

Methods

Starting point data

The study of Fernández-González and co-workers [20], in which the microbial communities of the root system of 36 cultivars from the World Olive Germplasm Collection (WOGC) were studied, was taken as a starting point, using the same OTUs tables. In addition, transcriptomics data published by Ramírez-Tejero and co-workers [55] were also considered. Thus, the trees surveyed in these two previous studies constituted the biological samples (individuals) considered for this present study. In fact, both soil and root samples for metabarcoding and transcriptomic analysis originated from the same sampling event. Indeed, all samples were collected on April 26th, 2017 during the flowering period, after a month without rainfall (Additional File 1: Table S1) or irrigation. The only difference between the two studies was that for the metabarcoding study three replicates were taken per cultivar, while for transcriptomics there were two replicates. Therefore, the two replicates that matched in both studies were the ones used in the correlation analyses between genes and OTUs.

Selection of olive cultivars according to the differential relative abundance of Actinophytocola

Of the 36 cultivars previously mentioned, 10 were selected for this study because of their differences on relative abundance of the bacterial genus Actinophytocola previously reported [20]. We particularly focused on differences found in the root endosphere compartment. Thus, ‘Koroneiki’, ‘Mavreya’, ‘Myrtolia’, ‘Uslu’ and ‘Verdial de Vélez-Málaga’ were preliminary selected as the extreme representatives of the Actinophytocola high-prevalence (AH) group while ‘Barri’, ‘Klon_14_1812’, ‘Maarri’, ‘Mastoidis’ and ‘Shengeh’ were chosen as the extreme members of the Actinophytocola low-prevalence (AL) group.

Gene expression analysis

To compare gene expression profiles between the two groups of cultivars, AH and AL, the gene expression study was conducted using the DNAStar (ArrayStar 16) Qseq software for RNA-seq analyses (www.dnastar.com). PE Illumina reads were mapped using as reference the complete transcriptome of Olea europaea cv. Picual [56]. Mapping was undertaken using high-stringency parameters k-mer = 63 and 95% of matches. The default normalization method of reads per kilobase per million mapped (RPKM) was implemented. Results and associated charts were processed with Microsoft Excel. A fold change (FC) ≥ 8 and a p-value ≤ 0.01 were considered to detect differentially expressed genes (DEGs) between AH and AL. Thus, comparisons between groups consisted in 10 biological replicates of each group, corresponding to five cultivars and two trees from each one.

Statistical analyses

To assess the differences in Actinophytocola relative abundance between the AH and AL groups, both groups were compared by the Student’s t-test. In addition, α-diversity indices (observed richness; Shannon and inverse of Simpson diversity; and Shannon evenness) were compared with the Mann-Whitney-Wilcoxon’s U test using the R function ‘wilcox.test’. These comparisons were split and the number of sequences was rarefied by compartment (i.e. root endosphere and rhizosphere) and by domain (i.e. bacteria and fungi). To analyze β-diversity, a normalization of the filtered OTUs sequence counts was performed using the ‘trimmed means of M’ (TMM) method with the BioConductor package ‘edgeR’ [57]. The normalized data were considered to perform Non-Metric MultiDimensional Scaling (NMDS) on Bray-Curtis dissimilarities to ordinate in two dimensions the variance of β-diversity between the AH and AL groups. Similarly as for α-diversity, these comparisons were split by compartment and by domain. Ordination analyses were performed using the R package ‘phyloseq’ [58]. To analyze the effect of the AH versus AL grouping on community dissimilarities, permutational analysis of variance (PERMANOVA) and permutational analysis of multivariate dispersions (BETADISPER) were performed using the functions ‘adonis’ and ‘betadisper’ in the ‘vegan package’ with 9,999 permutations [59]. All of the above-mentioned analyses that were carried out using R have been provided in a script (Additional File 2: Statistical_analyses_script.R).

Pearson’s correlations and heatmap construction

In order to perform the correlation analysis between DEGs and microbial OTUs in the root endosphere of all cultivars under study (i.e. AH and AL groups together), gene expression values (reads per kilobase per million mapped reads; RPKM) were multiplied by 1000. Therefore, setting the smallest value, bigger than 0, to 1 in the dataset. Then, the replicate from which only microbial but not gene information was available was removed (see section 2.1). Eventually, two replicates were kept for each cultivar. Thus, cultivar ‘Shengeh’ was excluded due to the low number of sequences in the fungal community of this cultivar (n = 18). Subsequently, the data were log-transformed to be able to compare them. To eliminate spurious correlations, only those OTUs found in at least half of the samples were retained (n = 308 OTUs in n = 9 replicates). The correlation between DEGs and OTUs was made with the Pearson correlation coefficient using the ‘corrplot’ [60] package and a heatmap was constructed with the ‘pheatmap’ [61] package, both available in R. Pearson’s correlation coefficient cut-off of 0.6 and p-value < 0.05 were stablished to consider significant and strong correlations. To sort the heatmap, hierarchical clustering dendrograms were carried out for DEGs and OTUs with a complete linkage clustering method and using Euclidean distances with the software MultiExperimentViewer v.4.9.0 [62].

Network construction, comparison and visualization

Microbial (bacterial and fungal) co-occurrence networks were separately constructed for each group (AH and AL) and each compartment (root endosphere and rhizosphere) by using MENAP [63] and according to our previous implementation [35]. Three replicates were used per cultivar (n = 15, see section 2.1), except of AL endosphere, for which “Shengeh” had to be removed because insufficient fungal reads. Moreover, 100 random networks were performed to each empirical network for comparing the major topological properties between groups with Student’s T tests, as earlier described [64], [65], [66]. All the networks were visualized by using Cytoscape v.3.7.1 [67]. The node properties tables were analyzed in Microsoft Excel to plot and highlight the keystone OTUs according to their within-module connections (Zi) and among-module connections (Pi) values.

To generate the co-occurrence network between DEGs and OTUs, the log-transformed and filtered dataset previously obtained to construct the heatmap were used (see section 2.5). By using Pearson's correlation coefficient, a network with an R square of power-law value of 0.861 was obtained; that is, it gave rise to a robust and far from random network. A lower value was obtained when building the network using the Spearman's correlation coefficient. Due to this result Pearson was used as the best correlation method in both the heatmap and the co-occurrence network. The network construction method of MENAP web software is based on the comparison of the correlation matrix against an RMT (Random Matrix Theory) model. In summary, the RMT-based method informs whether correlations of the matrix (between OTUs and DEGs in our study) follow a Gaussian distribution, confirming a true correlation and discarding a Poisson distribution suggestive of an artefactual correlation. Eventually, this methodological approach predicts which cut-off value for Pearson is the ideal one to filter out all non-significant correlations (parameter that in other methods of network construction must be chosen by the user). In our case, the algorithm determined that the most appropriate cut-off Pearson’s value was 0.8, and not 0.6 which was used as a cut-off value to consider a strong correlation when constructing the heatmap (section 2.5).

Results

Differential distribution of endophytic Actinophytocola among olive cultivars

We earlier described the composition of the bacterial and fungal endophytic communities inhabiting olive roots from 36 cultivars comprising a subset of a so-called “olive core collection” [20]. Remarkably, data showed that the bacterial genus Actinophytocola (and others like Flavitalea although to a lesser extent) was found, by far, at significantly higher relative abundance in the olive root endosphere, with differences among cultivars though [20]. From this starting point, and in order to assess whether the observed significant differences in Actinophytocola relative abundance among olive cultivars correlate with differential host transcriptome responses, two groups (five cultivars each) displaying the highest and the lowest Actinophytocola relative abundances were made ad hoc. The relative abundance of Actinophytocola was determined by the proportional number of sequences for each cultivar replicate in the total number of bacterial reads. Five cultivars with the highest Actinophytocola prevalence constituted the AH group while five cultivars with the lowest Actinophytocola prevalence formed the AL group. The AH group had an average relative abundance of Actinophytocola of 42.38 ± 4.52 (mean percentage ± SE). In contrast, the AL group showed a mean prevalence of 6.01 ± 0.86 (t-test p-value = 1.30 × 10-08). The AH group included cultivars Uslu (Turkey), Verdial de Vélez-Málaga (Spain), Koroneiki (Greece), Mavreya (Greece) and Myrtolia (Greece), while the AL group was comprised by cultivars Shengeh, (Iran), Mastoidis (Greece), Barri (Syria), Klon_14_1812 (Albania) and Maarri (Syria) (Fig. 1). No correlation was observed between the microbial communities of AH and AL cultivars and their geographical origins, nor with tolerance/susceptibility to biotic stresses such as Verticillium wilt [68] or Anthracnose (caused by Colletotrichum species) [69] (Fig. 1). Moreover, AH and AL cultivars distributed randomly between the two different phylogenies of cultivated olives recently reported [56].

Fig. 1

Relative abundance of the genus Actinophytocola in olive cultivars of the AL and AH groups. AH: Actinophytocola high-prevalence group; AL: Actinophytocola low-prevalence group. Tolerance (+) or susceptibiltiy (−) of olive cultivars to Verticillium dahliae and Colletotrichum spp. are indicated. When this information was unavailable, cultivars are marked with “?”.

The AH and AL olive cultivars show significantly different whole-transcriptome profiles

In order to assess whether the Actinophytocola prevalence in the AH and AL groups was associated to different transcriptome profiles in the host roots, the existence of DEGs among these groups was investigated. This analysis was performed using high stringency parameters to detect only highly consistent DEGs. As a result, differences in gene expression between the AH and AL groups at 8 FC and p-value ≤ 0.01 were found in 60 genes (Additional File 3: Fig. S1), 39 up-regulated in the AH group (Table 1 and Additional File 4: Fig. S2) and 21 up-regulated in the AL group (Table 2 and Additional File 5: Fig. S3). In the set of 39 up-regulated genes in the AH group and down-regulated in the AL cultivars three genes involved in plant defense to pathogens were found, one coding for a Germin-like protein subfamily 1 member 16 (Oleur061Scf2041g01041.1), and two tlp genes encoding Thaumatin-like proteins (Oleur061Scf0191g00025.1, Oleur061Scf0201g06012.1) with antifungal activity that are secreted to the apoplast. Additionally, one CIGR1 gene (Oleur061Scf2742g00025.1) coding for a chitin-inducible gibberellin (GA)-responsive protein 1, that may play a regulatory role in the early step of oligosaccharide elicitor response [70], and a response-to-stress ABA2 gene (Oleur061Scf2621g00015.1) encoding a zeaxanthin epoxidase, were also up-regulated in AH and down-regulated in AL cultivars. In addition, five genes (ATL22, ATL46, ATL8, PP2B11 and PCYOX1; corresponding to Oleur061Scf0789g04009.1, Oleur061Scf4176g04002.1, Oleur061Scf1016g05039.1, Oleur061Scf7684g00004.1, Oleur061Scf3403g03008.1, respectively) involved in protein ubiquitination and/or protein degradation were also identified in this set of genes. In addition, four genes encoding proteins involved in protein biosynthesis pathway (Utp12, L34Ae, pre-mRNA-splicing factor 38 and protein translation factor SUI1; corresponding to Oleur061Scf7124g01014.1, Oleur061Scf9068g00010.1, Oleur061Scf0172g00010.1, Oleur061Scf0579g00003.1, respectively) were also identified in this set, suggesting a higher protein turnover in AH roots. In this sense, gene expression was particularly activated in AH roots by two up-regulated transcription factors (TF), HAT5 and HOX6 (Oleur061Scf7420g01005.1, Oleur061Scf2236g05008.1), both coding for Homeobox-leucine zipper proteins. The HOX6 gene may be involved in drought response [71]. In addition, a XRN3 gene (Oleur061Scf1132g02001.1) coding for a suppressor of post-transcriptional gene silencing was up-regulated in AH roots as well. Another difference between AH and AL root transcriptomes was the expression of genes involved in the cell wall. Thus, in AH roots, two genes coding for Extensin proteins (Oleur061Scf3325g00004.1, Oleur061Scf2628g03012.1) that form a structural component of the primary cell wall and strengthen it, were found up-regulated. Related to this, a WAK2 gene (Oleur061Scf3771g00037.1) encoding a wall-associated receptor kinase 2, probably acting as a signaling receptor of extracellular matrix component and required for cell expansion [72], was also up-regulated in AH cultivars. Remarkably, two probable plant transposases (Ptta/En/Spm families) (Oleur061Scf6180g00005.1, Oleur061Scf1200g04004.1) were up regulated in AH cultivars, suggesting an elevated transposon activity of these families of transposable element in AH roots. Finally, other genes related to membrane traffic, respiration, oxidoreductase activity, development or just unknown function were up-regulated in the roots of the AH group (Table 1).

Table 1

Upregulated genes in the Actinophytocola high-prevalence (AH) group of olive cultivars.

Genes with higher expression in the AH group of cultivars at 8 FC and p-value ≤ 0.01. Genes that may be involved in the same process or same compartment are represented with the same color. If processes are not the same but related, the same color with different tonality was used. Colorless means that the process is not known.

Table 2

Upregulated genes in the Actinophytocola low-prevalence (AL) group of olive cultivars.

Genes with higher expression in the AH group of cultivars at 8 FC and p-value ≤ 0.01. Genes that may be involved in the same process or same compartment are represented with the same color. If processes are not the same but related, the same color with different tonality was used. Colorless means that the process is not known.

From the set of 21 genes up-regulated in AL roots and down-regulated in the AH group, two ethylene (ET)-responsive TF, encoded by ERF4 and ERF34 genes (Oleur061Scf9036g01047.1, Oleur061Scf1024g00024.1), that binds to the GCC-box pathogenesis-related promoter element, and one ET-responsive element, the ACO3 gene (Oleur061Scf1251g02004.1) coding for an aminocyclopropane-1-carboxylate oxidase involved in the biosynthesis of the plant hormone ET, are worth to be highlighted. Additionally, in this set of genes a TF (GAI1) (Oleur061Scf5473g04020.1) that acts as a repressor of the GA signaling pathway was identified. The SHR (SHORT-ROOT) gene (Oleur061Scf3289g01014.1), that encodes another TF, was found to be up-regulated in AL roots as well. This gene regulates the radial organization of the root in Arabidopsis [73]. The LBD30 gene (Oleur061Scf0088g04051.1), also up-regulated in AL cultivars, is involved in root development. In fact, overexpression of this gene produced shrunken root tips and disorganized columella cells in Arabidopsis [74]. Likewise, the GSO1 gene (Oleur061Scf3112g01022.1) that codes a LRR receptor-like serine/threonine-protein kinase is involved in root growth. This gene regulates root growth through control of cell division and cell fate specification in Arabidopsis [75]. Finally, the CSLC4 gene (Oleur061Scf2527g00016.1), coding for a xyloglucan glycosyltransferase 4 involved in synthesis of the non-cellulose component of the cell wall was also identified in this set of 21 up-regulated genes.

Effects of cultivar grouping in the distribution of the root microbiota

As mentioned above, the criterion to differentiate between AH and AL cultivars was the significant dominance of the genus Actinophytocola in the root endophytome of the AH group. This dominance was statistically proven when analyzing the α-diversity of the endophytic bacterial community. Interestingly enough, the same (i.e. lower richness and diversity indices in the AH group than in AL cultivars) was also observed for the endophytic fungal community. However, no differences were found in the rhizosphere microbial communities (Table 3, Additional File 6: Fig. S4).

Table 3

Comparisons of α and β-diversity indices between microbial communities of the Actinophytocola high-prevalence (AH) and low-prevalence (AL) groups.

	α-diversity				β-diversity
	Richness	Shannon	InvSimpson	Evenness	PERMANOVA	BETADISPER
Bacteria
Endosphere	2.46E-03	1.64E-02	2.35E-02	2.63E-02	1.00E-04	1.60E-02
Rhizosphere	5.57E-01	4.61E-01	9.99E-01	5.12E-01	1.00E-04	6.44E-01
Fungi
Endosphere	4.24E-03	1.73E-03	1.39E-02	1.27E-01	5.10E-03	4.69E-01
Rhizosphere	5.75E-01	3.05E-01	2.33E-01	2.85E-01	1.00E-04	2.70E-02

Mann-Whitney-Wilcoxon tests were performed to evaluate α-diversity. β-diversity was compared using Bray-Curtis dissimilarities. Significant p-values are highlighted in bold type and italics.

Surprisingly, when comparing the rhizosphere microbial (bacterial and fungal) composition by PERMANOVA and NMDS of the AH and AL cultivars (Table 3, Additional File 7: Fig. S5) both domains showed statistically different communities, an outcome a priori expected only for the endosphere. Therefore, the rhizosphere microbial communities of the AH and AL groups showed similar richness and α-diversity values but with statistically distinguishable microbial assemblies (β-diversity).

The belowground microbiota of the AH and AL groups show different co-occurrence network topologies

Concerning the rhizosphere microbiota, similar co-occurrence network properties were obtained between AH and AL groups. The most remarkable difference was the much larger intra-modular connectivity (formerly known as ‘centralization of stress centrality’ and shown in Table 4 as CS) of the AH network (3.082) compared to that of the AL network (0.579). That is, the AH network had greater number, and more interconnected, module hubs (Table 4 and Fig. 2).

Table 4

The major topological properties of co-occurrence networks of Actinophytocola high-prevalence (AH) and low-prevalence (AL) groups.

Compartment	Group	Total nodes	Total links	Percentage of positive edges (PEP)	Similarity threshold (St)	R square of power-law	Average degree (avgK)	Average clustering coefficient (avgCC)	Average path distance (GD)	Modularity (M)	Centralization of stress centrality (CS)	Transitivity (Trans)
Endosphere	AH	129	491	97.56%	0.72	0.849	7.612	0.305	2.890	12 (0.380)	1.013	0.286
	AL	233	272	65.81%	0.83	0.903	2.335	0.117	6.247	34 (0.834)	0.848	0.149
Rhizosphere	AH	845	984	90.75%	0.85	0.904	2.329	0.089	8.581	120 (0.861)	3.082	0.089
	AL	758	948	91.56%	0.85	0.894	2.501	0.099	8.722	115 (0.857)	0.579	0.150

Figures in bold type and italics show statistically-significant differences between AH and AL groups in each compartment (root endosphere and rhizsophere) by Student’s T-test.

Fig. 2

Co-occurrence networks of microbial communities of each compartment (root endosphere and rhizosphere) and each group (AH and AL). Each node corresponds to an OTU and were colored according to their phylum affiliation. Arrowheads represent module hubs and rhombus represent connectors. Lines represent positive (green) and negative (red) connections. AH: Actinophytocola high-prevalence group; AL: Actinophytocola low-prevalence group. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Regarding the root endosphere community two clearly different co-occurrence networks were obtained. Indeed, the AH network was more complex (higher avgCC and lower GD) and with much lower M value). As a whole, these indices were suggestive of a more interconnected network in the AH group. Additionally, the ‘percentage of positive edges’ (PEP) parameter was much lower in the AL group, indicating a higher proportion of negative interactions in this network (Table 4 and Fig. 2). Furthermore, an OTU belonging to Actinophytocola was a keystone node in the AH network, this OTU being the only module hub in the largest module of this network. In contrast, it was just a peripheral node in the AL network (Additional File 8: Fig. S6).

Connecting the olive root transcriptome with the endophytic microbial community

In light of the results so far presented, we investigated a potential correlation between DEGs and the different microbial profiles found in the root endosphere of AH and AL cultivars. Correlation analysis between the root endosphere microbial community and the set of DEGs previously identified in the AH and AL groups (8 FC, p-value ≤ 0.01; Table 1, Table 2) showed a strong (−0.6 ≤ r ≥ 0.6) and significant (p < 0.05) correlation with 111 OTUs. Two clearly different profiles could be observed (Fig. 3). On the one hand, a cluster with 82 OTUs belonging to nine different phyla, plus some unclassified bacteria and fungi, showed negative correlation with the up-regulated genes in the AH group (first row, yellow squares; Fig. 3) and positive correlation with the up-regulated genes in the AL group (first row, blue squares; Fig. 3). On the other hand, another cluster comprising 29 OTUs displayed the opposite correlation profile. Remarkably, 25 out of 29 OTUs in this cluster belonged to the phylum Actinobacteria, eight of them identified as Actinophytocola and eleven additional OTU unclassified at the genus level but assigned to the same family (i.e. Pseudonocardiaceae) as Actinophytocola.

Fig. 3

Heatmap with Pearson’s correlations between the selected DEGs and the strongly and significantly correlated endophytic OTUs. Both OTUs and DEGs (see Table 1, Table 2 for functional annotation) were clustered with a complete-linkage clustering method and using Euclidean distance. bOtu’s represent bacterial OTUs and fOtu’s represent fungal OTUs. Strong correlation was considered when absolute value of Pearson correlation coffecient was ≥0.6.

The Pearson’s correlation analysis enabled to draw a first approximation to the interaction taking place between the root host transcriptome and the endophytic community. To go one step further in the search for keystone OTUs and genes in such interaction, a co-occurrence network analysis between the 60 DEGs and all OTUs observed in the root endosphere was performed. This analysis showed significant interactions between 32 DEGs (18 up-regulated and 14 down-regulated in the AH group, yellow and blue circles, respectively) and 19 OTUs (Fig. 4). Among all interactions observed, some of them were worth mentioning. For instance, OTU00021, belonging to Actinophytocola, showed a positive correlation with three up-regulated genes in the AH group (Oleur061Scf2041g01041.1, Oleur061Scf7501g00024.1 and Oleur061Scf3782g00023.1). The first one is related to response to biotic stress (functionally annotated as ‘Germin-like protein subfamily 1 member 16′) while the other two with unknown functions (annotated as ‘unknown protein’ and ‘Leishmanolysin-like peptidase’) (Fig. 4 and Table 1). The same OTU showed a negative correlation with two down-regulated genes in the AH group (Oleur061Scf1251g02004.1 and Oleur061Scf0724g05015.1, functionally annotated as ‘1-aminocyclopropane-1-carboxylate oxidase homolog’ and ‘Zinc finger, RING/FYVE/PHD-type’, respectively), the first one being key in the plant hormone ET biosynthesis pathway (Fig. 4 and Table 2). Furthermore, another OTU belonging to Actinophytocola (OTU00124) showed a positive correlation with two upregulated genes in the AH group (Oleur061Scf0191g00025.1 and Oleur061Scf0201g06012.1, both functionally annotated as “Thaumatin-like protein”). These two genes are related to response to biotic stress and antifungal activity (Fig. 4 and Table 1). The overall scenario shown in Fig. 4 is a negative correlation between all up-regulated genes in the AH group and all OTUs, except those ones belonging to the genus Actinophytocola. In contrast, the down-regulated genes in AH group showed just the opposite pattern.

Fig. 4

Co-occurrence network showing the interactions between DEGs and endophytic OTUs. Only significant correlations were retained (−0.8 ≤ r ≥ 0.8) according to the Ramdom Matrix Theory (RMT) method (see co-occurrence network construction methodology). Circles represent up-regulated (yellow) and down-regulated (blue) genes in the AH group. Round squares and triangles represent bacterial and fungal OTUs, respectively, colored according to their deepest taxonomical annotation. Lines represent positive (green) and negative (red) connections. Gene-gene and OTU-OTU correlation lines were eliminated for clarification purpose. The main topological properties of the complete network were: Total nodes = 166, Total links = 551, PEP = 71.87%, St = 0.8, R2 of power-law = 0.861, avgK = 6.639, avgCC = 0.329, GD = 3.636, M = 17 (0.471), CS = 1.812, and Trans = 0.461. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Discussion

The interaction of the olive root transcriptome with the belowground microbiota has never been studied considering the holobiont context. Numerous studies have shown that the host genotype is decisive to explain the structure and composition of its resident microbial community [76], [77]. However, a range of soil and environmental parameters, as well as diverse stresses, also affect the holobiont. Indeed, external factors influence the host genetic responses what in turn may determine decisive structural and functional changes in the associated microbiome [78], [79]. Therefore, the global cross-talk established between the plant and its microbiome is not only crucial for the growth, development and fitness of the holobiont but also to efficiently cope with different (a)biotic constraints that must continuously face [15]. However, to our knowledge, this holistic approach has so far seldom accomplished, and available studies just focus either on the host side or on its associated microbiota [75], [76], [77], [78], [79]. In the case of olive trees, host whole-transcriptome responses have been analyzed under scenarios such abiotic [41], [80] or biotic [46] stresses, or to unravel their effects over specific organs or tissues [45]. Concerning the olive microbiome, several studies have aimed to address questions such as the link between host tolerance/susceptibility to biotic stress and the structure and functioning of the residing microbial community, or how the latter reacts against the presence of specific pathogens or insects [32], [33], [35]. One step further is to conduct dual or co-transcriptomic studies, but they mostly deal with a specific binomial as, for instance, the host transcriptome in response to an interaction with a given pathogen and vice versa [46], [47]. However, studies in which the entire microbiota is evaluated in connection with the host global genetic response are lacking. Recently, metatranscriptomics (based on data mining of preexisting RNA-seq studies) has been implemented to offer a more global overview of the interaction between the indigenous microbial communities and the host genetic response, including the case of olive [48]. While these studies are providing relevant information, results might be biased and some relevant information overlooked because of limitations on the methodological approach implemented [49], [81], [82].

In the present study, we connect for the first time the olive root transcriptome to the belowground microbiota, aiming to identify host genes and microbial OTUs potentially crucial in the adaptation of the olive holobiont to a specific natural scenario. Our starting point was the available information on the structure and composition of the olive-associated belowground microbiota at the WOGC [20]. In this latter study, we identified Actinophytocola as the most prevalent bacterial genus among the surveyed cultivars, which has been found linked to stress by water-supply deficiency. It must be emphasized that all cultivars in that study were grown under the same agronomical, pedological and climatic conditions, and that they underwent a prolonged period without water supply that might be related to the high prevalence of Actinophytocola genus in the roots (Additional File 1: Table S1). However, the Actinophytocola relative abundance showed large variations among cultivars. Under these conditions, the genotype was the determinant factor explaining the differences found in the composition of the belowground microbiota among cultivars. This led us to inquire whether cultivar-specific genetic responses (whole-transcriptome level) could be linked to the composition/structure of the microbial communities. Moreover, we also aimed to unravel differences in the network co-occurrence topologies of the outside (rhizosphere) and inside (root endosphere) communities of the two groups of cultivars selected. Finally, we investigated the feasibility to couple host gene expression profiles with the relative abundance of specific OTUs.

Our results demonstrated a correlation of the olive root microbiome with specific DEGs of the host. In fact, to our knowledge, this is the first study in which the interaction between specific host genes and precise root-inhabiting microbes is described. García and co-workers [81] carried out a similar study in the phyllosphere, although they used a different methodological approach (i.e. only transcriptomics data). Furthermore, a recent study used metatranscriptomics data to focus on the influence of the microbiome on differential gene expression in the tomato root endosphere [49]. Unfortunately, however, this work did not provide details on which genes and microorganisms were interacting. Nevertheless, our results agreed with those of García and co-workers [81] in that the phylum Actinobacteria seems to play a key role when the plant holobiont is confronting a stressing situation. According to them, Streptomyces (phylum Actinobacteria) is particularly abundant in plants subjected to drought stress. Although this genus was present in the correlation analysis and positively correlated with a stress response gene (see OTU00023 and OTU00102 and gene Oleur061Scf1329g05010.1 in Fig. 3), we did not find Streptomyces in the co-occurrence network (Fig. 4).

The roots of AH and AL cultivars displayed relevant differences in their global gene expression patterns, including genes related with plant defense response to stress. Strikingly, two Thaumatin-like coding genes, closely related with plant defense through apoplastic antifungal activity, showed higher expression in the AH group. Indeed, there was positive co-occurrence between these two genes and an Actinophytocola OTU (Fig. 4). This might be linked to a competitive relationship between fungi and bacteria, since the antifungal effect of thaumatin-like proteins might help Actinophytocola to flourish. Moreover, Actinomycetales are long known for their antagonistic behavior against root pathogenic fungi [83]. Nevertheless, it is very difficult to conclude whether the abundance of Actinophytocola is cause or consequence of the observed higher expression of Thaumatin-like genes. Interestingly, a higher number of genes involved in cell wall synthesis or strength were also differentially expressed in AH roots. On the one hand, AH roots showed a higher transcriptional activity, as suggested not only by the higher expression of genes encoding proteins linked to protein turnover, but also by the upregulation of two Homeobox-leucine zipper TF. These TF play significant roles in the response to abiotic stress [71], [84], [85]. On the other hand, AL roots overexpressed several genes involved in root growth or development, sugar and water transport, lipid metabolism and pigment synthesis.

Regarding the root microbiome side, it is worth mentioning that a synchrony between endophytic bacteria and fungi was observed. Indeed, the opposite trend is frequently reported in the literature; that is, bacterial and fungal communities associated to the plant roots behave independently, at least in terms of α-diversity, both in the rhizosphere [86] and the endosphere [87], [88]. In fact, this lack of synchrony was previously observed in the microbial communities of the olive root endosphere when comparing two cultivars that were challenged with V. dahliae [35]. Interestingly, and in contrast with most of the studies, we found that the endophytic bacterial and fungal communities displayed lower richness and diversity in AH roots than in AL cultivars. This difference, at the bacterial level, was mainly due to the overwhelming dominance of the genus Actinophytocola. However, in the case of fungi, no significant difference in the evenness index between the two groups was observed (Table 3), which led to a decrease in diversity solely due to a reduction in the richness of fungal OTUs in AH varieties. Regarding the rhizosphere, no differences in richness and α-diversity was found between the two groups of cultivars. Surprisingly, differences in β-diversity were observed in both the endosphere and the rhizosphere for the two microbial domains (Table 3, Additional File 7: Fig. S5). This result indicated that although the rhizosphere communities of both groups (AH and AL) were similar in richness and diversity, they showed significantly-different microbial compositions.

With regard to the interactions established within the microbial communities of each group of cultivars, it is worth mentioning the striking differences observed between network topologies of the endophytic communities (Fig. 2). The co-occurrence network of the AH group was more complex, with the individuals more interconnected (higher avgK and avgCC). In addition, these connections were stronger (lower GD and M) than in the AL network (Table 4). Considering the predominance of Actinophytocola in the AH group and the precedents on the adaptive capabilities of this genus to water stress conditions [50], [51], [53], the topological parameters found in this network suggested that the communities of each group responded to water stress (at least under the conditions found in the orchard at the time of sampling) in a very different way. In AH cultivars, the whole community seemed to work coordinately to adapt quickly to those stressing conditions. Accordingly to previous studies, higher and stronger connectivity (high avgK and avgCC linked to low GD) yielded networks with faster signal transmission among members of the community, and more rapid adaptation to any disturbance in the habitat [89], [90], [91]. Furthermore, this type of tightly-structured network (low GD and M) restricted the entry of new individuals with potential to alter it [92]. Therefore, such networks are generally more resistant and resilient to a/biotic stresses. In contrast, in the AL group a more heterogeneous topology was observed with larger distancing and compartmentalization of the community (high GD and M) [92], [93]. We previously found this kind of network when the effect of the stressor was greater and the community aimed to “confine” it [35]. In fact, some microbial groups (modules in the network) responded in the opposite way to stress (thus a decrease in PEP was observed, Table 4). Therefore, in this water stress scenario, results seemed to point to an increased competition among the members of the root endophytome in AL cultivars.

Conclusions

This study unveiled a significant interaction between specific genes up/down-regulated in olive roots and particular endophytes inhabiting this organ when trees were grown in the same agronomical and pedological scenario, and subjected to a period of water scarcity. Thus, a holistic view of the plant root-endophytic microbiome communication is here reported under specific environmental conditions. The roots of the two groups of olive cultivars examined showed significantly-different microbial profiles that strongly correlated with differential host gene expression patterns. The root transcriptome-microbiome “global communication” here described was solidly supported by both correlation and co-occurrence interaction analyses between host DEGs and microbial OTUs (Fig. 3, Fig. 4), and showed a close coordination between a group of olive genes that negatively interacted with most members of its root endophytome. Moreover, these genes significantly decreased their expression in cultivars with the lowest relative abundance of the genus Actinophytocola, which seemed to play a key role. Our findings also suggested that the gene expression profile in roots of AL cultivars “facilitated” a more permissive endosphere environment for microbial colonization, thereby leading to more diverse and richer microbiome than the one found in the AH group. However, it also generated a more competitive environment (Table 4 and Fig. 2) for the members of this microbiome under the stressing conditions found at the time of sampling. In contrast, the AH group showed more restrictive conditions in its root endosphere what could have favored, among others, the flourishing of Actinophytocola under the observed conditions. Overall, the methodological approach here implemented can be useful to generate a specific “host transcriptome-microbiome signature” for any plant holobiont subjected to a particular environmental scenario. In future studies, the identification of specific olive root exudates and/or the evaluation of structural/biochemical modifications of the roots under specific stressing situations, among other approaches, will provide a more precise view on how the plant genotype affects its associated microbiome.

Availability of data and materials

The datasets analyzed during the current study are available in the NCBI Sequence Read Archive (SRA) under the BioProject number PRJNA498945 (olive belowground microbiota), and in the Gene Expression Omnibus (GEO) repository, accession number GSE152236 (RNAseq).

Funding

Supported by grants AGL2016-75729-C2-1-R and AGL2016-75729-C2-2-R from the Spanish Ministerio de Economía, Industria y Competitividad/Agencia Estatal de Investigación, and grant PID2019-106283RB-I00 from the Spanish Ministerio de Ciencia e Innovación/Agencia Estatal de Investigación, all of them co-financed by the European Regional Development Fund (ERDF). These funding sources had no involvement in this work preparation.

Author’s contributions

JMB, FL and MFL conceived and designed the study. AJFG, FL, and MFL implemented the methodology used. AJFG, JRT, PNB and FL performed the bioinformatics analysis and analyzed the data. All authors contribute to the preparation of the original draft. AJFG, FL and JMB wrote, reviewed and edited the final version of the manuscript. Funding acquisition by FL, MFL and JMB. All authors read and approved the final manuscript.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.