Transcriptomic Analysis of Quinoa Reveals a Group of Germin-Like Proteins Induced by Trichoderma

Symbiotic strains of fungi in the genus Trichoderma affect growth and pathogen resistance of many plant species, but the interaction is not known in molecular detail. Here we describe the transcriptomic response of two cultivars of the crop Chenopodium quinoa to axenic co-cultivation with Trichoderma harzianum BOL-12 and Trichoderma afroharzianum T22. The response of C. quinoa roots to BOL-12 and T22 in the early phases of interaction was studied by RNA sequencing and RT-qPCR verification. Interaction with the two fungal strains induced partially overlapping gene expression responses. Comparing the two plant genotypes, a broad spectrum of putative quinoa defense genes were found activated in the cultivar Kurmi but not in the Real cultivar. In cultivar Kurmi, relatively small effects were observed for classical pathogen response pathways but instead a C. quinoa-specific clade of germin-like genes were activated. Germin-like genes were found to be more rapidly induced in cultivar Kurmi as compared to Real. The same germin-like genes were found to also be upregulated systemically in the leaves. No strong correlation was observed between any of the known hormone-mediated defense response pathways and any of the quinoa-Trichoderma interactions. The differences in responses are relevant for the capabilities of applying Trichoderma agents for crop protection of different cultivars of C. quinoa.

Symbiotic strains of fungi in the genus Trichoderma affect growth and pathogen resistance of many plant species, but the interaction is not known in molecular detail. Here we describe the transcriptomic response of two cultivars of the crop Chenopodium quinoa to axenic co-cultivation with Trichoderma harzianum BOL-12 and Trichoderma afroharzianum T22. The response of C. quinoa roots to BOL-12 and T22 in the early phases of interaction was studied by RNA sequencing and RT-qPCR verification. Interaction with the two fungal strains induced partially overlapping gene expression responses. Comparing the two plant genotypes, a broad spectrum of putative quinoa defense genes were found activated in the cultivar Kurmi but not in the Real cultivar. In cultivar Kurmi, relatively small effects were observed for classical pathogen response pathways but instead a C. quinoa-specific clade of germin-like genes were activated. Germin-like genes were found to be more rapidly induced in cultivar Kurmi as compared to Real. The same germin-like genes were found to also be upregulated systemically in the leaves. No strong correlation was observed between any of the known hormone-mediated defense response pathways and any of the quinoa-Trichoderma interactions. The differences in responses are relevant for the capabilities of applying Trichoderma agents for crop protection of different cultivars of C. quinoa.
Keywords: Chenopodium quinoa, transcriptomics, germin-like proteins, Trichoderma, plant-fungi interactions BACKGROUND Trichoderma is a genus of ascomycete fungi widely studied for its versatile interactions with other organisms. Trichoderma can feed or parasitize on other fungi, bacteria, oomycetes and nematodes (Harman et al., 2004;Verma et al., 2007). Several species of Trichoderma are also symbionts with plants and can promote plant growth by several, yet so far only partially known mechanisms. Strains of several symbiotic species in the Trichoderma harzianum species complex (e.g., Trichoderma afroharzianum T22 (Chaverri et al., 2015), previously called T. harzianum) are used commercially because they can substantially improve yields of several species of crops (Harman et al., 1989;Monte, 2001;Hermosa et al., 2012). The strain T22 can improve the soil nutrient availability to plants (Altomare et al., 1999), and several species of Trichoderma have also been shown to enhance plant growth through volatile compound emission (Lee et al., 2016;Jalali et al., 2017) and stimulate plant systemic defense responses (Harman et al., 2004;Mukherjee et al., 2013). Nevertheless, plants do not always benefit from these interactions as described for some maize cultivar in field trials and lab experiments (Harman, 2006). Plant growth inhibition by T. harzianum has also been observed in quinoa seedlings after six but not two days of axenic co-culture (Rollano-Peñaloza et al., 2018).
Quinoa (Chenopodium quinoa Willd.) is an emerging crop of great interest due to its nutritional values (Vega-Gálvez et al., 2010) and its resistance to hostile environmental conditions, especially salinity and drought (Ruiz et al., 2014;Bazile et al., 2016). Quinoa seeds are gluten-free, contain all essential amino acids and its composition (vitamins, antioxidants, fatty acids and minerals) is highly suitable for human nutrition (Repo-Carrasco et al., 2003). Quinoa has a high genetic diversity, e.g., >4,000 accessions have been registered by the Food and Agriculture Organization (Zurita-Silva et al., 2014). The high genetic diversity of cultivars is the result of many years of selection by the indigenous people of the Andean Altiplano, where quinoa may have been domesticated 7,000 years ago by pre-Columbian cultures (Dillehay et al., 2007).
Quinoa agricultural yields can be boosted by Trichoderma application, as previously described (Ortuño et al., 2013). However, the outcome of plant-Trichoderma interactions is not always beneficial. Plant genotype-specific growth inhibition by commercially available Trichoderma strains have been reported for lentils (Prashar and Vandenberg, 2017), tomato (Tucci et al., 2011) sugarbeet (Schmidt et al., 2020) and maize (Harman, 2006). Thus, the incompatibility of particular plants with particular biocontrol strains can lead to undesired agricultural losses. Therefore, there is a need to understand the genotype-specific mechanisms that determine beneficial plant growth effects upon treatment with biocontrol agents like T22. Among the beneficial plant growth mechanisms triggered by biocontrol agents and mycorrhiza, the activation of plant defenses is one of the most studied topics. Plant defenses activated by mycorrhiza have been shown to include expression of Germinlike proteins (GLPs) (Fiorilli et al., 2018), and plant-fungi protein-protein interaction studies have shown that GLPs interact with Trichoderma cellulases (Saravanakumar et al., 2018). However, no studies about effects of Trichoderma on GLPs expression have been reported.
GLPs are plant proteins that have been especially well studied in grains like barley. These proteins are characterized by being associated with several enzymatic activities including oxalate oxidase (OXO), superoxide dismutase (SOD), ADP-glucose pyrophosphatase/ phosphodiesterase (AGPPase) and polyphenol oxidase (PPO) (Zimmermann et al., 2006;Dunwell et al., 2008). A potential function of germin-like proteins (GLPs) is found in its OXO and SOD activities, which may play a key role in production of hydrogen peroxide (H 2 O 2 ) during plant defense (Ilyas et al., 2016). Because of a potential importance of GLP in protecting plant cells from superoxide toxicity produced under pathogen attacks, germin or GLP genes (i.e., HvOXO1) have been inserted into dicot plants like rapeseed or peanut to enhance their pathogen resistance (Livingstone et al., 2005;Dunwell et al., 2008;Liu et al., 2015).
In this work, we have studied the molecular response of two C. quinoa cultivars that have been shown to experience plant growth inhibition when treated over longer time with T. harzianum BOL-12 and T. afroharzianum T22 in axenic co-cultures. The response of quinoa to BOL-12 and T22 in the early phases of interaction was studied by transcriptomic analysis and RT-qPCR verification. Overall, we observed that upon interaction with the two fungal strains, a broad spectrum of putative quinoa defense genes were activated in Kurmi but not in the Real cultivar.

