Metatranscriptomic approach for microbiome characterization and host gene expression evaluation for “Hoja de malvón” disease in Vitis vinifera cv. Malbec CURRENT POSTED

Background: A particular GTD of enormous concern in Argentinean vitiviniculture is the complex etiology disease locally known as "Hoja de malvón" (HDM). At least four different fungi are involved in the disease, which complicates the diagnostic and the design of strategies for vineyard management. Similar to Esca grape disease, factors that make this disease difficult to control are the presence of pathogens that not always correlates with disease symptoms or physiological changes in the host. Also the abiotic stress on grapevine seems to favors the disease process. Based on this background, it is essential to have molecular tools that allow for simultaneous explorations of the host immunity status and the microbiome composition. Results: A metatranscriptomic approach was followed and different strategies for microbiome characterization, throughmolecular marker reconstruction or unique kmer counts, were evaluated. Malbec microbiome was mainly represented for Dothideomycetes and Actinobacteria. Higher Basidiomycota/Ascomycota ratio was found in symptomatic (SYM) than in asymptomatic (ASYM) plants, with the Basidiomycota Arrambaria destruens found in higher levels in SYM plants. Besides, from mRNA-derived reads, the host functional status and main microbial functions based on gene expression were evaluated. Stress-tolerance mechanisms were activated on ASYM plants, being spermidine synthesis one of the most important, while from the microbiome side, prokaryotic tricarboxylic acid pathways were the primary differential function in SYM plants. Conclusions: This is a pioneering work for the characterization of multi-kingdom endophytic microbiome in woody tissues of grapevine cv. Malbec. Several microorganisms with negative interaction with GTD pathogens were identified and can be further explored as biological control agents for HDM disease. The integral analysis employed suggests that measuring Basidiomycota/Ascomycota ratio and spermidine-associated gene expression would help to monitor the sanitary status of grapevine and the propensity to develop HDM disease. mapped on the grapevine genome and used for gene expression analysis. The unmapped reads (average of 42.7%; TableS2.xlsx) were used for microbial identification and quantification through a phylogenetic marker reconstruction approach or identifying and counting unique kmers.

Argentina's vitiviniculture has a big concern because of GTD, mainly with the Esca-like disease known as "Hoja de malvón" (HDM) [2]. The disease begins in one of the arms of the plant, affecting branches, leaves, and clusters. Leaves have a smaller size, become chlorotic, and with the inscribed edges curled down, taking the appearance of a "geranium" leaf, from which the disease takes the name. The shoots have less growth and shorter internodes than those in healthy plants. The clusters are smaller, lax, and without uniform grains. In the transverse section of the trunk, necrosis of different colorations and consistencies are observed that indicate the progression of the disease [2].
In the trunk of some plants, the presence of light brown colored symptoms is observed, corresponding to basidiocarps of the fungus Arambarria destruens (formerly Inocutis jamaicensis) [3]. In Chile, similar symptoms are observed for a disease named "Enrollamiento clorótico de la vid" [4]. Around the world in Esca disease, Phaeoacremonium spp. and P. chlamydospora are always isolated and are associated with different Basidiomycota-Hymenomycete fungi. In Europe and the Mediterranean grape-growing regions, Fomitiporia mediterranea is the prevalent species; in Chile Fomitiporella vitis and, in Australia, Fomitiporia australiensis [5]. In Argentina, the Basidiomycota was identified as Inonotus sensu lato (Inocutis Fiasson & Niemelä) [3], later reclassified as Inocutis jamaicensis [6] and recently, recognized as a new species named Arambarria destruens [7]. P. aleophilum and P. chlamydospora are frequently isolated from symptomatic tissues, but in lower frequency than the Hymenomycete A. destruens. HDM is commonly found in Argentinean vineyards in more than 20 years old plants, which correspond to 46.8% of the total vineyards of the country [8]. In Argentina, old vineyards are maintained because of the wine quality obtained from the plants and for the economic limitations to replant vineyards. This particularity obliges to reinforce disease prevention through periodic monitoring of the vineyard, reducing the opportunity to get infected plants with GTD.
The monitoring of GTD pathogens on vines is mainly done using conventional cultural isolations coupled to morphological characterization. Few molecular techniques for pathogen monitoring have been developed [9][10][11], but little has been done regarding HDM disease. Molecular tools are needed to understand the importance of each pathogen involved in HDM disease that correlates with the presence of symptoms and to understand which physiological change makes grapevine more susceptible. A non-destructive and easy to adopt technique will facilitate extensive analyses of the pathogens present in vineyards. Molecular techniques based on high throughput sequencing (HTS) are more adaptable than conventional cultural-based techniques to reach such aim.
Different techniques have been used to characterize the microbiome associated with Vitis vinifera.
The metatranscriptomics approach allows for the direct inspection of the whole sanitary status of a plant because of the simultaneous detection of several pathogens like viruses, bacteria, fungi, and archaea in a kingdom-independent manner. Also, it is possible to characterize the metabolically active microorganisms, avoiding the dead ones considered in a metagenomics approach. This approach could help to understand the interaction established among the microorganisms with their host, from a functional point of view.
Based on previous reports on GTD, abiotic stress predispose plants to fungal pathogen colonization, evidencing the importance of having a simultaneous exploration of the pathogen and the status of the host. Although several endophytic microbiome characterizations of different grapevine cultivars have been done worldwide on the main wine-producing regions [14], to the best of our knowledge, there is no previous exploration of the microbiome in Vitis vinifera cv. Malbec cultivated in the Mendoza region. Thus, this work aimed to identify differences in terms of microorganisms and plant gene expression that could help to identify molecular markers for monitoring and the early detection of HDM in vineyards. For that, a simultaneous exploration of both microbiome and host gene expression in symptomatic and asymptomatic plants of Vitis vinifera cv. Malbec was performed.

