Inferring co-expression networks of Arabidopsis thaliana genes during their interaction with Trichoderma spp.

Fungi of the Trichoderma genus are called "biostimulants" because they promote plant growth and development and induce disease resistance. We used conventional transcriptome and gene co-expression analyses to understand the molecular response of the plant Arabidopsis thaliana to inoculation with Trichoderma atroviride or Trichoderma virens. The transcriptional landscape of the plant during the interaction with these fungi showed a reduction in functions such as reactive oxygen species production, defense mechanisms against pathogens, and hormone signaling. T. virens, as opposed to T. atroviride, was more effective at downregulating genes related to terpenoid metabolism, root development, and chemical homeostasis. Through gene co-expression analysis, we found functional gene modules that closely link plant defense with hypoxia. Notably, we found a transcription factor (locus AT2G47520) with two functional domains of interest: a DNA-binding domain and an N-terminal cysteine needed for protein stability under hypoxia. We hypothesize that the transcription factor can bind to the promoter sequence of the GCC-box that is connected to pathogenesis by positioned weight matrix analysis.


Figure 1. Differential expression analysis. (A) MA-plot of the differential expression analysis
A.thaliana (T.atroviride) .On y-axis is show the log 2 fold changes versus x-axis represents log 10 the average expression signal.Red dots represent up-regulated genes; blue dots represent down-regulated genes.(B) MA-plot of the differential expression analysis A.thaliana (T.virens) .On y-axis is show the log 2 fold changes versus x-axis represents log 10  We observed that 138 genes increased their expression in response to both Trichoderma species.In addition, 118 and 125 genes were upregulated only in response to T. atroviride or T. virens, respectively (Figs. 1C,D; S2; Table S1).Using the logFC of the group of 138 upregulated genes in both interactions, applied a t-test to determine if any of the two Trichoderma species was more effective in modulating differential expression in the plant (Fig. 1C,D; Fig. S2a).The t-student test suggested that both fungi similarly modulate differential gene expression (p-val = 0.055, significance 95%).Similarly, we observed 37 genes that decreased their expression in response to both Trichoderma species.In addition, whereas 12 genes decreased their expression in response to T. atroviride and 89 in response to T. virens (Fig. 1B,D; Table S1; Fig. S2).The t-student test for downregulated genes suggested a difference in gene repression, with T. virens being more effective than T. atroviride (p-val = 0.02, significance 95%).
The Gene Ontology (GO) enrichment analysis showed the induction of phytohormone pathways such as ethylene, jasmonate, and salicylic acid (Fig. 2; Table S2).Additionally, the induction of mechanisms related to plant defense, such as acquired systemic resistance (GO:0009627), immune system process (GO:0002376), and response to molecule of bacterial origin (GO:0002237), among others, was observed (Fig. 2, Table S2).Mechanisms involved in the regulation of reactive oxygen species scavenging mechanisms, such as the glutathione metabolic process (GO:0006749) and response to hydrogen peroxide (GO:0042542), were also enriched.Finally, among the repressed GO terms: root development (GO:0048364), root morphogenesis (GO:0010015) and chemical homeostasis (GO:0048878) (Fig. 2, Table S2).

Differential network analysis
As expected, the differential network analysis showed changes in the topology of the interaction networks (network (T.atroviride) and network (T.virens) ) against their controls (Fig. 3A,B), probably because of the association of genes under different conditions.When the connectivity of the nodes using the degree as a metric was analyzed, the control networks exhibited a right asymmetry in both cases that extended to 150 and 300 degrees for the T. atroviride and the T. virens network, respectively.Although the control and the stimulus networks were built with the same nodes, we observed that the control networks had a much greater number of edges in both cases.The T. atroviride network contained 9826 edges, whereas its control presented 33,222.On the other hand, the T. virens network showed 43,784 edges, while its control had 116,682.
The evaluation of the weights of edges for each network and their respective controls showed that the edges of the interaction networks have greater weight than the edges of the control networks (Fig. 3C).The edges of the T. virens network have a correlation coefficient (r) greater than r s ≥0.95.Similarly, the edges of the T. atroviride network have an r s ≥0.9; finally, the edges of the control networks have an r s ⩾0.85.