Fungal Growth
T22 and BOL-12 were maintained on potato dextrose agar (BD-Difco, Detroit, USA) at 25 • C. To isolate conidiospore suspensions, one ml of sterile water was added to two-weekold Trichoderma cultures on potato dextrose agar and collected spores were filtered through a sterile piece of cotton wool. The spores were washed twice with sterile H 2 O (Milli-Q, Merck Millipore, Burlington, MA, USA) and pelleted at 3700g for 5 min at 4 • C in an Allegra X-12R centrifuge (Beckman, Brea, CA, USA). Spores were resuspended in sterile H 2 O and kept at 4 • C until experiments.
Germination of T22 and BOL-12 spores for C. quinoa treatment was performed as described by Yedidia, Benhamou (Yedidia et al., 1999) using 15 ml tubes shaken at 200 rpm for 18 h. The germinated spore suspension was washed twice by centrifugation as described above and finally resuspended in sterile H 2 O. The final spore concentration was adjusted to be 1 germinated spore/µl and verified by colony forming unit (CFU) counts in potato dextrose agar Petri dishes.
supplemented with 8 g/L agar. The Petri dishes were tilted 45 • during growth with the agar/air interface facing upwards and seedlings having the roots pointing toward the bottom part of the Petri dish. The seedlings were incubated at 24 • C for 4 h before treatment with T22 or BOL-12.
C. quinoa seedlings were treated by adding 10 µl [1 CFU/µl] of either T22 or BOL-12 germinated spore suspension on the neck of the primary root. Ten µl of sterile H 2 O was added to each seedling in the mock control group. After treatment, the seedlings were incubated at 24 • C in a 16 h light /8 h dark photoperiod. Cocultivation was done under fluorescent lights (Polylux XLr 30W, GE, Budapest, Hungary) at 50 µmol m −2 s −1 for 12 and 36 h.

Seedling Growth Analysis
Hypocotyl length was analyzed from images taken with a Digital Camera Canon EOS Rebel T3. Measurements from the photographs were done with the segmented line tool of ImageJ 1.49 (Abramoff et al., 2004).