Results
Grapevine endophytic microbial composition determined through culture isolation and morphological characteristics   was not adequate for fungal classification in low taxonomic levels (genus and species) because most of the taxonomic names assigned to the sequences deposited there are out of date, and most of the associated GTD pathogens are absent (which is the focus of this research). Therefore, fungal community characterization shown here is based on MATAM classification using Warcup-RDP-ITS, an alternative updated reference for fungi [28]. Ascomycota was almost the unique represented phylum (97-99% of all the identified fungi. Figure 1 and FigureS3.html). Basidiomycota represented only around 0-0.3% in ASYM plants and 0.4-0.9% in SYM plants ( Fig. 1 and Figure S3.html). From the Ascomycota, the Dothideomycetes were the most abundant (54-72% of all the identified fungi) at the class taxonomic level, followed by Sordariomycetes (7-19%), Eurotiomycetes (7-15%), Leotiomycetes
Compared with relative abundance determined through MATAM, the unique kmer counts coincided with the determination of Actinobacteria as the most abundant Phylum. On the other hand, Bacteroidetes were the second most abundant based on MATAM, but based on KrakenUniq output, they were recognized as relatively abundant (7%) but bellow of Firmicutes (11%) and Proteobacteria (42-43%). In MATAM output, Firmicutes were not recognized at all.
At the genus level, the differences with MATAM were even more evident; for instance, the relative abundance of Segetibacter sp. was around 0.03-0.04%, according to KrakenUniq, quite different from 10-18% obtained with MATAM. Similar differences were found for Pseudonocardia sp.  (Table I). The comparison of methods used for microbiome determination showed significant differences. For KrakenUniq, in general, the lower percentage for microbiome members was accompanied by higher diversity in terms of richness. However, there are some exceptions, i.e. Paraconiothyrium sp. and Harzia sp. that were quantified in higher relative abundance in KrakenUniq output than in MATAM output (Table I).

