Transcriptome changes induced by arbuscular mycorrhizal fungi in sunflower (Helianthus annuus L.) roots

Arbuscular mycorrhizal (AM) fungi are essential elements of soil fertility, plant nutrition and productivity, facilitating soil mineral nutrient uptake. Helianthus annuus is a non-model, widely cultivated species. Here we used an RNA-seq approach for evaluating gene expression variation at early and late stages of mycorrhizal establishment in sunflower roots colonized by the arbuscular fungus Rhizoglomus irregulare. mRNA was isolated from roots of plantlets at 4 and 16 days after inoculation with the fungus. cDNA libraries were built and sequenced with Illumina technology. Differential expression analysis was performed between control and inoculated plants. Overall 726 differentially expressed genes (DEGs) between inoculated and control plants were retrieved. The number of up-regulated DEGs greatly exceeded the number of down-regulated DEGs and this difference increased in later stages of colonization. Several DEGs were specifically involved in known mycorrhizal processes, such as membrane transport, cell wall shaping, and other. We also found previously unidentified mycorrhizal-induced transcripts. The most important DEGs were carefully described in order to hypothesize their roles in AM symbiosis. Our data add a valuable contribution for deciphering biological processes related to beneficial fungi and plant symbiosis, adding an Asteraceae, non-model species for future comparative functional genomics studies.

Differentially expressed genes at the early mycorrhizal stage. Comparing H. annuus control with mycorrhizal roots 4 days after inoculation (D4), we detected only 19 significantly over-expressed genes. No significant under-expressed transcripts were found in D4 versus control plants. Within biological processes, the most abundant GO-terms of D4 differentially expressed genes were "metabolic process" and "cellular process"; among molecular functions, the most frequent terms were "binding", "hydrolase activity" and "transporter activity"; cellular component class showed "membrane" and "cell wall" as major terms (Fig. 3).
Amongst differentially expressed genes of D4 plants we identified several transcripts encoding stress response proteins such as germin-like proteins and chitinase, and genes encoding polypeptides with protease activity (Subtilisin-like and carboxypeptidases). Other DEGs included those encoding a nrt1 protein of the ptr family, two UDP-glycosyltransferases, a lysM domain receptor-like kinase 3, an inorganic phosphate transporter, and a glutathione s-transferase (Supplementary Table 1).  Differentially expressed genes at the late mycorrhizal stage. At the second harvesting time point, 16 days after sunflower root inoculation with the AMF R. irregulare (D16), the comparison between control and mycorrhizal roots showed 694 over-and 13 under-expressed genes, respectively ( Supplementary Fig. 1). The greatly increased number of DEGs in D16, compared with D4, may be ascribed to the increasing level of AM root colonization. Furthermore, the number of over-expressed DEGs in D16 plants was considerably higher than that of under-expressed DEGs ( Supplementary Fig. 1).
GO terms analysis on D16 DEGs showed that "catalytic activity" and "hydrolase activity" ontologies were the most frequent in the Molecular Function class; the most represented Biological Process terms were "metabolic process" and "cellular processes"; "membrane" and "cell wall" were the two top terms in Cellular Component class (Fig. 3).
Some under-expressed genes were found in mycorrhizal D16 roots, although fewer than the over-expressed transcripts of the same root sample. Regarding Biological Processes, the most represented term was "metabolic process" (2 genes); within Cellular Components the top term was "membrane"; "transporter activity" was the most frequent term in the Molecular Function class (Fig. 3).
Gene ontology enrichment analysis showed that only the group of over-expressed genes in D16 was significantly enriched (P < 0.05). Using REVIGO, "hydrolases", "oxidoreductase" and "catalytic activity" categories were enriched in the late stage of colonization: "catalytic activity" were the most enriched among Molecular Functions, and "metabolic process" among Biological Processes (Fig. 4). The most enriched term in the Cellular Component class was "membrane" (Fig. 4).
In addition, among other over-expressed transcripts at D16 with logFC < 11 were those encoding for proteins involved in stress-defense responses (heat shock proteins, several hydroxyproline-rich glycoprotein, peroxidases, and several lectins) and putative proteins of root remodeling and growth (an acidic endochitinase, two auxin response factors and alpha-tubulins, CDPK, and several ABC transporters) (Supplementary Table 1).
Amongst under-expressed genes of H. annuus mycorrhizal roots at D16 we detected many transcripts encoding for SWEET-like bidirectional sugar transporters, one transcript belonging to a pathogen-related protein and a peroxidase gene. (Supplementary Table 1).
Noteworthy nineteen DEGs were found to be shared among mycorrhizal roots of H. annuus after 4 and 16 days of culture with R. irregulare. All shared transcripts showed an increasing log fold change between the first and the second time point (Table 2).
qRT-PCR validation analysis. In order to validate RNA-seq data, the expression of five genes selected among those differentially expressed between M4 and M16 samples were analysed by qRT-PCR experiments. As shown in Supplementary Figure 2, all genes were differentially expressed, confirming the accuracy of transcriptome profiling obtained with RNA-seq.

