S. meliloti transcriptome displays evidence for genotype-by-genotype interactions
Given the large number of phenotypes showing S. melilotistrain-specific differential response to Trichoderma species, we aimed to identify the transcriptome signatures (in terms of number and type of differentially expressed genes, DEGs) of such genotype-by-genotype interactions. The four S. meliloti strains were then exposed to 1W spent media for 24h and their transcriptome was evaluated by RNA sequencing.
The number of differentially expressed genes after treatment (here termed as stimulon) identified is shown in Table 1 . The list of DEGs can be found in Supplementary Datasets S2 and S3 . Overall, with respect to the four S. meliloti strains, the stimulon due toT. gamsii spent medium elicited the highest number of genes (average ca. 20.6%), while T. harzianum and T. tomentosumtreatments resulted in the lowest number of significant DEGs (from 6.0 to 10.0% of total genes) with limited differences between S. meliloti strains. T. velutinum , on the contrary, showed the ability to stimulate a high number of DEGs in 1021 and in the hybrid strain (ca. 23.8-25%), but in BL225C and AK83 the stimulon dropped to ca. 6.6-9.6%, indicating that a large part of the elicited transcriptome is strictly dependent upon the combination between givenS. meliloti genotypes and fungal genotypes. It is worth noting that the growth phenotypes under T. velutinum 1W spent medium resulted in a quite different response among the 4 S. melilotistrains (Figure 1 ).
Since S. meliloti harbours an open pangenome with a core set of genes shared by all strains and a dispensable set of genes present in a fraction only of strains we can expect that transcriptomic signatures of genotype-by-genotype interactions could be due to (mainly) the core set or the dispensable set. To clarify if and how much of a S. meliloti strain specific response is reflected in a transcriptional rewiring of core genes, we performed a principal component analysis (PCA) on the DEGs belonging to the core gene set (shared orthologs among the four S. meliloti strains) (Supplementary Dataset S4 ). Results (Figure 5 ) revealed that, for each strain, the transcriptional responses of the core genome to the different fungal spent media give rise to fungal-related grouping of samples, with the first principal component separating T. tomentosum /T. harzianum elicited transcriptomes from T. gamsii , the second principal component differentiating T. tomentosum from T. harzianum . This pattern resembles the separation based on spent media metabolome composition (Figure 2 ) and partial S. meliloti growth (Figure 1 ). Concerning up-regulated DEGs (Figure 5a ), the transcriptional response in all the strains treated with T. harzianum spent medium is very similar. This clustering is even tighter in the down-regulated genes (Figure 5b ). Similarly, T. tomentosum treatment elicited a similar response among the samples, both in up-regulated and down-regulated genes. However, AK83 and BL225C treated with T. velutinum spent medium were found in T. tomentosum clusters and, likewise, 1021 and the hybrid treated with T. velutinum were found clustering among samples exposed to T. gamsii spent medium, indicating that transcriptomic signatures of genotype-by-genotype interaction reside in the shared set of genes (core genome) also. Figure 6 reinforces the evidence that a large fraction of the transcriptome can be involved in the genotype-by-genotype interaction, in particular concerning strain-specific response to each fungal species. Indeed, most of DEGs identified are unique to the response to single fungal species (Figure 6a-d ) or at least can be shared among two fungal species (Figure 6 e-h ), while only a very limited fraction of DEGs is common to all fungal species. Interestingly, depending on the strain, the number of unique up-regulated genes belonging to the dispensable genome spanned from 40.4% to 69.9%, while for down-regulated genes, the dispensable range varied from 51% to 85%, suggesting the importance of the accessory genome in the strain- specific response. The lists of unique genes are reported inSupplementary Dataset S5 .
These data highlight the relevance of transcriptome variation in strain-specific molecular communication between soil fungi and rhizobia. Such variations also underlie the presence of regulatory interactions in the genome which are differentially affected by the presence of fungi. In this context, the transcriptional response of the hybrid S. meliloti strain (e.g. comparing the panels for 1021 and the hybrid strain in Figure 5 ) is different from that of the parental 1021 strain, though they share the same core genome (being different in the symbiotic megaplasmid pSymA only). We may consequently hypothesize the presence of epistatic interactions between genes on pSymA and those harboured by the chromosome and pSymB chromid related to the response to the presence of Trichoderma , similar to what has been found for the response to root exudates . This hypothesis broadens the relevance of pSymA megaplasmid, which could be related not just to the symbiotic interaction with the host plant (see for instance ), but also to the molecular dialogue with soil and rhizospheric fungi.
To inspect the type of genes functions affected by fungal spent media, we collapsed DEGs into Clusters of Orthologous Genes (COG) categories. Principal Component Analysis on this dataset (Figure 7 ) showed for the up-regulated genes, a relevant contribution of COG category K (transcription). This observation reinforces the hypothesis about the differential modulation of epistatic interaction in S. melilotigenomes triggered by fungal spent media, which may involve various regulons. Among down-regulated genes, an important contribution to variance was found for COG categories G (carbohydrate metabolism and transport) and J (translation), suggesting that fungal spent media may also have differential nutritional effects over S. melilotimetabolism. However, we should point out that the highest contribution to transcriptional variance was for genes not found in COG or with unknown function (COG category S), indicating that probably much of the genes relevant for the modulation of interaction taking place in the rhizosphere microbiota has still to be disclosed.
However, concerning specific functional genes with known relevance in plant-rhizobium interaction and which can explain part of the synergistic phenotypes observed (Figure 3 ), among the most highly expressed genes in S. meliloti 1021 and in thecis -hybrid strain after exposure to T . gamsii ,T . harzianum and T . velutinum , were components of the flagellar apparatus (flgA , C , D ,F , G , H , I , L and fliE ,G , I , K , L , M , N ), (Supplementary Dataset S2 ). It is worth noticing that the orthologs of these genes were not induced in BL225C or AK83. To reinforce the presence of epistasis and novel roles for pSymA megaplasmid, while for the cis -hybrid strain the same spent media strongly reduced root adhesion (Figure 1 ), in line with the expectation of an active flagellar apparatus,, the same was not true for 1021, which formed an abundant biofilm on the roots when treated withT. velutinum spent medium. However, regarding the general differences among S. meliloti strains, we must consider that the transcriptome was evaluated after 24h incubation, while the root biofilm after 48h and in the presence of the plant root system. Indeed, the motility-to-biofilm transition often includes an early phase with active flagellar apparatus, followed by an inhibition of flagellar gene transcription . We may speculate that the kinetics of this transition can be differentially modulated by Trichoderma species among the tested S. meliloti strains, however direct observations of biofilm formation under co-inoculation are needed to sustain such hypothesis.
No evidence was found for changes in the expression of galactoglucan (EPS-II) and succinoglucan (EPS-I) biosynthesis genes, but exoDgene was overexpressed in all the strains. This gene has been demonstrated to be needed for the invasion and development of alfalfa nodules .
Regarding genes involved in nodulation and nitrogen fixation, unexpectedly nodA expression was found induced in BL225C and 1021 in the presence of T. gamsii and T. velutinum , respectively, while nodD3 expression was induced in 1021 in the presence of T. tomentosum and T. velutinum . The gene was also induced in the hybrid strain when treated with T. velutinumand in AK83 in presence of T. gamsii . NodD3 does not require specific plant compounds to activate nod gene expression , but its expression is activated by a LysR family transcriptional regulator,syrM . This latter gene was indeed found overexpressed in the same conditions. This gene also activates syrA , which when overexpressed, causes an increase in EPS-I production. Strikingly,nifN and nifH expression was induced in AK83 in the presence of T. velutinum spent medium, questioning the possibility of activation of part of the nitrogenase complex in nonsymbiotic conditions. Concerning other genes which could be relevant for rhizobium phenotypes connected with symbiosis, among the up-regulated genes under spent media treatment those encoding for proteins involved in conjugal transfer (tra genes) were retrieved, suggesting that presence of the fungus may promote the mobilization of the pSymA megaplasmid. However, the rctA gene, the repressor of the conjugation machinery was found differentially expressed in 1021 and Bl225C strains under T. gamsii spent medium treatment (Log2 Fold change ~ 2), suggesting repression of the megaplasmid conjugation. Interestingly, the same gene was not differentially expressed in the cis -hybrid strain, indicating the presence of epistatic interaction among genes on the megaplasmid and those residing on the chromosome and the pSymB chromid.
o quantitatively estimates the number of DEGs whose expression changes among strains and treatments are explained by the rhizobial strain they belong to, by the fungal species and by the interaction between strain and species, a nested LRT model was constructed. Overall results are reported in Table 2 and the list of genes inSupplementary Dataset S6. A total of 327 genes were shared as DEGs among stimulons, which allowed to train a model for a total of 981 strain, fungus, and strain x fungus combinations. Within these genes, a large fraction (ca. 2/3 of total DEGs) showed a significant effect of strain, fungus, or strain x fungus, only 34.1% of their variability in expression not being supported by the models. For 146 genes changes in relative expression levels were modeled with respect to the rhizobial strain, for 270 to Trichoderma (giving a total of 416 combinations for both rhizobium and Trichoderma ), and for 230 genes to rhizobium x fungus interaction. Within this latter list of genes, several are related to redox balance (as sox gst , and nuo genes), suggesting that part of the genotype-species specific response could be related to a general differential response to stressful conditions.