Interaction Between Microorganisms Based On The Correlation Matrix
Without taking into consideration the origin of the sample, a matrix of correlations was done for the fungal community to identify putative biological control agents (BCAs) interacting negatively with GTD pathogens. Focusing on the pathogens A. destruens and P. minimum, the primary negatively interacting fungi were found in the order Pleosporales, which includes Alternaria alternata, Alternaria solani, Parastagonospora nodorum and Didymella pinodes (Fig. 3). Also, it was clear that there are mainly inverse correlations between Basidiomycota and Ascomycota fungi ( Fig. 3  despite having a sparse microbial composition were more similar between them than sample m16, which shows a clear enrichment in the Basidiomycota phylum ( Fig. 2A). The symptoms from m16 sampled plants were more pronounced, as a branch was severely affected than those obtained from m11 and m26, indicating that the increase in Basidiomycota corresponds to the decay of vine tissues, typical behavior of fungi involved in white rot. This fact was taken into consideration for further statistical analyses, and only the unique kmers counts obtained from c23 + c25 in ASYM plants and m11 + m26 in SYM were used. The comparison showed that three genera of Archaea have significant relative abundance differences (FDR < 0.15). Halococcus sp., Haloquadratum sp., and Halalkalicoccus sp. presented higher relative abundance in ASYM plants ( Fig. 4 and TableS4.xlsx). In contrast, A.
destruens showed a higher relative abundance in SYM plants ( Fig. 4

Grapevine Gene Expression And Functional Enrichment Analysis
Using counts of reads mapped in the Pinot Noir PN40024 grapevine reference genome, the differential gene expression was calculated. Fifty-seven up-regulated genes, and 159 down-regulated genes were identified (Using a threshold of FDR < 0.05). The grapevine Gene Ontology (GO) term enrichment analysis in SYM vs. ASYM plants, GO:0008295 'spermidine biosynthesis process' was enriched among downregulated genes when considering FDR < 0.05 threshold, and 34 GO terms were enriched when considering a p-value < 0.005 as a threshold ( Fig. 6  The approach used here, based on the sequencing of ribodepleted RNA, was chosen to enrich the target host mRNA molecules and the rRNA from a microbial community. To the best of our knowledge, this is the first time this strategy is used for the study of the GTDs-grapevine pathosystem. One of the Interestingly, Phaeoacremonium spp. and P. chlamydospora, regularly isolated from HDM plants [2,6], were not identified in differential abundance, even by high sensitive molecular methods. Despite not differential, Phaeoacremonium minimum (formerly P. aleophilum) was present at low but detectable levels (> 0.0001 relative abundance) in all the samples,. This result, is in agreement with a previous report of P. minimum as a common isolate from SYM HDM plants [2]. P. chlamydospora was not detected by cultural isolation nor by molecular methods. P. chlamydospora has been mainly associated with young plants less than eight years old [30]. Since the plants sampled were around 23 years old, this could explain why P. chlamydospora was not detected at all. In any case, more work needs to be done to understand the role of P. chlamydospora in HDM disease.
For the prediction of antagonistic microorganisms to HDM pathogens, first, we suggested those with a negative correlation based on fungal presence/absence (Fig. 3) [34]. Our data showed several Trichoderma species but in low proportions. Thrichoderma reesei was more represented in ASYM plants c23 (> 0,001), suggesting it could be a good candidate for further biological control ability evaluation. Another Ascomycota previously found as BCA candidate for the GTD pathogen Diplodia seriata is Aureobasidium pullulans [35]. In the present work, A. pullulans was less represented (< 0.001), but A. namibiae and A. subglaciale were found in levels above 0.001, therefore could be further evaluated for the capacity to control GTD.
On the other hand, the prokaryotic microbiome indicates that Actinobacterias were differentially represented in symptomatic vines. Bruez et al. [22] and Elena et al. [36] explored the bacterial microbiome in Esca, but did not identify any BCA candidates. A role of Actinobacterias as BCA has been reported [37], but contrasting to our results, Streptomyces spp. (actinobacterias) were found antagonist to GTD pathogens, while our data showed they were mainly represented in SYM plants; therefore, a promoting disease role was assumed. Recently, it was proposed that the Actinobacteria/Bacteroidetes ratio in rhizosphere increases in domesticated plants with low tolerance to biotic and abiotic stress [38], suggesting that this ratio can help for monitoring the sanitary status of grapevine. Flavobacterium commune was the only bacteria identified with a significantly higher abundance in ASYM plants (Fig. 4). Flavobacterium sp. was reported to inhibit the germination of the GTD pathogen Lasiodiplodia theobromae [39]; therefore, a role as a BCA is hypothesized.
Grapevine gene expression and functional enrichment indicate that processes related to stress, as responses to unfolded protein, heat, drought, and nutrient deficiencies, are induced in symptomatic grapevines. This global response indicates that SYM plants are under stress, but it is not possible to determine if such stress is the result of the management of the vineyard or is due to a change in a microbiome that favors the colonization of A. destruens. Despite the triggering cause, it is clear that stress has a significant role in the appearance of symptoms. This was previously addressed [40], however to the best of our knowledge, it was never found using an integrative "in field" approach.
Previous works evaluated the effect of stress on GTDs pathosystem using only one or two GTDs pathogens [41][42][43]; thus, this is the first time that all the complexity of agents involved in HDM disease in Vitis vinifera cv. Malbec is addressed in Argentinean field conditions. Results highlight the importance of avoiding stress on plants of local vineyards.
In Botryosphaeriaceae, abiotic stresses such as drought and heat were suggested as triggering factors for the endophyte-pathogen transition [44][45]. A recent review illustrates the hypothesis for GTD occurrence [40] based on those with direct effect (abiotic stress on the plant, pathogen, or both) and those with indirect effect (abiotic stress perturbing microbiome interactions). In this study we explored the microbiome and grapevine gene expression simultaneously. Polyamines are required for plants to deal with abiotic and biotic stress; therefore, it is difficult to understand why grapevine is repressing such functions in symptomatic plants. Also, microorganismscould synthesize polyamines.
Therefore a pathogen hijacking strategy must not be discarded [46]. Probably in ASYM plants, polyamine synthesis is the main defense that restricts pathogen colonization.
From the microorganism's functional side, the TCA and TCA-glyoxylate bypass cycles are the main pathways activated on the microbial population hosted in SYM plants. The glyoxylate bypass allows the cell to synthesize glucose uniquely from acetate or dicarboxylic acid derived from fatty acid degradation [47]. This strategy could be fundamental to grow in a tissue with a considerable limitation of carbohydrates, such as grapevine vascular tissue, especially under stressful conditions because the plant redirects sugars toward the defense metabolic pathway. On the other side, a recent study indicates that microbes involved in hydrocarbon degradation are under oxidative stress, because of the sub-production of hydrogen peroxide during the oxidation of such compounds. In response to that environment, microbes upregulate glyoxylate shunt pathways to keep gluconeogenesis and growth working under such stressful conditions [47]. The degradation of aromatic compounds is required to deal with several components of grapevine woods such as pectin and lignin; therefore, it is quite probable that microorganisms are obligated to activate that pathway for survival.
Besides studying the host gene expression, we aimed to identify markers for the early detection of HDM disease. Due to the association of the disease with stress, it is essential to simultaneously quantify general stress marker expression and specific microorganism abundance, as a strategy to prevent HDM disease. Such markers will be beneficial, since periodic monitoring of vineyards along with preventative managing practices, would extend the production time of grapevines. We believe that one of the main findings of this study is the increased ratio of Basidiomycota/Ascomycota in SYM plants and the significantly increased polyamine biosynthesis in ASYM plants. Therefore specific quantification of such molecular markers must be developed and validated to be used as an indicator of early HDM stages or to detect levels of susceptibility to develop this disease.
Since only two plants were used for comparison of host gene expression and microbiome, conclusions must be taken with caution. Of course, more plants need to be analyzed to have higher confidence.
However, the originality of this work has sustenance in the determination of microbiome from RNA derived reads allowing a higher resolution in terms of taxonomical "deepness" of classification than classical metabarcoding analyses.
We hope the approach here presented serves as a road-map to develop a more specific evaluation of the microbiome-pathogens-grapevine interactions in GTD.

Conclusions
The integral analysis of the multi-kingdommicrobiome with grapevine gene expression indicates that stress-responsive genes (through polyamine and trehalose synthesis) are the main factor explaining differences between symptomatic and asymptomatic plants. These findings support the previous suggestions to manage GTD using preventive practices that guaranteed the minimization of stress in the vines.

Sample collection
Plants Plates were incubated at 24 °C, monitored daily, and after one week, evaluated for morphological characterization. Observed bacterial or fungal growth was replicated onto fresh MEA medium and incubated at 24 °C again. Cultures grown in MEA were initially classified based on the macroscopic mycelium characteristics, including colony shape, texture, color, and microscopic features, including shape and color of mycelia, and shape of the conidia. The identified microorganisms were counted per sample and per isolates, considering the pure cultures to estimate the abundance.

Rna Extraction And Sequencing
Three SYM and three ASYM plants were selected for RNAsEq. A pool of grapevine woodchip was obtained and ground with a mortar and liquid nitrogen. Around 400 mg of the resulting powder was used for RNA extraction, adapted from Chang et al., [48] in 15 mL centrifuge tubes. Briefly, the grapevine powder was incubated with cetyl trimethyl ammonium bromide (CTAB) extraction buffer after incubation at 65 °C for 10 min, a Chloroform:Isopropyl alcohol mixture (24:1) was added and after centrifugation, the aqueous phase was recovered. Then a precipitation step using 3M pH 8 sodium acetate was done, and the RNA was precipitated with 10M lithium chloride. Finally, two washing steps with 75% ethanol were included. Total RNA obtained was treated with RQ1 Promega DNAse to eliminate contaminant DNA, and RNA cleaning was done with R1015 Zymo Clean and Concentration kit following the manufacturer's recommendations.
Total RNA was quantified in a Nanodrop2000, and its integrity was determined in 1.5% agarose gel. The mapped read counts matrix was used for edgeR v3.20.9 analysis. Based on the similarity of microbial composition samples c23 and c25 were used as "controls" while samples m11 and m26 were used as "symptomatic". When only those samples were used, several differentially expressed genes (degs) were found. GO terms enrichment was evaluated for degs list through Goseq v1.32.0 [53]; the output of enrichment was used to make plots through Goplot v1.0.2 [54].

Microbiome Analysis
Because of the recent development of bioinformatics tools to process metatranscriptomic data, there are no clear delineations or consensus regarding what is the best tool for characterization of the microbiome composition from metatranscriptomic derived data. In this work, three alternatives were carried out. The first was using MATAM v1.5.3 [60] to identify all the rRNA derived reads (except those from plants that were previously removed) and reconstruct the phylogenetic marker gene through a reference-guided assembly strategy. The second was through the use of KrakenUniq v0.5.7 [61] that takes all the trimmed and filtered reads for k-mer indexing as input and then map the k-mer in nucleotide microbial DB and use the unique k-mer counts for relative abundance estimations. The last was through the use of Humann2 v0.11.2 [62] that use marker DB as a reference, and later construct pangenomes reference only with the microorganisms previously identified for a more in-depth taxonomic characterization in a more manageable microbial part of 'nt' NCBI DB.

Matam Assembly With Prokaryotic-16s-silva And Fungal Its-rdp Model
The unmapped reads on Vitis vinifera cv. Pinot Noir genome filtered to remove remaining plant rRNA and organelle derived reads, as explained before, were used for mapping and reconstruction of ribosomal regions and comparison on SILVA-SSU-DB [63] for prokaryotic, SILVA-LSU-DB [63] for eukaryotic, and ITS-RDP [29] for fungal classification, through MATAM v1.5.3 [60]. In the case of the Warcup_RDP_ITS DB [28], the threshold for clustering reference sequences was set to 200 bp when running matam_db_preprocessing.py script, and when running the assembly, the minimal scaffold length was set to 200 bp.The training model for the assembly was specified as "fungalits_warcup".
Unique k-mers used for the whole classification and relative abundance quantification KrakenUniq v0.5.7 [60] was used for splitting the grapevine genome unmapped reads in kmers of 31 bp length and directly mapping on the microbial part of NCBI 'nt' DB [29].  [62] and then used as a reference for mapping in the second step, while the unmapped reads in the first step were processed through diamond against UniRef90 DB [67]. At the end, the counts of reads were grouped by gene families (GF) or metabolic pathways (through the relationship GF-MetaCyc DB [68]

Consent for publication
Not applicable.

Availability of data and materials
The dataset supporting the conclusions of this article is available in the ENA-EMBL repository, [PRJEB31098 https://www.ebi.ac.uk/ena/data/view/PRJEB31098]

Competing interest
The authors declare that they have no competing interests

Funding
This research was supported by Instituto Nacional de Tecnología Agropecuaria through the 'Proyecto  Interactions between members of the fungal community of the grapevine microbiome.
Positive (trends to blue) and negative (trends to red) correlations among highly abundant fungi identified through unique kmers using KrakenUniq software. The "Hoja de malvón" pathogen Arambarria destruens have negative interactions with three members of Pleosporales order: Alternaria alternata, Parastagonospora nodorum, and Didymella pinodes. Another pathogen Phaeoacremonium minimum has negative interactions with the same fungal community as negatively interacts with Arambarria destruens, Leptosphaeria biglobosa, and Alternaria solani. Microbial functional pathways in SYM vs. ASYM comparison. Using humann2 the main abundant microorganisms were identified through a multi-marker strategy, a pangenome (only having the identified microorganisms) was constructed to be used as a reference for mapping the raw reads of each sample. Counts of mapped reads on genes grouped through the MetaCyc pathways database were used to calculate significant differences through the ANOVA Tukey-Kramer test (p-value<0.05). A) Heatmap showing the main differential pathways. B), C), D), and E) corresponds to bar plots from TCA cycles-associated pathways.
The bypass of the glyoxylate cycle might indicate that microorganisms with great functional activities activate a process of anabolism from lipids and/or amino acids that could be available because of host cell disruption.