Gene co-expression network analysis
We generated two gene co-expression networks of A. thaliana when stimulated with T. atroviride or T. virens.The resulting networks consisted of 1094 and 2011 nodes for Arabidopsis (T.atroviride) and Arabidopsis (T.virens) , respectively (Fig. S4a,b).The networks were simplified, as explained in Fig. S4.
The community identification analysis found 104 communities for the A. thaliana (T.atroviride) network and 256 for the A.thaliana (T.virens) network (Fig. S5).We focused only on those communities formed by at least 100 nodes and generated a functional enrichment analysis for communities Ca9, Ca16, and Ca18 for A.thaliana (T.atroviride) and Cv5, Cv14, Cv22 and Cv41 for A.thaliana (T.virens) (Figs.4B & 5B; Table S2).
Finally, we selected the HUB nodes using a score based on degree, authority score, and closeness to refine the network (Fig. 4A) and selected the best 50 nodes according to this score from each of the communities (Ca9, Ca16, Ca18 (150 total nodes) and Cv5, Cv14, Cv22 y Cv41(200 total nodes).When we removed the HUB nodes from the network (Fig. S6a), we observed that the network becomes less cohesive and an increase in the geodesic path between nodes.On the other hand, when the nodes considered non-relevant are eliminated, their general cohesion measures increased, indicating that our methodology to identify relevant nodes was adequate.
In the HUB nodes, we identified the AT2G47520 locus in both networks (Arabidopsis (T.atroviride) and Arabidopsis (T.virens) ), which encodes a transcription factor necessary to cope with hypoxia in A. thaliana 28 (Fig. S7).This transcription factor has two functional domains (AP2/ERF and Cys2) that suggest a dual role: mitigating hypoxia and coping with pathogen stress.The AT2G47520 locus interacted with 22 genes in the Arabidopsis (T.atroviride) network, while in the Arabidopsis (T.virens) network, it directly connected with 35 genes.Genes related to hypoxia, jasmonate, and ethylene hormonal pathways, defense mechanisms, and uncharacterized proteins were identified within these groups.

Experimental validation
We focused on validating our transcriptome and hypotheses using the A.thaliana-T.atroviride system.However, considering that both networks have the same statistical parameters, we expected them to behave similarly between the high throughput sequencing and the experimental validation.To validate our transcriptomic analysis, we selected 12 genes and compared the similarities between our differential expression results obtained by transcriptome and RT-qPCR (Fig. 6).During this analysis, we observed a similar behavior between the expression profiles obtained in the transcriptome and those obtained by RT-qPCR.As expected, in some cases, the same trend was observed for the transcriptomic data and the data obtained by RT-qPCR, as in the cases of the AT5G44420 (ρ s = 0.98), AT2G47520 (ρ s = 0.95), AT3G45140 (ρ s = 0.98) loci, while in other cases the trend was not exact.Inherently, minimal, and imperceptible changes between the experimental biological replicas explain this phenomenon 29 .
Even though the magnitudes differed since the quantification methods have different principles 30 , we observed similar trends.

Differential expression analysis
To obtain a general overview of the biological processes involved in A. thaliana when in interaction with Trichoderma spp., we generated a heat map using differential expression contrasts and GO enrichment analysis for each time point (48, 72, and 96 hpi) (Fig. 2; Table S2).Figure 2 shows the biological processes that were differentially expressed in A. thaliana during its interaction with T. atroviride or T. virens, including mechanisms related to hormonal pathways such as camalexin metabolic process (GO:0052317), response to salicylic acid (GO:0009751), camalexin biosynthetic process (GO:0010120), and response to jasmonic acid (GO:0009753).In addition, defense mechanisms to other organisms including systemic acquired resistance (GO:0009627), response to molecule of bacterial origin (GO:0002237), response to oomycetes (GO:0002239), defense response to fungus (GO:0050832), immune system process (GO:0002376), defense response to nematode (GO:0002215), innate immune response (GO:0045087), and jasmonic acid and ethylene-dependent systemic resistance (GO:0009861).Finally, mechanisms that help alleviate reactive oxygen species were enriched (glutathione metabolic process (GO:0006749), response to hydrogen peroxide (GO:0042542).These results explain well the phenotype observed in the plants after their interaction with Trichoderma previously published by different groups 12,[31][32][33] .It has been reported that different species of Trichoderma affect the mechanisms of induced systemic resistance (ISR) in plants, which is related to the jasmonic acid/ethylene signaling pathways 12 .To contrast the induction of ISR by T. virens or T. atroviride, we visualized (Fig. S3) the differential expression of some relevant genes for this biological process.The LOX2 and PDF1.2 genes are ISR marker genes.Despite having no statistical significance (p-val = 0.1 & p-val = 0.2, respectively), when we performed a Wilcoxon test, we observed that T. virens was more efficient for the induction of marker genes (Fig. 3S).Besides, the WRKY51 (AT5G64810), a JA-regulated gene, showed a lower induction by T. virens than that by T. atroviride (Fig. S3).The transcription factor WRKY51 was reported by Yan et al. 32 as a negative regulator of the jasmonate pathway, which forms a protein complex with the JAV1 and JAZ8 proteins.This protein complex (JAV1-JAZ8-WRKY51) controls the biosynthesis of JA so that the plant can defend itself from insect attack 32 .Our results show that T. virens is more effective in inducing ISR than T. atroviride.Our results agree with those reported by Salas-Marina et al 15 , where they observed that T. virens had an enhanced effect on the induction of resistance against the necrotrophic fungus Botrytis cinerea in tomato plants (Solanum lycopersicum).
On the contrary, the ontological term development of the root (GO:0048364) was repressed in both cases after 96 hpi; these data agree with what was reported by Garnica-Vergara et al. 33 , where they observed an inhibition of primary root growth and induction of lateral root formation in A. thaliana after its treatment with T. atroviride.Garnica-Vergara mentions that the plant modifies the architecture of the root due to a blockage in auxin signaling, highlighting the participation of indole-3-acetic acid (IAA).Among the genes found in this group was AT5G56320, which codes for a plant expansin that causes the expansion of cell walls in plants and is www.nature.com/scientificreports/known to promote the growth of lateral roots.Additionally, we observed repression of the transcription factor LRL3 (AT5G58010) involved in root hair development.These facts show the complex interaction between plant hormones and their repercussions on the phenotype 33 .

Differential network analysis
The plant stress response is characterized by an intricate network of signaling cascades that trigger transcriptional reprogramming 34 .To analyze the behavior of genes in A. thaliana under control conditions and under stimulus conditions, we performed differential network analysis (Fig. 3A,B).The differential network analysis showed a change in the topology of the interaction networks (network (T.atroviride) and network (T.virens) ) against their controls.This fact represents the association behavior of genes under different conditions.Besides, the connections between the nodes were reconfigured, as seen in the histograms in Fig. 3A,B, where the control networks had more edges (33,222 and 116,682 for control-network (T.atroviride) and control-network (T.virens) respectively) than the interaction networks which had 9826 and 43,784 edges for network (T.atroviride) and network (T.virens) respectively.This fact suggests that the association between A. thaliana genes changes when Trichoderma colonizes the plant.The control network is denser and more connected with edge density = 0.06.In contrast, the T. atroviride and T. virens networks had edge densities of 0.008 and 0.007, respectively, when we filtered the correlation matrix with a p-val = 0.001.Even though the interaction networks have fewer edges than the control network they showed the highest correlation values higher than 0.9 and 0.95 for T. atroviride and T. virens, respectively.In contrast, the control network has correlation values higher than 0.8 (Fig. 3C).Taken results show that although interaction networks have fewer connections than the control, they are more robust and have higher correlations than their control networks.Thus suggesting a reconfiguration of the genetic transcriptional landscape of the plant when subjected to an external stimulus.

Gene co-expression network analysis
Here, we generated two gene co-expression networks of A. thaliana when stimulated with T. atroviride or T. virens, and these were simplified as explained below (Fig. 4S).From the transcriptome data, we generated the correlation matrices.The resulting networks consisted of 1094 and 2011 nodes for Arabidopsis (T.atroviride) and Arabidopsis (T.virens) , respectively (Fig. S4a,b).
An important step in network analysis is to know how the nodes group in communities.Newman and Girvan 35 defined a community as natural divisions of network nodes into densely connected subgroups.For this purpose, we used the algorithm developed by Pons & Latapi 36 .Our analysis identified 104 communities for the A. thaliana (T.atroviride) network and 265 for the A.thaliana (T.virens) network (Fig. S5).However, only communities Ca9, Ca16, and Ca18 for Arabidopsis (T.atroviride) and Cv5, Cv14, Cv22, and Cv41 for Arabidopsis (T.virens) were formed with a large number of nodes (Figs. 4, 5).
Due to the nature of the gene co-expression phenomenon, it is important that the genes form a compact and densely interconnected structure 37 .An adequate synergy between the functioning of the genes or their products (proteins) allows the cell to respond and adapt better to stress.Therefore, we focused only on those communities formed by at least 100 nodes (Figs. 4, 5).
Starting from the premise that communities within gene co-expression networks represent biologically functional modules, we generated a functional enrichment analysis for communities Ca9, Ca16, and Ca18 for A.thaliana (T.atroviride) and Cv5, Cv14, Cv22, and Cv41 for A.thaliana (T.virens) (Figs.4B, 5B; Table S2).According to graph theory, the most relevant nodes for a network are those that give the network the most stability 38 .To verify if our methodology to select relevant nodes was adequate, we used the networks formed by 1094 and 2011 nodes (Fig. S6).When we eliminated the relevant nodes, their general cohesion measures decreased with a diameter from 19.14 to 15.85 and 17.29 to 10.85, resulting in a much more compact and stable network.These results showed that our methodology to identify relevant nodes was adequate.
The GO results by community analysis (Figs. 4, 5), such as chorismate biosynthetic process (GO:0009423), defense response to fungus (GO:0050832) or aerobic electron transport chain (GO:0019646), agree with what has been reported.However, studies on the role of hypoxia during the plant-Trichoderma interaction are limited.In this sense, plants induce hypoxia-related genes during immersion, flooding, or when growing at high altitudes as the partial pressure of molecular oxygen decreases; however, little is known about the significance of hypoxia triggered by biotic stress.Thus, we focused on this biological process and its possible role during the Arabidopsis-Trichoderma interaction.
Despite the scarce studies on hypoxia during plant-microorganism interactions, there is experimental evidence indicating that hypoxia is a common phenomenon during the process.Establishing a symbiotic relationship with nitrogen-fixing rhizobia in legumes requires hypoxic conditions since nitrogenase activity, which can fix N 2 , is inactivated by free O 2 .In other plant species, hypoxia-sensitive genes are induced in stems infected by Agrobacterium tumefaciens, suggesting that hypoxic conditions are generated at sites of infection by these microorganisms 39,40 Although the plant-microorganism interaction largely depends on the lifestyle of the microbe, these studies show that hypoxia exists during the plant colonization event.However, to our knowledge, Trichoderma has not been reported to generate hypoxia in its hosts.Based on our results, we hypothesize that Trichoderma triggers hypoxia in the host and enhances defense mechanisms against pathogens by stabilizing a family VII transcription factor (ERF-VII) that responds to ethylene.
The association between communities Ca9 and Ca18 makes biological sense since it is known that plants obtain their energy during hypoxia from ATP alcoholic fermentation pathway 41 .On the other hand, the association between communities Cv14 and Cv22 in the A.thaliana (T.virens) network supports what was reported by 28 , who observed a functional relationship between the mitochondrial electron transport chain and hypoxia after evaluating the functioning of mitochondrial oxidases in different oxygen concentrations.Besides, the network www.nature.com/scientificreports/and functional enrichment analyses suggest that hypoxia may play a relevant role in the defense response against fungi (Figs. 4, 5).
To retrieve information, we searched for all the genes with an ontological term related to hypoxia within the general network (1094 and 2011 nodes; Fig. S4).We identified the AT2G47520 locus in both networks (Arabidopsis (T.atroviride) and Arabidopsis (T.virens) ), which codes for a transcription factor necessary to cope with hypoxia in A. thaliana 28 .Interestingly, this transcription factor has an Ethylene-inducible DNA binding domain (AP2/ERF position 49-106 a.a) (Fig. S7a).In addition, it has a cysteine (Cys2) prone to oxidation in its N-terminal (Fig. S7a,c).Therefore, this transcription factor is more stable during hypoxia than under normoxic conditions.It is known that the AP2/ERF domain can interact with the GCC-box promoter regions upstream of genes encoding pathogenesis-related proteins 42 .Therefore, the presence of these two domains (AP2/ERF and Cys2) in AT2G47520 suggests a dual role, first being more stable in hypoxic conditions to activate defense genes further.We analyzed positional weight matrices (PWM) and identified that the transcription factor can bind to the GCCbox (Fig. S7b,d).Brown et al. 43 , observed that deletion or point mutations in the core sequence of the GCC box substantially reduced the plant's ability to respond to JA.On the other hand, Epple & Bohlmann 44 observed that plants that overexpress the Thi2.1 gene, a gene with the GCC-box sequence in its promoter, improve resistance against the phytopathogen Fusarium oxysporum.
The AT2G47520 locus interacted with 22 genes in the Arabidopsis (T.atroviride) network.AT3G03270, AT4G24110, and AT3G27220 respond to hypoxia.AT1G43800 encodes a protease that participates in ISR, among others.In the Arabidopsis (T.virens) network, it was directly connected with 35 genes.Within this set of genes is AT2G26020, which encodes the plant defensin PDF1.2b, which is involved in defense mechanisms against fungi; AT5G07580, which codes for an ethylene response transcription factor; and AT1G66350, a protein with a role in jasmonic acid signaling.This fact suggests a functional relationship between hypoxia and the plant's defense mechanisms, as well as the participation of hormones in this genetic association.
The idea that plants induce their hypoxic genes and defense mechanisms against pathogens is consistent if we consider that plants suffer constant biotic stress after floods.On the one hand, there is dispersion of pathogens by rain, in addition to the fact that plants exhibit morphological changes after floods that facilitate infection by pathogens, such as change in intracellular pH, change in stomata behavior, and a decrease in hydraulic conductance 45 .These observations suggest that despite being a stressful stage for the plant, hypoxia promotes resistance against certain pathogens.For these reasons, we hypothesize that hypoxia and plant defense mechanisms are interconnected through signaling pathways that have not yet been explored and could involve signaling pathways related to the synthesis of JA.Yuan et al. 46 reported that after hypoxia, there is a rapid accumulation of (JA) and an induction of genes involved in the synthesis of JA.In addition, applying exogenous methyl jasmonate (MeJA) improved tolerance to posthypoxic stress in A. thaliana plants.A recent work 47 mentions that these two stresses occur sequentially or simultaneously in plants.During evolution, plants developed efficient strategies to respond to stimuli with the minimum energy expenditure through multiple functions of a single gene.

Experimental validation
To validate our transcriptomic analysis, we selected 12 genes, and we collected samples at times 48, 72, and 96 hpi.In this sense,10 genes were chosen from the main network, selecting those that will help us to solve future investigations and hypotheses and that represent different communities of the general network (1094 nodes).Additionally, 2 extra genes PR-1a (AT2G14610) and LOX2 (AT3G45140) were selected, which were not in the network but are commonly used as marker genes to validate the induction of systemic responses against pathogens.We compared the similarities between the differential expression results obtained by transcriptome and RT-qPCR (Fig. 6).During this analysis we observed a similar trend between the expression profiles obtained in the transcriptome, RNA-seq, and those obtained by RT-qPCR.In some cases, the exact same trend is observed for the transcriptomic data and the data obtained by RT-qPCR, as in the cases of the AT5G44420 (ρ s = 0.98), AT2G47520 (ρ s = 0.95), AT3G45140 (ρ s = 0.98) loci, while in other cases the trend is not exact (for example, AT2G16060 ρ s = 0.58, AT4G34200 with ρ s = 0.48 and AT1G05680 ρ s = 0.77).This fact is completely natural and expected, inherently there are minimal and imperceptible changes between the experimental biological replicas that explain this phenomenon 29 .

Conclusions
Our current study shows gene co-expression network analysis as a methodology that enables us to observe how genes associate into biologically functional modules in Arabidopsis that respond to an interaction with mutualistic fungi.
During plant-fungus interaction, there is a functional interconnection between hypoxic mechanisms and the plant defense mechanisms.However, we have carefully mentioned that we have only generated hypotheses requiring further experimental evidence.In the future, we will experimentally study the possible interconnection between hypoxia and defense mechanisms during the plant-fungus interaction.These interconnections will be identified in mutant plants in selected genes of our analysis.These experiments will allow us to clarify this complex mechanism that has led plants to respond appropriately to different types of stress by expending the minimum energy and anticipating future challenges that could end in plant death.
The AT2G47520 locus encodes a transcription factor with protein domains (AP2/ERF position 49-106 a.a. and Cys2) with dual function to counteract hypoxia and pathogens.

Quality analysis and mapping of the reads to the reference genome
Langebio provided a quality analysis of the FASTQ files.All reads presented a Phred quality score higher than 30, which was considered as high quality and used for mapping.Reads were mapped to a concatenated reference genome (A.thaliana-T.atroviride-T.virens).The concatenated genome was constructed taking the reference genomes reported by TAIR (The Arabidopsis Information Resource) version TAIR10 and for T. atroviride and T. virens reported by Ensembl Fungi (https:// fungi.ensem bl.org) version IMI_206040_v2.0 and ASM17099_v1 respectively.Mapping of the fragments was carried out using the HISAT2 software with default parameters for paired-end reads and a bash script 51 .Mapped fragments were quantified with the featureCounts tool from the SubRead package 52 .With these quantifications, a counts matrix was generated with an in-house R script.

Differential expression analysis
The differential gene expression analysis was performed with the edgeR package 53 using the count matrix.Normalization factors and variability were calculated with the calcNormFactors(x) and estimateDisp(x) functions, respectively.Subsequently, the data were fitted with the function glmFit(x).Contrasts were performed with a likelihood ratio test, glmQLFTest(x) function, to identify changes in gene expression under control and plantmicroorganism interaction conditions for each time (48, 72, and 96 hpi), with which a differential expression matrix was generated for each contrast.We keep differential expressed genes with a cutoff of FDR < 0.001 and a logFC > 1.5 or logFC < -1.5.

Differential network analysis
Two gene co-expression networks were generated (Arabidopsis (T.virens) and Arabidopsis (T.atroviride) ).All differential expression contrasts (Arabidopsis (T.virens) and Arabidopsis (T.atroviride) ) results were used to keep differentially expressed genes, using a false discovery rate (FDR) of FDR = 0.001.The control networks were built using the same genes that were present in the interaction networks (Arabidopsis (T.atroviride) ) 1094 and Arabidopsis (T.virens) , 2011 nodes) and the control libraries were used (A.thaliana without interaction).The control networks were filtered with an FDR = 0.001.A correlation test was performed with the cor.test(x) functions from the R base package.The correlations were filtered using a threshold of FDR = 0.001, and the correlation significance threshold was ρ s = 0.8.
Figure 1.Differential expression analysis.(A) MA-plot of the differential expression analysis A.thaliana (T.atroviride) .On y-axis is show the log 2 fold changes versus x-axis represents log 10 the average expression signal.Red dots represent up-regulated genes; blue dots represent down-regulated genes.(B) MA-plot of the differential expression analysis A.thaliana (T.virens) .On y-axis is show the log 2 fold changes versus x-axis represents log 10 the average expression signal Red dots represent up-regulated genes; blue dots represent down-regulated genes.(C) Comparison of genes up-regulated in A. thaliana by T. atroviride vs T. virens.Left Venn diagram showing up regulated genes by T. atroviride, T. virens or both.Middle.Comparison of frequency distribution of up-regulated genes by T. atroviride and T. virens Right.Results T-test of upregulated genes.The yellow dots represent the genes induced by T. atroviride while the pink dots were induced by T. virens.(D) Comparison of genes down-regulated in A. thaliana by T. atroviride vs T. virens.Left.Venn diagram showing down regulated gens by T. atroviride, T. virens or both.Middle.Comparison of frequency distribution of down regulated genes by T. atroviride and T. virens Right.Results T-test of downregulated genes.The yellow dots represent the genes repressed by T. atroviride while the pink dots were repressed by T. virens. https://doi.org/10.1038/s41598-023-48332-w

Figure 2 .
Figure 2.Gene ontology enrichment analysis.Biological process that were differentially regulated in Arabidopsis thaliana during its interaction with Trichoderma spp.Each column represents a contrast of differential expression, and each row represents the gene ontologies in terms of the biological processes involved.The color gradient goes from blue to red for downregulated and upregulated biological processes, respectively.

Figure 3 .
Figure 3. Differential network analysis.(A) Co-expression networks under control conditions (green edges), and in interaction with T.atroviride (yellow edges) and bottom histogram of degree.(B) Co-expression networks under control conditions (green edges), and in interaction with T. virens (pink edges) and bottom histogram of degree.(C) Correlation distribution analysis for the control, T. atroviride and T. virens networks.The color gradient goes from black to light yellow for low values and high values of correlation respectively.The graph presents the absolute value of correlation.

Figure 4 .
Figure 4. (A) Communities of co-expression genes network A.thaliana (T.atroviride) .Communities were identified using the cluster_walktrap(x) algorithm of the igraph package.Each community is represented by a color.(B) Functional enrichment analysis.The dark-brown color represents the community 9 (Ca9), the brown color represents the community 16 (Ca16) and the color light-brown represents the community 18 (Ca18).On the y-axis the enriched gene ontologies are shown on the x axis the fold enrichment is observed.

Figure 5 .
Figure 5. (A) Communities of co-expression genes network A.thaliana (T.virens) .Communities were identified using the cluster_walktrap(x) algorithm of the igraph package.(B) Functional enrichment analysis.The dark-brown and light-brown bars represent the Cv5 and Cv14 respectively.The dark-pink and light-pink bars represent the Cv22 and Cv41 respectively.On the y-axis the enriched gene ontologies are shown on the x-axis the fold enrichment is observed.

Figure 6 .
Figure 6.Transcriptome validation.Experimental validation of the transcriptome was performed by RT-qPCR at 48, 72 and 96 hpi.The validation experiment was performed by duplicate.Blue color values obtained of differential expression from the transcriptome; red color values obtained by RT-qPCR.The bars represent the standard deviation (sd) of the average of a triplicate.The correlation coefficient of Sperman (ρ s ). https://doi.org/10.1038/s41598-023-48332-w