Sample Collection and RNA Extraction
For RNA extraction quinoa seedling were sampled 12 and 36 h after Trichoderma treatment. Each treatment enclosed five plates containing five seedlings. To reduce inter-plate variability, one root from each of the five plates was pooled for a biological replicate into pre-weighed aluminum foil envelopes. Roots were excised at the root-hypocotyl interface with a scalpel. The envelopes were weighed on a precision balance and shockfrozen in liquid nitrogen. Frozen samples were either processed immediately or stored at −80 • C until RNA extraction. Roots and shoots were pooled separately.
Total RNA was extracted using the RNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA), with the following modification: Root tissue samples preserved in liquid nitrogen were placed in a precooled mortar containing liquid nitrogen followed by thoroughly grinding without letting the samples thaw. Then 450 µL of Buffer RLT (Qiagen, Valencia, CA, USA) supplemented with B-mercaptoethanol (1%) was added. Grinding continued until samples thawed and were transferred to a 1.5 ml microcentrifuge tube. The rest of the procedure was followed according to Qiagen instructions. Total RNA quantity and quality was determined with a NanoDrop spectrophotometer. DNase treatment was performed with the DNA-free kit (Ambion, Carlsbad, CA, USA), following the instructions of the manufacturer. The integrity and quality of the RNA was determined as follows: 500 ng of DNase-treated RNA were dissolved in 8 µl of sterile water and split in two aliquots, one placed on ice and the other placed at 37 • . After 20 min incubation, 2 µl of loading buffer was added to each sample and both aliquots were loaded on an agarose gel (2%) stained with ethidium bromide. The gel was run at 80 V for 30 min and visualized in an UV-transilluminator. Samples with sharp 18S and 28S rRNA and showing no evidence of degradation were retained.