Discussion
A large number of terrestrial plant species establish symbioses with AMF, leading to morphological and functional changes in roots. Next generation sequencing techniques can be conveniently used for analysis of transcriptome profiling during plant/fungus interactions. In this work, RNA-Seq experiments in roots of sunflower after R. irregulare inoculation gave a comprehensive view of changes in gene expression, confirming the activation of peculiar AMF regulated genes in sunflower and identifying new genes so far not associated with the establishment of AMF symbioses. Here, we detected a total of 19 and 694 genes over-expressed in mycorrhizal roots at early and late stage of colonization, respectively. Number of DEG is lower to the one retrieved during AM symbiosis in Lotus japonicus and Oryza sativa 10,14 . Differences may be related to slightly lower percentage of hyphae colonization in sunflower roots compared to these studies, to different experimental conditions used and to species specific responses. Among new genes which were not described so far as involved in AMF symbiosis, we identified, for example, those encoding several BAHD-acyltransferases, an oxoglutarate dependent dioxygenase, and a stellacyanin-like protein (Supplementary Table 1).
Overall GO analysis of DEGs showed a number of ontologies: "biosynthetic process", "hydrolase activity" and "membrane" were the most frequent functional classes. In particular, the GO term "membrane" reflects the major cell membrane rearrangement associated with fungal invasion, "biosynthethic process" suggests host transcriptomic, metabolic and transport activity alteration during symbiosis, and "hydrolase activity" the wide protein recycling associated to cell wall modification (Figs 3 and 4).
Global gene expression data show that D16 mycorrhizal roots were far more active in transcription than D4 root samples. Furthermore, we noticed an unbalanced number of over-expressed genes in comparison with  under-expressed genes, as described for other mycorrhizal plant species, for example Lotus japonicus 14,49 . Indeed, very few transcripts were significantly under-expressed in sunflower AM roots.
One of the major obstacles in studying AMF symbiosis is the difficulty in synchronizing the initial developmental events leading to the establishment of the symbiosis 26 . At the early stage time-point (D4), low AMF colonization, mainly represented by entry points with a few developed arbuscules, occurred in inoculated sunflower roots. Consequently, it can be suggested that DEGs in D4 roots could represent genes regulated in the early phases of fungal colonization. During the late phase of symbiosis establishment (D16), when half the length of host roots were colonized by R. irregulare, a higher number of DEGs were found compared to D4, reflecting the advanced stage of root colonization. DEGs were differentially expressed at both time-points, although they were modulated between D4 and D16. This finding is in accordance with the asynchronous nature of the mycorrhizal colonization process, as D16 samples showed fungal structures representative of both early and late stages of the symbiosis, as new entry points may be produced by extraradical hyphae, giving rise to secondary colonization 26 .
Differentially expressed genes of H. annuus roots in symbiosis with R. irregulare, shared between different time points, encompass genes of symbiosis development, maintenance and functioning and belong to a large number of GO terms. Many of these GO terms and specific genes will be discussed hereafter in relation to their modulation, expression and function.
Differentially over-expressed genes in D4 roots mainly belong to the functional classes "response to stress", "cell wall" and "transporter activity", which are largely represented during early interactions between plants and AMF 11,28 . Interestingly, all D4 DEGs were differentially expressed also in D16 plants, with increasing expression levels at the latter time point (Fig. 3). Indeed, among DEGs involved in plant defense, 3 and 15 transcripts of Germin-like protein were detected in D4 and D16 plants ( Table 2) respectively, suggesting a modulation of this multi-gene family as AMF colonization proceeds. Similarly, the number of induced protease-encoding genes increased from 4 to 27 (including genes encoding subtilisins and carboxypeptidases) between the two time points, confirming the involvement of this family in maintaining cell functions and recycling cell wall proteins during mycorrhizal symbioses 14,50 .
A similar induction profile was observed also for genes involved in transporter activity, which increased from 4 to 60 between early and late stage of colonization (Fig. 3, Table 2), confirming the importance of transporter genes modulation during the symbiotic interaction 11 . Among these DEGs, those encoding an NRT1 PTR protein, involved in NO 3 -and auxin transport, an inorganic phosphate transporter and the phosphate and water transporter nodulin-26-like protein were confirmed to be expressed at both early and late stages of symbiosis establishment 11,51,52 .
A LysM-domain receptor kinase encoding gene is the only signal transduction gene identified amongst D4 DEGs. The number of DEGs of this family increased to 7 in D16 roots, including some leucin rich receptors ( Table 2, Supplementary Table 1). LysM is required for the earliest cellular and physiological symbiotic responses and its expression was detected in M. truncatula roots 8 days post inoculation 28,53 .
UDP-glycosyltransferase is one of the stress-responsive gene belonging to the GO term "transferase activity" and was over-expressed even in D4 AM roots, increasing its logFC between D4 and D16 ( Table 2). Its activity should combat oxidative stress and pathogen infection. Genes of this family were reported to be differentially expressed in AM plants 14,20,54 , but they were never previously described in association with the early steps of mycorrhizal establishment.  DEGs observed only at D16, we will discuss them in relation to their functional classes. The GO term "response to stress" contains several genes encoding transcripts involved in plant stress and defense, and is one of the most represented among DEGs (Figs 3, 4). It has been established that common defense pathways are activated in both pathogenic and beneficial plant-fungal interactions, as a result of fungal elicitors recognised by the plant 28,53 .
New transcripts implied in fungus-plant relationship are BAHD-acyltransferases, which are here described for the first time in AM symbioses, and a transferase super-family including Vinorine synthases, which play a central role in the endogenous formation of specialized monoterpenoid secondary metabolites, functioning as defensive compounds upon pathogenic infections 55,56 . Accumulation of terpenoids belonging to phytoalexins was also detected in mycorrhizal roots of Solanum tuberosum and Solanum lycopersicum 33,57 .
Another plant defense related gene, encoding Glutathione-S-transferase, is active in M16 sunflower roots. Its activation was already observed in Medicago truncatula 58 .
Other over-expressed transcripts encode phytyltransferases; these genes play a role in the isoprenoid metabolic pathway. Usually, isoprenoids are highly expressed as secondary metabolites during AM formation, conferring resistance to stress with activation of plant defense system [59][60][61] .
"Catalytic activity" is another wide functional class involved in plant defense, cell wall loosening and oxidative stress defense. DEGs of this functional class encode for germin-like proteins, the monocopper oxidase SKU5 and several peroxidases 62,63 (Supplementary Table 1), all known to be induced during AM symbiosis 13,29,49,64,65 .
Other DEGs encoding defense and stress-responsive proteins were found in mycorrhizal sunflower roots after 16 days of culture. These genes encode pathogen-related proteins, such as chitinases and endoglucanases, and several peroxidases, which were shown to be differentially expressed in other mycorrhizal plant species [66][67][68][69] .
Other defense-related putatively encoded lectins and enzymes involved in isoprenoid and terpenoid metabolism were reported to accumulate specifically in mycorrhizal roots 23 . In connection with isoprenoid and terpenoid metabolism, an over-expressed gene in D16 mycorrhizal roots encoded a beta-amyrin 28-oxidase that catalyzes the oxidation of beta-amyrin to oleanolic acid 70 .
Another new transcript involved in symbiosis events is an oxoglutarate dependent dioxygenase encoding gene. This gene, for the first time associated with AM symbiosis, was over-expressed in D16 mycorrhizal roots. Oxoglutarate dependent dioxygenases are involved in alkaloid production, whose content (hyoscyamine and scopolamine) was previously reported to increase in mycorrhizal roots of Datura stramonium 71 .
It is interesting to note that some pathogen-related and disease-resistance proteins encoding genes were significantly down-regulated in sunflower AM roots 16 days after inoculation (Supplementary Table 1). Although a partial overlap between symbiotic and pathogenic signaling has been reported 72 , transcript levels of some of these genes were much lower in plant roots colonized by beneficial fungi than in those infected by root pathogens. This reduced expression was associated with a very restricted localization of transcripts in arbuscule containing cells 73 . The repression of transcription of some plant defence genes suggests that AMF are able also to trigger mechanisms in plant roots that suppress defence response allowing colonisation to proceed.
Reception of molecular signals and signal transduction are crucial biological events for symbiosis establishment. DEGs with high logFC related to signal formation and transduction encode P450 oxidases (Supplementary  Table 1), which are known to be involved in the production of strigolactones and co-expressed both in the plant and in the fungus during AM symbiosis 14,[74][75][76] .
Signals released by the fungus activate transcription factors in root cells. Among highly over-expressed putative transcription factors, we detected genes encoding a PHR1, several DELLA proteins involved in gibberellin signaling regulation in AM [77][78][79][80] (Supplementary Table 1), and many ethylene responsive factors (ERFs). PHR1 is a transcription factor belonging to the MYB family that binds to promoters of most genes positively or negatively affected by Pi starvation 81,82 , hence it may be considered a candidate gene product involved in Pi transport regulation. Ethylene-responsive factors and the ERF Required for Nodulation (ERN) transcription factor are induced in response to AMF, and ERN-like proteins encoding genes showed specific expression in arbusculated cells 23 .
Some over-expressed transcripts are related to fatty acid metabolism, as for example the GDSL esterase/lipase. The role of these proteins is not fully known, but their remarkable expression during AM symbiosis, mainly in arbusculated cells, has been correlated with the production of signal molecules and membrane components involved in the synthesis of periarbuscular membranes 14 .
A transcript encoding a calcium-dependent protein kinase 1-like protein was over-expressed in sunflower D16 mycorrhizal roots, compared with controls (Supplementary Table 1). Calcium-dependent protein kinases (CPKs or CDPKs) are important Ca 2+ sensors in signaling processes during the establishment of symbioses between host plants and both AMF and rhizobia 83,84 .
As reported above, previous studies suggested that hormonal balance of mycorrhizal roots is tightly linked to the development of fungal colonization and root structural changes. The expression pattern of genes putatively encoding for auxin and gibberellin biosynthesis showed that transcripts encoding GA oxidases and auxin response factor were over-expressed in D16 mycorrhizal roots (Supplementary Table 1).
Gibberellin content increased in the roots of mycorrhizal plants compared with uninoculated controls in parallel with the accumulation of transcripts involved in their biosynthesis [85][86][87] .
The role of auxins in mycorrhizal symbioses is not clear, as biochemical and molecular studies showed variable auxins and auxin-related molecules levels, although data on auxin mutants suggested that auxin perception is required for arbuscule development and mycorrhizal symbiosis establishment, at least by modulating strigolactone levels 79,88,89 .
A transcript homologous to zeatin O-glucosyltransferase differentially accumulated in D16 mycorrhizal roots (Supplementary Table 1 Many over-expressed transcripts putatively encoding transporter proteins were found in both D4 and D16 mycorrhizal roots, including an inorganic P transporter, a multidrug resistant ABC transporter, nitrate and ammonium transporters, and aquaporins (Supplementary Table 1). These genes are known to be highly expressed during mycorrhizal interactions and play a pivotal role in symbiotic signaling and nutrition pathways [90][91][92][93] .
ABC transporter genes are required for arbuscule development 94 and were over-expressed in D16 mycorrhizal roots samples.
In this work, genes encoding inorganic P, ammonium and nitrate transporters and NIP class of aquaporins were over-expressed in D4 and D16 mycorrhizal roots, and a larger number of highly expressed isoforms was detected at the late stage of the symbiosis, confirming previous data 16,29,95,96 . The nitrate transporter Nrt1 was shown to play a pivotal role also in auxin transport and consequently in lateral root formation 52 .
Another newly identified transcript is related to an aquaporin-like encoding gene over-expressed in D16 AM roots that encodes for a stellacyanin-like protein, exhibiting sequence similarity with nodulin-like transporters 97 (Supplementary Table 1). A gene encoding for a putative stellacyanin-like protein, CASLP1, was strongly expressed in Capsicum annuum plants after inoculation with different pathogenic bacteria and fungi, possibly due to water losses or to a role of this protein in limiting pathogen spread 98 . Although at a lower level of over-expression, a probable tonoplast intrinsic protein (TIP) aquaporin was also detected in D16 mycorrhizal roots.
A 9-fold over-expressed sequence recovered from D16 mycorrhizal roots showed homology to a bidirectional SWEET16-like sugar transporter (Supplementary Table 1), characterized as an exporter of sucrose and monosaccharides 99 , and could play a role in regulating the transfer of plant carbon resources at the plant-fungus interface. Potato SWEETs showed differential expression regulation in response to AMF, although further studies are needed to clarify SWEETs' role in sugar partitioning during mycorrhizal colonization process 100 .
Plant cell walls form a dynamic extracellular matrix that controls growth and development and, at the same time, mediates most plant/fungal interactions. During AM symbiosis the cortical layers are reached by fungal hyphae branches that, after contact, lead to unique structures called arbuscules, and this process requires the activation of several genes 101 . Many transcripts involved in cell wall and membrane shaping, functional to hyphal penetration and arbuscule development, were highly expressed, such as that encoding SKU5 protein, which exhibits sequence similarity with multiple copper ascorbate oxidase and laccase 102 (Supplementary Table 1). Other over-regulated transcripts encode carboxypeptidases, subtilisin-like proteases and a type 1 ATPase, all belonging to the GO term Hydrolase and involved in cell remodeling and membrane transport. Here a considerably higher logFC of the expression of some subtilisins was detected in D16 mycorrhizal roots, compared with D4. This result is similar to that reported in roots of L. japonicus inoculated with R. intraradices 50 .
In D16 mycorrhizal roots, transcripts encoding proteins functioning in cell wall and membrane shaping were detected, such as a ripening-related protein, a Clathrin vesicles protein, some vacuolar processing enzymes and Ankyrin-repeat containing proteins, as previously reported 52,103,104 (Supplementary Table 1).
EG45 domain is typical of the expansin family 105 , a superfamily of proteins that play a pivotal role in cell wall loosening, thus increasing root cell wall plasticity. Expression of these proteins is triggered by fungal contact or by a diffusible fungal factor and is a prerequisite to accommodate the fungus in plant root cells 27,106 . Expansin-encoding genes showed up-regulation in mycorrhizal roots of M. truncatula and C. sativus, and their products were detected both in the interface matrix and in mycorrhizal root cell walls 28,107 . A gene encoding a protein similar to Prunus armeniaca major allergen, over-expressed in D16 mycorrhizal roots, showed a high similarity to a β-expansin 108,109 (Supplementary Table 1).
Many D16 over-expressed genes belong to the GO term "carbohydrate metabolic process" (Fig. 2). Genes included in this GO, such as those encoding endoglucanases, galacturonases, hydroxyproline-rich glycoproteins, cellulose synthase and ascorbate oxidase, are mainly related to root cell wall degradation and remodeling and were expressed in root cells containing fungal structures 101,110,111 .
Tubulin encoding genes were over-expressed in D16 mycorrhizal roots (Supplementary Table 1), confirming previous findings on tubulin upregulation in S. lycopersicum roots colonized by R. intraradices, compared with uncolonized roots 112 .
Plant-fungal symbiosis stimulates new adventitious roots contributing to lateral root growth as showed in symbioses of Allium porrum and Prunus cerasifera with F. mosseae and R. intraradices 113,114 . Auxin response factors (ARF) are transcription factors involved in adventitious root formation and regulate auxin responsive genes. The auxin response factor 18-like transcript detected here (Supplementary Table 1) was not reported before in mycorrhizal plants, whereas other transcripts encoding putative auxin response factors were induced during colonization of O. sativa and M. truncatula by R. irregulare 18,65 . Such an expression pattern might be related, as in other ARFs, to lateral root formation and growth 115 .

Conclusions
In this work, a comprehensive transcriptome analysis of roots of sunflower in symbiosis with the AM fungal symbiont R. irregulare was carried out adding data from a plant species belonging to Asteraceae, the largest Angiosperm family, which was not investigated yet in this respect. Hence, this is a valuable contribution to deciphering gene expression related to plant symbiosis, for future comparative functional genomics analyses. Results confirmed the importance of AM-modulated transcription of genes involved in plant defence, signal transduction, hormonal balance, nutrient transport and host cell shape in the activation of mechanisms controlling AMF establishment and spread in plant roots.
The repertoire of genes regulated at both early and late stages of the symbiosis allowed us to confirm AMF-regulated genes as well as to identify new genes which were not previously described as involved in AMF symbiosis, for example those encoding a BAHD-acyltransferase, an oxoglutarate dependent dioxygenase, and a stellacyanin-like protein. Such results can guide application research targeted at exploiting metabolic changes occurring in commercially-valuable crops when colonized by arbuscular mycorrhizal symbionts in field conditions, in pursuit of sustainable and resilient agricultural processes.
Further studies are in progress to elucidate gene expression in the development of AM symbiosis by the comparison of transcriptome profiling in cultivated and wild genotypes of sunflower, showing different susceptibility to AMF.

Materials and Methods
Plant material and treatments. Fourty seeds of H. annuus HA412-HO inbred line (USDA acc. number PI 642777) were germinated on moistened filter paper in Petri plates. The root system of one week old plantlets was individually placed between a 90-mm diameter membrane (cellulose acetate and cellulose nitrate mixture, 0.45 µm pore diameter size, MF-Millipore ™ ) and a 100-cm 2 nylon net (41 µm mesh, Millipore ™ ). Spores, mycelium and finely cut colonised roots of R. irregulare (formerly Rhizophagus irregularis) IMA6, obtained from pot-culture soil after wet sieving through a 100-µm-mesh size sieve, were spread onto all the nylon net surface of 20 plantlets. Uninoculated control plantlets received 10 ml of a filtrate, obtained by sieving the same pot-culture soil through a 50 μm pore diameter sieve and through Whatman N. 1 paper, to ensure a common AM fungal-associated microflora. In both inoculated plants and uninoculated controls another 90-mm membrane was placed on the nylon net, and these "sandwich systems" were placed in sterile 150-mm Petri plates containing steam-sterilized quartz grit. After sealing with Parafilm M, the lower half of each plate was wrapped in aluminium foil and plants were maintained in a growth chamber at 24 °C. Plants were supplied weekly with 6 ml half strength Hoagland's solution. The progression of mycorrhizal colonization was monitored across time-points 48 hours apart on 3 replicate plants, by localizing succinate dehydrogenase activity 116 and by staining roots with Trypan blue 117 . Percentage of AM fungal root colonization was determined by the gridline intersect method 118 .

RNA extraction.
On the basis of root colonization data, RNA-seq analyses were carried out on roots collected from 3 replicate H. annuus mycorrhizal and control plantlets harvested after 4 and 16 days of culture.
Whole roots were ground in liquid nitrogen and total RNA was isolated using the Logemann procedure 119,120 . Purification from genomic DNA was performed by digestion with DNaseI (Roche). Finally, RNA was purified with phenol/chloroform and precipitated with standard procedures. RNA quality was evaluated by electrophoretic, spectrophotometric and bioanalyzer analysis using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA). High quality RNA-Seq reads were mapped using the CLC genomics workbench (v. 9.5.3). Traces of ribosomal RNA contamination were removed from all libraries by mapping against sunflower ribosomal sequences downloaded from the NCBI database repository. The parameters used for ribosomal filtering were: length fraction = 0.5 and similarity fraction = 0.8. cDNA sequences were mapped to the reference transcriptome of Helianthus annus inbred line XRQ 38 (https:// www.heliagene.org/HanXRQ-SUNRISE/) using stringent parameters (length fraction = 0.8 and similarity fraction = 0.8) but retaining mismatch penalties = 2 and gap penalties = 3 as set by CLC default parameters. Sequence reads are deposited in the Sequence Repository Archive (SRA) under the ID SRP124596. Differential expression analysis. Gene expression level was calculated as reads per kilobase per million reads mapped (RPKM) as described in Mortazavi et al., (2008). Besides unique reads, reads that occurred up to ten times were also included in the mapping analysis, because this strategy correctly estimates the expression of paralogous genes 122 . Gene expression was filtered to ensure RPKM > 1 in at least one library.

RNA sequencing and mapping procedures.
Raw counts deriving from CLC aligner were analyzed with the R statistical package EdgeR 123 . A pairwise comparison test was performed between mycorrhizal (M) and control libraries (C) after 4 and 16 days of inoculation with fungus (C4 vs. M4 and C16 vs. M16). The resulting P values were corrected for the FDR 124 and significant genes identified with an FDR corrected P value < 0.05. The fold changes (FC) between controls and treatments were considered significant when the expression value of a gene was at least two-fold higher or lower than the control, thus splitting transcripts in two groups: up-regulated and down-regulated.
Transcriptome functional analysis. GO terms for differentially expressed annotated transcripts were retrieved using a Blast2GO file provided by HanXRQ portal (https://www.heliagene.org/HanXRQ-SUNRISE/). GO enrichment analysis with Fisher's exact test on differential expressed genes versus whole predicted transcriptome was performed with Blast2GO analysis tools 125 using FDR corrected P value < 0.05. Overall enriched GO terms were summarized using REVIGO with allowed similarity of 0.4 126 .
GO slim (Blast2GO tool) was performed in order to reduce the complexity of GO terms for functional analysis of differentially expressed transcripts.
Scientific RePoRTS | (2018) 8:4 | DOI:10.1038/s41598-017-18445-0 qRT-PCR analysis. We used qRT-PCR to validate RNA-seq expression profiles of five genes selected from those differentially expressed between M4 and M16 samples. One gene encoding E3 ubiquitin-ligase BRE1-like 1 (genome ID HanXRQChr05g0138451 and one gene encoding mitochondrial ubiquitin ligase activator of NFKB 1 (genome ID HanXRQChr10g0282971) were selected among genes showing stable expression profiles in RNA-seq analyses and used as internal reference genes in qRT-PCR analyses. Complementary DNA from 400 ng total purified RNA was prepared using iScript cDNA synthesis kit (BioRad, Hercules, CA, USA) according to the manufacturer's protocol.
The synthesised cDNA was used for quantitative real-time polymerase chain reaction (qRT-PCR). Primers were designed using Primer3 software (Applied Biosystems) and are reported in Supplementary Table 2. qRT-PCR reactions (20 μL) were carried out with 10 ng of cDNA, 250 nM primers and 1X Fast SYBR Green Master Mix (Applied Biosystems, Foster City, CA, USA) following the manufacturer's instructions. PCR was run in a StepOne realtime PCR System (Applied Biosystems, Foster City, CA, USA) using the recommended thermal cycling conditions (hold 95 °C, 20 s; 40 cycles 95 °C, 3 s; 60 °C, 30 s).
The qRT-PCR results were obtained from three biological replicates and three technical repeats for each treatment. The relative abundance of transcripts was calculated by using the 2−ΔΔCt method 127 . For statistical analysis a two tailed t-test was used.