RNA-Seq Library Construction and Sequencing
Total RNA treated with DNase was sent to IGA technology services (IGA, Udine, Italy; http://www.igatechnology.com) for poly(A) + mRNA purification, strand-specific cDNA synthesis, library construction (Truseq stranded mRNA-seq) and sequencing using a HiSeq2500 (Illumina Inc., San Diego, CA, USA) in paired-end mode with a read length of 125 bp. Raw sequences have been deposited at the National Center of Biotechnology Information (NCBI) under project accession number: PRJNA720675.

Transcriptomic Analysis
RNA-seq reads were checked for quality by FastQC (v.0.9.0) and mapped on the quinoa genome "Kd" (Yasui et al., 2016a) and to the QQ74 coastal genome (Jarvis et al., 2017) by Tophat2 (v.2.2.9). Transcript abundances were assessed with HTSeq (v.0.9.1) with "intersection-nonempty" mode. Genes that had a minimum of 1 read mapped in each of the samples considered for analysis were included. Gene expression levels were measured as counts per million (CPM) (Anders et al., 2013). Library size normalization was performed using the trimmed mean of M-values (TMM) within the R package edgeR (v.3.14.0) (Robinson et al., 2010;Dillies et al., 2013). CPM were TMMnormalized in order to compensate for library size differences. Differential gene expression analysis comparing mock-treated samples with samples treated with Trichoderma was performed using edgeR with TMM normalized libraries (Anders et al., 2013) with a false discovery rate (FDR) of 5% (q < 0.05) (Benjamini and Hochberg, 1995).

Functional Annotation of Differentially Expressed Genes
Gene ontology (GO) term enrichment for sets of differentially expressed genes were estimated with Argot2 through sequence function prediction (Falda et al., 2012). Singular enrichment analysis (SEA) for biological processes was performed with AgriGO v2.0 (Du et al., 2010). The statistical test for SEA was Fisher's exact test and for false discovery rate the Yekutieli method was applied (Tian et al., 2017).

cDNA Synthesis and Gene Expression by qRT-PCR Analysis
Synthesis of cDNA was carried out with 500 ng of total RNA added to each 20 µl reaction of the RevertAid H Minus Reverse Transcription Kit (Thermo Scientific). The cDNA samples were stored at −20 • C for downstream analysis. qRT-PCR of plant RNA was performed in a CFX384 Touch Real-Time PCR system (Bio-Rad, Hercules, CA, United States) using Maxima SYBR Green qPCR Master Mix (Thermo Scientific) supplemented with 0.25 µM of each specific primer and 10 ng of cDNA as template in a total volume of 10 µl/reaction. The PCR program had the following conditions: 1 cycle of: 95 • C, 20 s; 30 cycles of: (95 • C, 15 s; 60 • C, 20s; 72 • C, 20s). The specificity of each PCR amplification was determined by melting curve analysis and by analysis in 2% agarose gels. The primer sequences can be found in Supplementary Table 7. The relative transcript expression was calculated by the Pfaffl algorithm, using CqACT2A and CqMON1 as endogenous genes for normalization. CqACT2A and CqMON1 were selected based on their stable expression as reference genes in Arabidopsis (Wallström et al., 2012). Tenfold dilutions of cDNA template were used to determine the amplification efficiency for each gene (Pfaffl, 2001). Primer pairs were designed using Perlprimer (Marshall, 2004) so that one of the primers in each pair spanned an exon-exon border, and the primer pairs were additionally checked using Netprimer (premierbiosoft.com) to avoid primerprimer interactions.

Evolutionary Analysis of GLP
Protein alignments were made using Muscle (Edgar, 2004). Aminoacid substitution models were evaluated by MEGA X (Kumar et al., 2018). The protein evolutionary tree was performed by maximum likelihood using LG model (Le and Gascuel, 2008) with gamma distribution (LG + G) and 95% limit for partial gaps. Total positions in the final dataset were 191. Bootstrap testing was conducted with 1000 replicates.

RESULTS
Transcriptome Sequencing of C. quinoa in Axenic Co-culture With Trichoderma Two cultivars of quinoa were grown axenically and subjected to 36 h treatments with the Trichoderma strains T22 and BOL-12 for gene expression analyses to detect molecular responses. During this short-term treatment, the effect of Trichoderma on seedlings was not measurable, consistent with previous observations (Rollano-Peñaloza et al., 2018).
Five RNA samples from quinoa roots treated with Trichoderma for 12 and 36 h were collected. The three RNA samples at 12 h post inoculation (hpi) showing best quantity and quality parameters were selected for RNA-seq. Similarly, three RNA samples at 36 hpi were selected to be analyzed with qRT-PCR. Sequencing was carried out in paired-end mode, and the final number of reads that passed the quality control varied between 10.2 and 23.1 million paired-end reads of 125 bp per sample (Table 1). Reads were mapped to the draft quinoa genome of cultivar Kd as well as to the chromosome-level assembly of the quinoa genome cultivar QQ74 (Table 1). On average, the proportion of mapped reads was substantially increased when reads were mapped to the QQ74 quinoa genome (93.9 %), as compared to the draft Kd quinoa genome (71.8%). Therefore, all downstream analyses were performed with data mapped to the QQ74 genome.

Differential Gene Expression in Quinoa in Response to Trichoderma Treatment
The differential gene expression analysis considered only reads that mapped to unique locations in the QQ74 genome. The average number of reads that were mapped to unique locations in the QQ74 genome was 89.8% ( Table 1). The remaining reads (10.2%) producing multiple alignments were discarded. Further, only quinoa genes with at least one read in each of the samples analyzed were considered ( Table 2).
Quinoa roots in general induced more genes than they repressed upon interaction with Trichoderma, with the exception of Kurmi interacting with T22 where more genes were repressed than induced ( Table 2). Kurmi treated with T22 showed 16 times more differentially expressed genes than in the treatment with BOL-12. Similarly, quinoa cv. Real treated with T22, compared to the mock-controls, had 5.5 times more differentially expressed genes than Kurmi in the treatment with BOL-12 ( Table 2). Regarding communal effects by both Trichoderma strains, we observed more genes differentially expressed in cv. Real (141 genes) than in cv. Kurmi (75 genes) (Supplementary Tables 1-3). Among the quinoa genes up-or downregulated under one or several conditions, only 19 were communally differentially expressed in all experimental combinations, and all were induced (Figure 1, Table 3). That is, they were significantly induced during the interaction of each quinoa cultivar with each Trichoderma strain. The group of 19 differentially expressed genes were not significantly associated with any functional GO term upon analysis by SEA. However, the structural GO analysis suggests that 13 of the 19 gene products are localized outside the cytoplasm. This suggests an association of the response pattern to functions of the plasma membrane, the cell wall and the extracellular compartment, indicating functions relating to interactions with external stimuli ( Table 3).

Quinoa Genes Differentially Expressed Unique to Each Cultivar
We decided to analyze genes that were induced by Trichoderma and were uniquely expressed in each cultivar. The Kurmi cultivar upon interaction with either BOL-12 or T22, expressed 75 genes that were communally differentially expressed (DE) but were not differentially expressed upon either Trichoderma interaction in cv. Real (Figures 1, 2, Supplementary Table 1). The expression profiles of these genes were clustered by Euclidean distance and are shown in Figure 2. From the 75 DE genes in cv. Kurmi by both strains of Trichoderma, 59 genes were induced (Supplementary Table 1), whereas 16 DE genes were repressed (Supplementary Table 2). The 75 DE genes expressed in cv. Kurmi are expressed in cv. Real but are not responsive to the treatment with either of the Trichoderma strains (Figure 2).
Analysis of the 59 significantly induced genes revealed that 17 genes (CqGLPs) are highly expressed and share a high protein sequence identity (90%). These genes encode proteins that belong to the germin-like protein family (GLPs) (Figure 2;  Supplementary Table 1). Further, several genes involved in flavonoid biosynthesis were specifically responsive in Kurmi. We identified 9 genes whose orthologs in Arabidopsis thaliana are described to be involved in the flavonoid biosynthesis pathway. These differentially expressed genes are orthologs to four out of five enzyme-coding genes necessary for production of flavonol glycosides from naringenin, also known as chalcone (Dao et al., 2011) (Figure 2, Supplementary Table 1). FIGURE 1 | Venn diagram of quinoa genes differentially expressed in response to Trichoderma. Quinoa genes differentially expressed were grouped according to the cultivars and Trichoderma strains studied. The black circle indicates genes differentially expressed in both quinoa cultivars by each of the Trichoderma strains tested. The numbers in parenthesis indicate the number of genes differentially expressed in each of the quinoa-Trichoderma interactions studied.
The Real cultivar had 141 genes differentially expressed common to both Trichoderma strains tested (Figure 1,  Supplementary Table 3). The cv. Real response to both Trichoderma strains showed mostly activation of transcription factors and enzymes without a significant match to a known pathway. Among the genes that were differentially expressed there are four ethylene-responsive transcription factors (ERFs), nine probable WRKY transcription factors and three chitinases (Supplementary Table 3).
Genes differentially expressed related to biotic interactions were observed in a larger proportion in quinoa cv. Kurmi than Real. Therefore, the focus of this study was on the response of the Kurmi cultivar.

Functional Annotation of Differentially Expressed Genes
To assess the function of the differentially expressed genes we annotated the differentially expressed genes of all combinations with Gene Ontology (GO) terms for biological processes. The quinoa genome has 44 776 annotated genes (Jarvis et al., 2017) but the annotation with Argot only assigned GO terms to 50.5 % of the genes (i.e., 22 650 genes annotated with GO terms). Despite the low percentage of GO terms assigned, GO annotation for the All genes commonly differentially expressed in all plant-Trichoderma combinations were induced. a Genes annotated in the Quinoa QQ74 genome (Jarvis et al., 2017) curated with information from their closest ortholog in A. thaliana. b Gene codes from the Arabidopsis Information portal Araport.
*CqWRKY33 was included in the list of genes because significant difference respective to the control was confirmed by qRT-PCR ( Figure 5).
Frontiers in Fungal Biology | www.frontiersin.org FIGURE 2 | Heatmap profile of genes that respond to to Trichoderma BOL-12 and T22 in cv. Kurmi but not in C. quinoa cv. Real. Genes differentially expressed in cv. Kurmi in response to either Trichoderma strain, but not significantly responsive in cv. Real, were analyzed. Clustering by Euclidean distance shows the similarity in expressional change upon Trichoderma treatment. Brown arrows indicate CqGLPs.
Frontiers in Fungal Biology | www.frontiersin.org FIGURE 3 | Correlation of RNA-seq and qRT-PCR gene expression data. Ten differentially expressed genes and two reference genes from the RNA-seq dataset were evaluated by qRT-PCR. Gene expression by qRT-PCR was normalized to the CqAct2 reference gene. Fold change was measured by comparing samples treated with Trichoderma against mock-treated. The selected genes were assessed in all quinoa-Trichoderma combinations as averages of triple biological replicates. The Pearson correlation coefficient between the RNA-seq and qRT-PCR data was 0.921 (For data see Supplementary Table 5).
biological process category in Kurmi plants treated with BOL-12 revealed defense response (GO:0006952) and response to biotic stimulus (GO:0009607) as the main and only processes associated to Trichoderma BOL-12 treatment (Supplementary Table 4). In contrast, the interaction between Kurmi and T22 did not show any significant GO term for biological processes. Quinoa plants of the Real cultivar had more genes associated to GO terms than Kurmi. However, no specific association to a cluster of similar GO terms were observed (Supplementary Table 4). Specifically, Real treated with T22 showed 38 genes that were annotated as stress response genes (GO:0006950) and 6 genes that were annotated to chitin catabolic processes (GO:0006032) and associated redundant GO terms (Supplementary Table 4). Further, Real treated with BOL-12 did not show any GO terms directly associated to defense response or response to stress, yet highest significance was observed to GO terms for cell wall related processes (Supplementary Table 4). Nonetheless, in the interaction between Real and each strain of Trichoderma we observed several genes related to defense being commonly activated. Among them were WRKY transcription factors (nine differentially expressed genes), ethylene-responsive genes (4 differentially expressed genes) and chitinases (three differentially expressed genes) (Supplementary Table 3).

Validation of RNA-Seq With qRT-PCR
Quinoa root transcriptomes have not previously been analyzed. We therefore validated the gene expression data obtained by RNA-seq by performing qRT-PCR for 10 selected genes, including induced, repressed and stably expressed genes (Figure 3, Supplementary Table 5). A log-log linear model analysis of the RNA-seq data and the qRT-PCR data showed a strong correlation (R 2 ) of 0.848. The correlation was higher when the different quinoa-Trichoderma interaction pairs were assessed independently (Supplementary Table 5). Real root samples treated with Trichoderma BOL-12 and T22 strains were assessed at 36 hpi by qRT-PCR. mRNA levels were normalized to the CqAct2 reference gene. Fold changes (mean ± SE) were determined by comparing samples treated with Trichoderma against mock-treated ones. Asterisks denote significant changes in the gene expression as compared to the control treatment by qRT-PCR at 36 hpi (p < 0.05). Symbol π denotes significant changes in the gene expression as compared to the control treatment, using RNA-seq (p < 0.05) and confirmed by qRT-PCR (p < 0.05) at 12 hpi.

Changes in Root Gene Expression at 36 hpi
Time changes in the expression of quinoa genes by Trichoderma treatment were assessed by qRT-PCR at 36 hpi. We followed time-dependent changes in two highly induced genes (CqGLP1 and CqGLP10) representing the GLP family and one gene that was induced in all quinoa-Trichoderma interactions (CqHSP83). The gene expression of CqGLP1 and CqGLP10 was reduced at 36 hpi compared to 12 hpi in the Kurmi cultivar but its expression was still higher than the mock-treatment. In contrast, the gene expression of CqGLP1 and CqGLP10 in the Real cultivar was higher at 36 hpi than at 12 hpi, being statistically significant in the Real -BOL-12 interaction (Figure 4). The CqHSP83 gene maintained its level of gene expression between 12 and 36 hpi by application of T22 in both cultivars. In contrast, the application of BOL-12 to both cultivars downregulated CqHSP83 gene expression at 36 hpi as compared to 12 hpi. However, the downregulation was only significant in the Kurmi cultivar (Figure 4). Overall, the results suggest that the induction of the analyzed genes is slower in cv. Real than in cv. Kurmi. FIGURE 5 | Gene expression changes in quinoa shoots after treatment of roots with Trichoderma. Quinoa shoot samples were assessed by qRT-PCR 36 h after treatment with Trichoderma to the roots. Gene expression was normalized to the CqAct2 reference gene and is shown as mean ± SE. Asterisks denote significant changes in the gene expression as compared to the control treatment (p < 0.05). Two CqGLP genes that were induced in Kurmi roots at 12 hpi as well as the CqHSP83 gene, which was induced in roots in all quinoa-Trichoderma interactions at 12 hpi are shown.

Shoot Gene Expression
We investigated changes in the quinoa shoot gene expression 36 h after Trichoderma application to the root neck (Figure 5). Ten out of 12 genes investigated were also expressed in the shoots, CqPER39 and CqPR1C gene expression was not detected in the shoots in any of the combinations studied (Supplementary Table 6). Trichoderma-induced gene expression changes in the shoots (Figure 5) showed a generally similar pattern of gene expression as observed in the roots at 36 hpi (Figure 4). CqGLP1 and CqGLP10 are significantly expressed in both cultivars upon interaction with BOL-12 but not with T22. Likewise, CqHSP83 is significantly expressed in both cultivars when interacting with T22 but not when interacting with BOL-12 ( Figure 5). The other genes did not show a significant correlation in the shoot-root expression in any of the quinoa-Trichoderma interactions studied (Supplementary Table 6).

Evolutionary Analysis of the Germin-Like Proteins
Plant germins were first investigated and have been characterized in most functional detail in cereals (Ilyas et al., 2016). To investigate the coincidental induction of germin-like proteins in cv. Kurmi (Supplementary Table 1), we carried out BLAST searches to identify all germin and germin-like homologs in C. quinoa, Beta vulgaris, A. thaliana and Hordeum vulgare and performed alignments and evolutionary analyses. We found that 16 of the 17 quinoa GLPs induced by Trichoderma (highlighted in green) in cv. Kurmi belong to a single (98% bootstrap) C. quinoa-specific clade of 29 homologs (Figure 6, Supplementary Table 8). The remaining GLP induced by Trichoderma (CqGLP20) is grouped in an unresolved putative clade, which contains four non-expressed quinoa GLPs and one sugarbeet GLP (Figure 6). Quinoa germin-like proteins were significantly associated with specific homologs in B. vulgaris. Species-specific gene groups were also observed for B. vulgaris and A. thaliana. The result suggests that recent expansions of gene groups have occurred independently in the amaranth family species C. quinoa and B. vulgaris.
Predicted structure of quinoa GLP1 shows a highly conserved structure compared to the crystal structure of Barley Germin (PDB: 1FI2), which has oxalate oxidase and manganese superoxide dismutase activities. The model suggest that the protein active site with aminoacid residues that bind to manganese is conserved and thus CqGLP1 might have a potential enzymatic activity (Figure 7).

DISCUSSION
The outcome of plant-Trichoderma interactions with respect to both growth and physiological changes has been shown to be genotype-specific regarding both the plant and the Trichoderma biomaterial (Harman et al., 2004;Tucci et al., 2011). Here, we have observed a small set of quinoa genes being responsive in all combinations of Trichoderma strains and quinoa cultivars studied. However, we have found many more genes that are differentially expressed by a specific quinoa cultivar in response to either or both of two Trichoderma strains (Figure 1).
The set of 19 genes that showed significant responses in all cultivar-strain combinations mainly include genes connected to biotic stress response and cell wall modification ( Table 3). Orthologs of these genes are known to be involved in defense response. For example, the polygalacturonase inhibitor protein AtPGIP1 (AT5G06860) in A. thaliana is thought to inhibit cell wall pectin degrading enzymes, commonly produced by fungal pathogens (Lorenzo and Ferrari, 2002). Xyloglucan endotransglucosylase/hydrolases are cell wall repairing enzymes, many of which are induced by fungal infection (Miedes and Lorences, 2007). Further, two highly similar (protein sequence identity of 94.3%) heat-shock proteins (CqHSP83A and CqHSP83B) annotated to be involved in general stress responses (GO:0006950) were upregulated in all quinoa-Trichoderma interactions. Large heat-shock proteins (70-90 kDa) are known to be involved in plant defense response through the stabilization of protective plant proteins (Hubert et al., 2003;Takahashi et al., 2003;Zhang et al., 2015). The heat-shock proteins expressed by quinoa might have been induced to contribute to the stabilization of defense proteins that would prevent or counteract damages induced by Trichoderma (Rollano-Peñaloza et al., 2018).
Several differences in the defense response patterns between the quinoa cultivars were observed. Especially, the Kurmi cultivar displays specific activation of several homologs to biotic stressassociated plant genes (Supplementary Tables 1, 2). In contrast, the responsive genes in cv. Real were mostly involved in general cellular processes (Supplementary Table 3), and to a lesser extent involving defense response genes. The defense response gene set induced in cv. Real was also completely different from the one activated in Kurmi (Figure 1, Table 2, and Supplementary Tables 1-3). Surprisingly, in neither case an obvious association to known major pathogen response pathways like jasmonic acid, salicylic acid or ethylene pathways (Pieterse et al., 2009;Salas-Marina et al., 2011;Mathys et al., 2012) could be observed in the GO analysis at 12 hpi. The low association of the quinoa differentially expressed genes with these known pathways could, however, be caused by the relatively low level of GO annotation observed for the quinoa genome, which limits the statistical discovery power of the GO analysis for quinoa genes.
The quinoa genes specifically induced in the Kurmi cultivar upon interaction with Trichoderma (Supplementary Table 1) resemble a set of defense response genes observed in plantpathogenic interactions (Dixon, 2001). Several of these induced genes (chalcone synthase, chalcone isomerase, flavonol synthase, UDP glycosyl transferase and cis-zeatin O-glucosyltransferase) belong to the flavonoid biosynthetic pathway (Dao et al., 2011). Flavonoids have an important role in plant defense (Dixon, 2001;Treutter, 2006). For example, some A. thaliana mutants lacking the UDP glycosyl transferase gene (AtUGT74F1) are more susceptible to Pseudomonas syringae infection than the wild-type (Boachon et al., 2014). Further, a QTL analysis for pathogen resistance in soybean identified two UDP glycosyl transferase genes as the candidate genes responsible for resistance to Fusarium (Cheng et al., 2017). Thus, the Kurmi cultivar might be producing flavonol glycosides in order to prevent damage from Trichoderma overgrowth.
The Kurmi cultivar specifically induced several plant defensins that belongs to the germin and GLP family (Figure 6 and Supplementary Table 1). The Trichoderma-responsive GLPs form a majority (16 genes out of 29) of a recently expanded quinoa-specific clade (Figure 6, Supplementary Table 8), which is thus strongly connected to Trichoderma interaction. Two of the GLP genes were further tested, and were found to be also induced in leaves (Figure 5). The timing of the induction further indicated that GLPs are induced in both Kurmi and Real, albeit more slowly in Real (Figure 4). Given that several quinoa GLPs were expressed upon interaction with Trichoderma, it is very likely that these GLP defensins have an important role in the plant immune response of quinoa. Resistance to microbe attacks has been previously connected to the speed of response in different cultivars (Pritsch et al., 2000;Venisse et al., 2002). Thus, the rapid induction of a cluster of GLPs in both roots and leaves of Kurmi compared to Real (Supplementary Tables 2, 3;  Figure 4) makes these genes potential candidates for breeding to increase the tolerance to microbial attacks in quinoa plants.
The set of defense-related genes (peroxidases, chitinases, ERF and WRKY transcription factors) induced in quinoa cv. Real upon interaction with Trichoderma (Supplementary Table 3) has been observed in L. japonicus roots upon 1 h of incubation with chitin oligosaccharides (Nakagawa et al., 2011). However, FIGURE 7 | I-TASSER model of the C. quinoa GLP1. Whole quinoa germin-like protein (red) aligned to the template crystal structure of barley germin (oxalate oxidase; pdb:1FI2; cyan). The manganese ion of the active site is shown as a purple sphere. Inset shows the alignment integrated in the hexameric germin structure with neighboring germin monomers colored yellow.
the levels of such defense-related genes returned to normal after 7 h in L. japonicus whereas in quinoa remained induced after 12 h, possibly due to the persistance of interaction with the living Trichoderma agent as compared to the transient nature of the elicitor. Similar to the L. japonicus system, the 24 h interaction of Trichoderma with A. thaliana resulted in the induction of the same defense-related genes as seen here in quinoa (Brotman et al., 2013). This set of genes could thus be a basal gene response of plants after recognition of beneficial fungi like Trichoderma through chitin. In contrast, the Kurmi cultivar might have a different set of receptors that helps the plant to perceive a possible negative effect from Trichoderma and thus rapidly activate a different set of defense-related genes.
The plant root transcriptomic response to Trichoderma has been poorly studied compared to aerial parts (Vos et al., 2015). Nevertheless, it has been observed that in tomato roots the recognition of Trichoderma at 24 hpi activates ROS signaling, SA responses, cell wall modifications (Palma et al., 2019), JA responses and induction of plant defenses (Brotman et al., 2013;Ho et al., 2016). In our study, we have observed a similar pattern for ROS signaling, cell wall modifications and induction of plant defenses (Table 3, Supplementary Tables 1-3), confirming that the first response of root plants to beneficial fungi like Trichoderma is to activate defenses. Further, our study reveals that the defense response against beneficial fungi is variable between cultivars (Figure 2). The variable molecular response between cultivars to Trichoderma could help to create molecular markers of compatibility between certain plant cultivars and certain strains of Trichoderma.
In conclusion, our study suggests that Trichoderma triggers a defense response in quinoa plants. Comparing the defense response of two quinoa cultivars we can observe that the Kurmi cultivar mainly induced a set of genes involved in plant defense. In contrast, the Real cultivar did not have a clear response because most of the changes mediated by Trichoderma were related to general stress and regulation of biological processes. The Kurmi cultivar might have higher tolerance to microbe attacks due to the expression of genes involved in the biosynthesis of flavonol glycosides and a clade of GLP-defensins unique to quinoa. These genes are thus candidates for selection of quinoa cultivars with higher resistance to microbe attacks.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/bioproject/PRJNA720675.

AUTHOR CONTRIBUTIONS
AR and OR-P designed the study, analyzed the data, and wrote the paper with input from all authors who reviewed and approved the final manuscript. OR-P performed the experiments. All authors contributed to the article and approved the submitted version.

FUNDING
This work was funded by the Swedish International Development Agency (SIDA) in a strategic collaboration between Universidad Mayor de San Andrés (UMSA) (Bolivia) and Lund University (Sweden), under the SIDA-UMSA Cooperation Agreement 75000554; the Lars Hiertas Minne Stiftelsen (OR) and the Jörgen Lindströms Stipendiefond (OR). Genome mapping, read count and blastn computations were performed on resources provided by SNIC (Swedish National Infrastructure for Computing) through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX, Uppsala, Sweden) under Project SNIC b2015011/b2016265 along with resources provided by UMSA (Universidad Mayor de Sán Andres) through Centro Nacional de Computación Avanzada para Bioinformática y Genómica (CnCaBo, La Paz, Bolivia).