Identification of Populus Small RNAs Responsive to Mutualistic Interactions With Mycorrhizal Fungi, Laccaria bicolor and Rhizophagus irregularis.

Ecto- and endo-mycorrhizal colonization of Populus roots have a positive impact on the overall tree health and growth. A complete molecular understanding of these interactions will have important implications for increasing agricultural or forestry sustainability using plant:microbe-based strategies. These beneficial associations entail extensive morphological changes orchestrated by the genetic reprogramming in both organisms. In this study, we performed a comparative analysis of two Populus species (Populus deltoides and P. trichocarpa) that were colonized by either an arbuscular mycorrhizal fungus (AmF), Rhizophagus irregularis or an ectomycorrhizal fungus (EmF), Laccaria bicolor, to describe the small RNA (sRNA) landscape including small open reading frames (sORFs) and micro RNAs (miRNAs) involved in these mutualistic interactions. We identified differential expression of sRNAs that were, to a large extent, (1) within the genomic regions lacking annotated genes in the Populus genome and (2) distinct for each fungal interaction. These sRNAs may be a source of novel sORFs within a genome, and in this regard, we identified potential sORFs encoded by the sRNAs. We predicted a higher number of differentially-expressed miRNAs in P. trichocarpa (4 times more) than in P. deltoides (conserved and novel). In addition, 44 miRNAs were common in P. trichocarpa between the EmF and AmF treatments, and only 4 miRNAs were common in P. deltoides between the treatments. Root colonization by either fungus was more effective in P. trichocarpa than in P. deltoides, thus the relatively few differentially-expressed miRNAs predicted in P. deltoides might reflect the extent of the symbiosis. Finally, we predicted several genes targets for the plant miRNAs identified here, including potential fungal gene targets. Our findings shed light on additional molecular tiers with a role in Populus-fungal mutualistic associations and provides a set of potential molecular targets for future enhancement.

Introduction

Mycorrhizal fungi are a diverse group of beneficial organisms that colonize the roots of more than 90% of higher plant species. Typically, two main types of fungi are described: (1) arbuscular mycorrhizal fungi (AmF), developing inside the root cell walls, e.g., Rhizophagus irregularis, and (2) ectomycorrhizal fungi (EmF), colonizing the intercellular spaces of the roots, e.g., Laccaria bicolor (Bonfante and Genre, 2010). These fungi play an important role in the maintenance of the plant health and growth by promoting water cycling, nutrient exchange and enhanced tolerance/resistance to biotic and abiotic stresses, while in exchange, the fungi receive plant-fixed carbon (Smith and Read, 2008; Bonfante and Genre, 2010). Several studies have shown that field application of mycorrhizal fungi improves the overall productivity of a number of crops including cereals, legumes, fruits and trees (Abbott and Robson, 1977; Brundrett et al., 1996; Al-Karaki et al., 2004; Powell, 2018). To address the challenge to food and energy security caused by increases in the global population, and decreases in agricultural and forest land, it is important to gain a deeper understanding of the molecular mechanism underlying beneficial symbiosis between plant and fungi to effectively design and develop plant:microbe-based strategies to enhance forestry and agriculture health and sustainability (Martin et al., 2017).

Much progress has been made in understanding the establishment and maintenance of these mutualistic associations (Bonfante and Genre, 2010; Plett and Martin, 2011). Many studies support the hypothesis that fungi-derived protein signals or effectors facilitate and/or maintain the symbiotic interactions (Daguerre et al., 2017). For example, the genome of L. bicolor encodes a large number of mycorrhizal-induced small secreted proteins (MiSSPs), many of which are expressed and accumulated in the fungal hyphae during colonization (Martin et al., 2008). Plett et al. (2011) reported that the effector protein of L. bicolor, MiSSP7 could enter the nucleus of Populus roots cells to affect transcription and promote symbiosis. MiSSP7 protects the Populus jasmonate zim-domain protein 6 (PtJAZ6), which is a negative regulator of jasmonic acid (JA)-induced gene regulation in Populus, and promotes symbiosis by blocking the action of jasmonic acid (Plett et al., 2014). Similarly, a secreted protein of R. irregularis, SP7, interacts with ERF19, a pathogenesis-related transcription factor that counteracts the plant immune response and facilitates root colonization by this AmF (Kloppholz et al., 2011).

On the plant side, colonization by fungi requires that the plant distinguish between beneficial and pathogenic fungi (el Zahar Haichar et al., 2014). Plant root exudates have been shown to (1) mediate chemotaxis signaling that facilitate the colonization of roots by flagellated bacteria (Scharf et al., 2016) and (2) contain secondary metabolites, strigolactones, which promote hyphal branching and stimulates fungal spore germination (Akiyama et al., 2005). Moreover, transcriptome analysis of Populus trichocarpa roots colonized by L. bicolor has revealed 417 putative plant-encoded small secreted proteins (SSPs) with 39% of them appearing to be specific to Populus, where several of these SSPs were able to enter the hyphae and affect the growth of L. bicolor (Plett et al., 2017). These studies suggest that the genetic contributions from a plant in mutualistic association may be more complex than our current understanding and may involve several levels of regulation. It is unclear if this molecular toolbox for symbiosis, i.e., set of molecular determinants (e.g., protein-encoding genes, non-coding RNAs) are shared across different plant species when colonized by the same fungi or alternatively, the same plant species colonized by different types of symbiotic fungi.

In recent years, the role of small non-coding RNAs (sRNAs), broadly defined as regulatory RNA molecules ranging in size from 20 to 300 nucleotides (Großhans and Filipowicz, 2008), have become apparent in biotic stresses and regulation of plant development and physiology (Mallory and Vaucheret, 2006; Großhans and Filipowicz, 2008; Ruiz-Ferrer and Voinnet, 2009; Chen, 2012; Zhang and Chen, 2013). These regulatory RNA molecules include small interfering RNAs (siRNAs), microRNAs (miRNAs), piRNAs (Piwi-associated RNAs), and long non-coding RNAs (lncRNAs), which may originate from intergenic, intronic, or antisense transcripts. Several detailed reviews of molecular mechanism of these different population of non-coding RNAs (ncRNAs) were published recently (Ruiz-Ferrer and Voinnet, 2009; Chekanova, 2015; Mohanta and Bae, 2015; Huang et al., 2016). miRNAs, typically between 20 and 25 nucleotides, are processed from single-stranded RNA to form imperfect base-paired hairpin secondary structures, and generally negatively regulate their targets including mRNAs (Chen, 2008; Lanet et al., 2009) and ncRNAs such as TAS RNAs (Vaucheret, 2006). Many lines of evidence now confirm that miRNAs are necessary for plant association with AmF (Branscheid et al., 2010; Devers et al., 2011; Lauressergues et al., 2012; Etemadi et al., 2014). For example, colonization of Medicago truncatula roots by R. irregularis was reduced when miR171h was overexpressed (Lauressergues et al., 2012; Hofferek et al., 2014). In plants, sRNAs have also been shown to move between cells via the plasmodesmata, and transverse long distance via the phloem (Brosnan and Voinnet, 2011). Furthermore, sRNA can transverse into interacting fungi through the haustorium, to trigger gene silencing (Molnar et al., 2011; Huang et al., 2016). Consequently, sRNAs are now viewed as an important part of the genetic response network for evolutionary adaption (Amaral et al., 2013).

Studies have reported that sRNAs can also have protein-coding potential and encode small open reading frames (sORFs) which typically contain between 30 and 100 amino acids (Ulveling et al., 2011; Andrews and Rothnagel, 2014; Li et al., 2014; Ruiz-Orera et al., 2014; Kumari and Sampath, 2015). The portfolio of organisms with functional sORFs include bacteria, yeast, Drosophila, plants and human (Oyama et al., 2004; Kastenmayer et al., 2006; Hanada et al., 2007; Ladoukakis et al., 2011; Yang et al., 2011; Aspden et al., 2014; Ruiz-Orera et al., 2014; Shell et al., 2015; Ma et al., 2016). sORFs are partially classified based on their genomic position, e.g., upstream or downstream annotated genes, within annotated genes but out-of-frame, antisense to annotated genes, within novel transcripts and within ncRNAs. The biological significance for these types of sORFs varies. For example, the tarsel-less (tal) gene from Drosophila was classified as ncRNA but the transcript does encode four sORFs (11–32 AAs) that affect Drosophila gene expression and development (Galindo et al., 2007). The intergenic region of Arabidopsis has been predicted to encode approximately 33,809 sORFs (Hanada et al., 2007). Additionally, Castellana et al. (2008) identified ∼5,000 small peptides in Arabidopsis in a proteomic study, several of which were novel and/or identified in the aforementioned study by Hanada et al. (2007). Despite an unclear molecular mechanism, Hanada et al. (2013) showed varying morphological changes caused by overexpression of sORFs in Arabidopsis. Thus, it is plausible that sORFs derived from ncRNAs may provide another level of developmental regulation in mutualistic symbiosis.

In this study, we (1) identify and compare the Populus-derived sRNAs response to two different mutualistic symbionts (i.e., R. irregularis and L. bicolor) and (2) predict and compare the encoded sORFs, as well as their miRNAs and putative gene targets potentially involved in formation and maintenance of the beneficial associations. To achieve this goal, we profiled the sRNA response of two Populus species (P. deltoides and P. trichocarpa) that were each colonized by either of the two different fungi (Tagu et al., 2001; Labbé et al., 2011). Our study provides insight into the common and differential regulation that may be involved in both the endo- and ecto-mycorrhizal associations among members of the genus Populus.

Materials and Methods

Plant Material and Inoculation

Internode cutting of P. trichocarpa (genotype 101–74) and P. deltoides (genotype 73028-62) of the 54B F1 pedigree (Jorge et al., 2005) parental lines were rooted and individually inoculated with a 1/9 (v/v) mixture of L. bicolor inoculum and calcined clay (Oil Dri US Special, Damolin, Denmark). Eight replicates per genotype and treatment were carried out. L. bicolor S238N inoculum was produced by growing mycelium on an autoclaved peat-vermiculite nutrient mixture in 1-L glass jars for 2 months in the dark at 25°C and kept at 4°C before use (Duponnois and Garbaye, 1991). R. irregularis inoculum consisted of spores (200 spores per litter of calcinate clay). Inoculated cuttings were grown for 3 months post-inoculation in a greenhouse at 15–28°C, 12 h photoperiod at INRA-Nancy, France during summer and then harvested. One portion of the sampled root (about 25%) was used to assess successful root colonization, which was confirmed under a stereomicroscope by counting ectomycorrhizas for L. bicolor and for R. irregularis by counting arbuscules and measuring how much of the root length is colonized after Trypan blue staining procedure as described by Koske and Gemma (1989). The remaining lateral root samples (about 75%) were immediately frozen in liquid nitrogen at -80°C for RNA isolation.

RNA Isolation and Sequencing

Total RNA was extracted from the root system of three of the eight replicates per genotype per treatment by combining a CTAB extraction method and the SpectrumTM Total Plant RNA extraction kit (Sigma) as reported in Bryan et al. (2016). Prior to synthesizing sequencing libraries, RNA templates were quantified using a NanoDrop Spectrophotometer (Thermo Scientific) and RNA quality was determined using BioRad ExperionTM. Samples with a 260/280 reading between 2.0 and 2.2 and a RIN > 7 were determined to be of high enough quality and the high-quality RNA samples were used for sRNA-sequencing library preparation with an Illumina TruSeq small RNA sample prep kit (Illumina, CA, United States). As specified in the Illumina protocol, small RNA fragments were selected after PCR amplification by gel purification to capture RNA approximately 20–170 bp in length. The Illumina TruSeq small RNA library prep workflow is considered to be stranded because adapters are ligated directionally, facilitating maintenance of strand information during analysis. A total of 18 sRNA-sequencing libraries (three biological replicates of the control P. deltoides and control P. trichocarpa, and from the treatments, EmF-P. deltoides, EmF- P. trichocarpa, and AmF-P. deltoides and AmF-P. trichocarpa) were constructed. The libraries were sequenced on an Illumina MiSeq system (Illumina, CA, United States) to generate singe-end reads (1 × 150 bp).

sRNA-Seq Data Analysis

In this study, small RNAs are defined as the RNA molecules with size ranging from 20 to 300 nucleotides, consistent with the definition by Großhans and Filipowicz (2008). The raw RNA-Seq data in Fastq format were quality checked and trimmed for low-quality regions and adaptor sequences using Trimmomatic v0.35 (Bolger et al., 2014), parameters of “2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:20.” The trimmed sequencing reads from P. trichocarpa and P. deltoides were mapped to the P. trichocarpa genome v3.0 and P. deltoides WV94 v2.1[1], respectively, using TopHat2 (Kim et al., 2013). Novel genes were assembled and transcript abundance in RPKM (Reads Per Kilobase of transcript per Million mapped reads) was estimated using Cufflinks (Trapnell et al., 2012). Analysis of differential gene expression was performed using Cuffdiff (Trapnell et al., 2012). The fold-changes of differential expression was calculated by log2(RPKM ratio), with a pseudo-count of 1 added to each RPKM value. To identify the homologs to the significantly differentially expressed transcripts, the sequences were used as query for BLASTN search against the publicly available genomes in the Phytozome[2] using an E-value cutoff of 1.0E-3.

sORF Annotation

ORFfinder (Wheeler et al., 2003) was used to annotate putative ORFs within the differentially expressed sRNAs. We searched for 6 translational frames with the standard genetic code between 30 and 300 nucleotides for each prediction and selected the longest full-length sORF (including start and stop codon). The subcellular localization was predicted for each sORF using Loctree3 (Goldberg et al., 2014).

miRNA Identification

The small RNA sequencing reads were filtered to remove adaptors and low-quality bases by FastQC (Andrews, 2010). The sequencing reads were processed to remove adaptors and cleaned by Q30 value. Reads with over 20% bases less than Q30, and N base more than 10% were filtered. The clean reads were further selected by length between 18 and 30 bp for miRNA identification process. The miRDP1.3 package (Yang and Li, 2011) was used for miRNA precursor prediction with size set at 250 nt. To validate the miRNAs, the secondary structure and reads distribution were evaluated to generate the final set of miRNAs (Yin et al., 2016). The precursors and mature sequences were aligned with the searching toolbox in PNRD database for annotation of miRNAs (Yi et al., 2015). To define the conservation of miRNAs, a 17 bp match in the mature sequences were used to identify conserved miRNA families, allowing two mismatches. To quantify the abundance of miRNA, transcripts per million (TPM) value was defined as “counts of read mapped to miRNA∗1,000,000”/”reads mapped to reference genome” (Fahlgren et al., 2007). To identify differentially expressed miRNA DESeq2 software was used, and expression folder change greater than 2 and p-value (Benjamini-Hochberg FDR corrected) less than 0.01 were selected (Love et al., 2014). The miRNA expression data were normalized by row z-score method. Hierarchical clustering of gene expression was performed by clustergram function in MATLAB Bioinformatics toolbox with minor changes.

Target Prediction and Degradome Analysis

The miRNA targets were annotated by standard settings of psRNATarget (Dai and Zhao, 2011) with maximum expectation value 2.0. The transcripts of P. deltoides (WV94_445 v2.1) and P. trichocarpa (Nisqually-1 v3.0) were downloaded from Phytozome 12[3] while the transcripts for L. bicolor (20110203) and R. irregularis (Gloin1_GeneCatalog_transcripts_20120510.nt) was downloaded from MycoCosm[4]. To validate the targets, the degradome datasets in P. trichocarpa were downloaded from NCBI Short Read Archive (SRR4010497, SRR4010498). The clean reads from degradome library were obtained for target site identification. The CleaveLand4.0 pipeline (Addo-Quaye et al., 2009) was used for target scanning. The reads were mapped to transcript datasets and the alignment scores and p-value were calculated according the signatures (abundances of potential slicing end based on reads distributions). The t-plots were generated to visualize the miRNA directed slicing to targets and p-value of less than 0.05 was used to identify highly confident targets. Gene ontology enrichment for biological processes and molecular functions were performed online at PopGenIE[5] and visualized with REVIGO using default settings (Sjödin et al., 2009; Supek et al., 2011).

Results

Differential Expression of sRNAs During Populus-AmF/EmF Colonization

The AmF, R. irregularis, and the EmF, L. bicolor, were able to colonize both species of Populus, i.e., P. deltoides and P. trichocarpa, though colonization rates were higher for P. trichocarpa (Table 1). To identify differentially expressed sRNAs in response to the AmF and EmF treatments, we performed RNA-sequencing (RNA-Seq) of RNA fragments less than 300 nt isolated from roots of P. deltoides colonized with the AmF R. irregularis or the EmF L. bicolor (PDA and PDE, respectively) and from roots of P. trichocarpa colonized with R. irregularis or L. bicolor (PTA and PTE, respectively, Figure 1A). Roots without fungal inoculation from P. deltoides and P. trichocarpa were used as controls (PTC and PDC, respectively). Thus, a total of 18 sRNA-sequencing libraries (three biological replicates per treatment PDA, PDE, PTA, PTE, PDC, and PTC) were constructed and sequenced. Following trimming, the RNA-Seq reads were mapped to the respective P. trichocarpa v3.0 and P. deltoides WV94 v2.1 genome sequences (Table 2). Based on a corrected p ≤ 0.05 and > 2-fold change in expression compared to the controls, we found 81 and 59 differentially expressed transcripts in P. deltoides and P. trichocarpa, respectively, relative to the control (Supplementary Tables S1, S2 and Supplementary Files S1, S2). Among the 81 fungal-responsive transcripts in P. deltoides, 12 and 69 transcripts were differentially expressed in the PDA and PDE interactions, respectively. Furthermore, among the 59 fungal-responsive transcripts in P. trichocarpa, 32 and 27 transcripts were differentially expressed in the PTA and PTE interactions, respectively. Interestingly, the differential transcripts were mostly specific to each plant-fungal combination, except for 7 transcripts that were shared between the PDA and PDE (all down-regulated except Podel.CUFF.90 which was up-regulated) and 3 transcripts shared between PTA and PTE (two down-regulated while Potri.013G036500 showed an opposite trend in each treatment; Figure 1B). The expression fold-change trend is shown in Figure 1C and Supplementary Tables S1, S2. Across species comparison showed that 5 homologous transcripts could be detected between species, of which three homologous transcript pairs showed the opposite expression trend between PDE and PTA and the other two homologous transcripts pairs were up-regulated in PDE and PTE (Figure 1D). These results suggest that, to a large extent, the expression of the sRNAs were distinct for each plant species within corresponding fungal treatments despite both fungi having a mutualistic relationship with these plants.

Table 1

Frequency of root colonization by arbuscular (AmF) and ectomycorrhizal (EmF) mycorrhizal fungus in Populus deltoides (PD) and Populus trichocarpa (PT).

Treatments	AmF colonization (%)		EmF colonization (%)
	PD	PT	PD	PT
Control	0	0	0	0
Inoculated	72a (+/-5)	87a (+/-5)	34b (+/-5)	57b (+/-5)

Figure 1

sRNA landscape of Populus species in response to the arbuscular mycorrhizal fungus (AmF), Rhizophagus irregularis and the ectomycorrhizal fungus (EmF), Laccaria bicolor. (A) Populus species and corresponding treatments used in the current study; (B) Venn diagram showing the transcript overlap from the different treatments in P. deltoides and P. trichocarpa; (C) Significant differentially expressed transcripts from the controls (p ≤ 0.05) and greater/less than two-fold change in the expression level; and (D) comparison of fold-change and trend for conserved transcripts across species. Transcript IDs and fold change can be found in Supplementary Tables S1, S2 and transcript sequences can be found in Supplementary Files S1, S2.

Table 2

Total number of RNA-Seq reads (in the range of three biological replicates per genotype/treatment), following trimming, that were mapped to Populus deltoides WV94 v2.1 and Populus trichocarpa genome v3.0 available on Phytozome 12.

Tissue	Treatment	Total reads	Reads mapped to respective genome	Percent mapped (%)
P. deltoides	Control-1	5749345	2611186	45.4
P. deltoides	Control-2	3014126	1940271	64.4
P. deltoides	Control-3	3075042	2066300	67.2
P. deltoides	Rhizophagus irregularis-1	4980677	2539132	51.0
P. deltoides	Rhizophagus irregularis-2	2993283	1524495	50.9
P. deltoides	Rhizophagus irregularis-3	2206148	1244202	56.4
P. deltoides	Laccaria bicolor-1	3564011	666296	18.7
P. deltoides	Laccaria bicolor-2	2080970	712191	34.2
P. deltoides	Laccaria bicolor-3	2232290	603916	27.1
P. trichocarpa	Control-1	2612622	1478334	56.6
P. trichocarpa	Control-2	1928877	982300	50.9
P. trichocarpa	Control-3	2744908	1280051	46.6
P. trichocarpa	R. irregularis-1	4234040	2492972	58.9
P. trichocarpa	R. irregularis-2	3012095	1482320	49.2
P. trichocarpa	R. irregularis-3	1968986	775833	39.4
P. trichocarpa	L. bicolor-1	4053717	781250	19.3
P. trichocarpa	L. bicolor-2	2063636	516891	25.0
P. trichocarpa	L. bicolor-3	1486150	268589	18.1

Apart from the 17 transcripts that were previously annotated within the respective genomes (Supplementary Table S3), the remaining significantly differentially expressed fungus-responsive plant sRNAs (with name starting as “CUFF”) are located in the genomic regions without previously annotated genes. From the 17 previously annotated genes, 6 were proteins of known function (e.g., MYB-LIKE DNA-BINDING PROTEIN and AMMONIUM TRANSPORTER 1 MEMBER), while the remaining 11 genes were annotated as proteins of unknown function (Supplementary Table S3). We assessed the conservation of the differentially expressed sRNA across 63 plant genomes available on Phytozome 12.0 and found homology to 47 P. deltoides and 20 P. trichocarpa sRNAs in at least one genome analyzed (E-value ≤ 03; Supplementary Table S4).

Novel sORFs in Populus in Response to Mycorrhizal Fungi

Since sRNAs can potentially encode sORFs (Li et al., 2014; Ruiz-Orera et al., 2014), we used ORFfinder (Wheeler et al., 2003) to scan the differentially expressed sRNAs for potential sORF translations. In total, we could predict 22 and 6 full-length sORFs from the treatments in P. deltoides and P. trichocarpa, respectively (Supplementary Tables S5, S6). The shortest predicted sORF was 39 bp (12 AA) while the longest was 273 bp (90 AA). The transcript, coding, and protein sequences of these sORFs can be found in Supplementary Files S1, S2, S5, S6. Prediction of subcellular localization showed that 26 sORFs are secreted extracellularly while a single sORF is targeted to the nucleus and another to the chloroplast (Supplementary Table S7).

Identification of miRNAs in Populus in Response to Mycorrhizal Fungi

To gain insight into the role of miRNAs in mutualistic interactions, we analyzed the sRNA datasets to identify miRNAs that were responsive to the different mycorrhizal treatments. Using the miRDeep-P pipeline (Yang and Li, 2011), we identified 287 putative miRNAs in P. deltoides and 357 putative miRNAs in P. trichocarpa, including precursor and mature sequences (Supplementary Files S3, S4 and Supplementary Tables S8, S9). In P. deltoides and P. trichocarpa, 174 and 130 miRNAs were classified into “Plant conserved” miRNA families, respectively, and 113 and 227 miRNAs were classified as “Plant novel” miRNAs, respectively (Figure 2A and Supplementary Tables S8, S9). Note that the corresponding known families for the miRNAs identified in this study can found in Supplementary Tables S8, S9 while the miRNA abundance expressed as TPM for each treatment can be found in Supplementary Tables S10, S11. To refine the list of miRNAs, we focused on those that were significantly differentially expressed (negative binomial test, FDR corrected p ≤ 0.05; fold-change > 2; Soneson and Delorenzi, 2013) between the treatment and control. In this regard, we only considered miRNAs with greater than 50 counts in total and found 34 and 130 unique miRNAs with significant differential expression in P. deltoides and P. trichocarpa, respectively, including Plant Novel miRNAs (Supplementary Tables S12, S13). Hierarchical clustering of differentially expressed miRNAs showed distinct differences in expression pattern in response to EmF and AmF in both P. deltoides and P. trichocarpa (Figure 2B).

Figure 2

miRNA landscape Populus deltoides and P. trichocarpa roots colonized with Rhizophagus irregularis (PDA and PTA, respectively) and Laccaria bicolor (PDE and PTE, respectively). (A) miRNAs identified were classified as “Plant conserved” or “Plant novel”; (B) Hierarchical clustering of differential expression in transcripts per million of identified miRNAs; (C) Venn showing significant differentially expressed miRNAs that were common and unique miRNAs to each treatment; and (D) Fold-change and the significance (p ≤ 0.05) of differentially expressed miRNAs.

We next investigated whether each species shared a common set of miRNAs during the different mutualistic symbiosis. We found that in P. deltoides 4 miRNAs were common between PDA and PDE and in P. trichocarpa 44 miRNAs were common between PTA and PTE (Figure 2C and Supplementary Tables S12, S13). In P. deltoides, all 4 common miRNAs were down-regulated relative to the control while in P. trichocarpa, 16 shared miRNAs were up-regulated, 23 were down-regulated and 5 showed the opposite fold-change trend (Table 3). Further, this analysis also revealed that, within species, 3 and 27 miRNAs were unique to PDA and PDE, respectively, while 78 and 8 miRNAs were unique to PTA and PTE, respectively (Figure 2C and Supplementary Tables S12, S13). We further compared the Best-Hit of the significantly differentially expressed miRNAs across species, i.e., PDA vs. PTA and PDE vs. PTE (Supplementary Tables S12, S13). We did not find any miRNAs in common between different species colonized by a common fungus (data not shown). Also, we compared the differential transcripts listed in DATA 1 and 2 with the miRNA precursor sequences using the CAP3 program implemented in http://biosrv.cab.unina.it/webcap3/ and identified that five differential transcripts were miRNA precursors (Supplementary Table S14). The miRNA expression fold-change trend can be seen in Figure 2D and Supplementary Tables S12, S13. In P. trichocarpa, we found that 23 and 7 miRNAs were only detected in the PTA and PTE treatments, respectively, and not in the controls.

Table 3

Expression fold-change trend of miRNAs in Populus deltoides and P. trichocarpa in response to Rhizophagus irregularis (PDA and PTA, respectively) and Laccaria bicolor (PTA and PTE, respectively).

		Down-regulated miRNAs in P. deltoides
miRNA ID	miRNA name	PDA	PDE
Pde_miRNA_219	pde-MIRf11885b	-1.57	-5.28
Pde_miRNA_228	pde-MIR393f	-1.58	-5.6
Pde_miRNA_155	pde-MIR396a	-1.47	-5.08
Pde_miRNA_156	pde-MIRf12503b	-1.43	-4.93
		Up-regulated miRNAs in P. trichocarpa
		PTA	PTE
Ptr_miRNA_178	ptr-miR482a	7.6	5.31
Ptr_miRNA_277	ptr-miR164c	5.29	12.95
Ptr_miRNA_299	ptr-miR166g	14.54	9.39
Ptr_miRNA_248	ptr-miR162a	5.06	7.8
Ptr_miRNA_206	ptr-miR167g	6.2	7.77
Ptr_miRNA_62	ptr-miR2912x	3.58	4.46
Ptr_miRNA_179	ptr-miR167a	5.19	8.56
Ptr_miRNA_298	ptr-miR159b	10.31	5.83
Ptr_miRNA_15	ptr-miR812a	4.43	12.96
Ptr_miRNA_342	ptr-miR156h	Not detected in control, induced in PTA	Not detected in control, induced in PTE
Ptr_miRNA_339	ptr-miR171i	Not detected in control, induced in PTA	Not detected in control, induced in PTE
Ptr_miRNA_330	ptr-miR482a	Not detected in control, induced in PTA	Not detected in control, induced in PTE
Ptr_miRNA_282	ptr-miR319a	Not detected in control, induced in PTA	Not detected in control, induced in PTE
Ptr_miRNA_247	ptr-miR159a	Not detected in control, induced in PTA	Not detected in control, induced in PTE
Ptr_miRNA_113	ptr-miR166d	Not detected in control, induced in PTA	Not detected in control, induced in PTE
Ptr_miRNA_131	ptr-miR156g	Not detected in control, induced in PTA	Not detected in control, induced in PTE
		Down-regulated miRNAs in P. trichocarpa
		PTA	PTE
Ptr_miRNA_155	ptr-miR2912g	-4.65	-4.84
Ptr_miRNA_473	ptr-miR2912w	-1.58	-2.67
Ptr_miRNA_44	ptr-miR2912v	-6.92	-10.09
Ptr_miRNA_160	Ptr_miRNA_160	-12.58	-12.07
Ptr_miRNA_210	ptr-miR2912m	-3.37	-4.3
Ptr_miRNA_311	ptr-miR2111b	-1.19	-6
Ptr_miRNA_357	ptr-miR156j	-4.59	-7.08
Ptr_miRNA_16	ptr-miR2912i	-8.4	-15.39
Ptr_miRNA_255	ptr-miR2912n	-1.08	-4.25
Ptr_miRNA_209	Ptr_miRNA_209	-3.99	-6.1
Ptr_miRNA_105	ptr-miR482d	-6.09	-1.43
Ptr_miRNA_328	ptr-miR839	-9.57	-11.79
Ptr_miRNA_329	ptr-miR156d	-5.59	-11.13
Ptr_miRNA_135	ptr-miRf11683	-1.58	-3.39
Ptr_miRNA_29	ptr-miR6427	-1.54	-3.14
Ptr_miRNA_110	ptr-miR2912b	-2.31	-6.18
Ptr_miRNA_156	ptr-miR2912h	-4.33	-6.07
Ptr_miRNA_369	ptr-miR476	-0.82	-1.56
Ptr_miRNA_370	ptr-miR2912s	-3.26	-3.4
Ptr_miRNA_187	ptr-miR172b	-2.15	-3.95
Ptr_miRNA_106	Ptr_miRNA_106	-6.07	-7.17
Ptr_miRNA_21	Ptr_miRNA_21	-2.01	-5.36
Ptr_miRNA_232	Ptr_miRNA_232	-5.29	-8.41
		miRNAs showing opposite expression trend in P. trichocarpa treatments
		PTA	PTE
Ptr_miRNA_191	ptr-miR475b	-2.05	7.58
Ptr_miRNA_134	ptr-miR2912e	-0.81	5.14
Ptr_miRNA_275	ptr-miR166f	11.17	-2.37
Ptr_miRNA_103	ptr-miR166b	10.58	-1.38
Ptr_miRNA_337	ptr-miR396c	5.18	-4.97

Identification of Target Genes of miRNAs in Populus in Response to Mycorrhizal Fungi

To determine the biological function of all identified miRNAs, we searched the “degradome” libraries of public databases for potential gene targets in Populus. While we identified 61 gene targets to 41 unique miRNAs in P. deltoides and 117 gene targets to 46 unique miRNAs in P. trichocarpa (p ≤ 0.05; Supplementary Tables S15, S16), there were only 12 and 59 unique gene targets (including alleles) to 7 and 25 unique differentially expressed P. deltoides and P. trichocarpa miRNAs, respectively (Supplementary Table S17). We also identified putative targets to the “Plant novel” miRNAs category, i.e., 10 miRNA-gene targets in P. deltoides and 18 miRNA-gene targets in P. trichocarpa (Table 4). In both P. deltoides and P. trichocarpa, gene ontology enrichment analysis of the miRNA targets (Supplementary Tables S15, S16) for biological processes included terms associated with regulation of metabolism, response to hormones and endogenous stimulus and regulation of transcription (Figure 3A,B) while molecular functions in P. trichocarpa were enriched in transcription factor activity, nucleic acid binding and DNA binding (Figure 3C). There were no enrichments for molecular functions in P. deltoides.

Table 4

Gene targets predicted for the “plant novel” Populus miRNAs in response to Rhizophagus irregularis and Laccaria bicolor.

miRNA ID	miRNA name	Target gene	Annotation of target gene
Pde_miRNA_73	pde-MIR2628b	Podel.02G280600.1	Photosystem I subunit IV (psaE)
Pde_miRNA_85	pde-MIR857	Podel.03G089000.1	Beta-amylase/Saccharogen amylase
Pde_miRNA_115	pde-MIRf12523d	Podel.03G183100.1	CALMODULIN-BINDING TRANSCRIPTION ACTIVATOR
Pde_miRNA_21	pde-MIR395a	Podel.09G062800.1	ZINC FINGER CCCH DOMAIN-CONTAINING PROTEIN 56
Pde_miRNA_187	pde-MIRf10448b	Podel.07G012300.1	MATE EFFLUX FAMILY PROTEIN
Pde_miRNA_133	pde-MIR2912h	Podel.04G165900.1	FRUCTOSE-BISPHOSPHATE ALDOLASE 2, CHLOROPLASTIC-RELATED
Pde_miRNA_149	pde-MIRf12506i	Podel.01G087300.1	Mediator of RNA polymerase II transcription subunit 10
Pde_miRNA_225	pde-MIRf10282	Podel.19G057400.1	D-XYLOSE-PROTON SYMPORTER-LIKE 1-RELATED
Pde_miRNA_135	pde-MIR5715a	Podel.09G129400.1	FRUCTOSE-BISPHOSPHATE ALDOLASE 2, CHLOROPLASTIC-RELATED
Pde_miRNA_53	pde-MIRf10816c	Podel.11G131500.1	POLYUBIQUITIN 3
Ptr_miRNA_169	ptr-miR166a	Potri.009G014500.2	Homeobox-leucine zipper family protein/lipid-binding START domain-containing protein
Ptr_miRNA_221	ptr-miRf10117a	Potri.005G091700.1	Serine carboxypeptidase-like 34
Ptr_miRNA_283	ptr-miR847	Potri.014G071000.1	Protein of unknown function
Ptr_miRNA_269	ptr-miRf10610b	Potri.002G154600.2	Ribosomal L5P family protein
Ptr_miRNA_105	ptr-miR482d	Potri.018G138500.1	NB-ARC domain-containing disease resistance protein
Ptr_miRNA_335	ptr-miR7748a	Potri.008G181000.5	Response regulator 11
Ptr_miRNA_409	ptr-miRf10535	Potri.001G038100.1	Serine-domain containing serine and sphingolipid biosynthesis protein
Ptr_miRNA_223	ptr-miRf10117d	Potri.001G234400.1	Nudix hydrolase homolog 9
Ptr_miRNA_85	ptr-miRf10350	Potri.007G013100.1	Ribosomal protein L2 family
Ptr_miRNA_215	ptr-miRf10010b	Potri.002G124200.5	Poly(A) binding protein 2
Ptr_miRNA_467	Ptr_miRNA_467	Potri.006G138600.1	Chaperonin 20
Ptr_miRNA_259	ptr-miRf10010c	Potri.002G124200.3	Poly(A) binding protein 2
Ptr_miRNA_71	ptr-miRf10010d	Potri.002G124200.5	Poly(A) binding protein 2
Ptr_miRNA_161	ptr-miRf10116a	Potri.005G091700.3	Serine carboxypeptidase-like 34
Ptr_miRNA_421	Ptr_miRNA_421	Potri.005G192300.1	Tubby-like protein 5
Ptr_miRNA_417	Ptr_miRNA_417	Potri.001G189700.1	NADP-malic enzyme 4
Ptr_miRNA_51	Ptr_miRNA_51	Potri.001G001500.1	RING/U-box superfamily protein
Ptr_miRNA_285	ptr-miR169h	Potri.010G218700.2	GTP binding Elongation factor 1 alpha family protein

Figure 3

Gene ontology enrichment of the target genes of the miRNAs identified in Populus deltoides and P. trichocarpa roots colonized with the arbuscular mycorrhizal fungus (AmF), Rhizophagus irregularis, and the ectomycorrhizal fungus (EmF), Laccaria bicolor. (A) Gene ontology biological processes for P. deltoides with AmF and EmF; (B) Gene ontology biological processes for P. trichocarpa with AmF and EmF; and (C) Gene ontology molecular processes for P. trichocarpa with AmF and EmF. Enrichment was performed online at PopGenIE (http://popgenie.org/) and visualized with REVIGO (Sjödin et al., 2009; Supek et al., 2011).

Identification of Putative Target Genes for Populus-Derived miRNAs in R. irregularis and L. bicolor

We predicted potential miRNA-gene targets in the respective fungi since small RNA may transverse between symbionts via extracellular vesicles to mediate RNA interference and effect physiological changes across species (Lefebvre and Lécuyer, 2017). We found 5 gene targets in L. bicolor for 14 significantly differentially expressed P. deltoides miRNAs, 4 gene targets in L. bicolor for 6 significantly differentially expressed P. trichocarpa miRNAs and 5 gene targets in R. irregularis for 5 significantly differentially expressed P. trichocarpa miRNAs (Supplementary Table S18). We did not find any genes targets in R. irregularis to the significantly differentially expressed P. deltoides miRNAs. The majority of the predicted gene targets within the respective fungi were hypothetical protein/predicted proteins while a transport protein particle complex subunit (TRAPP 20 K subunit) was a common target for several P. deltoides and P. trichocarpa miRNAs. In addition, despite an annotation as a hypothetical protein, the target genes of two miRNAs had a KOG description as transcription factors.

Discussion

In the current study, we applied high-throughput RNA sequencing to describe the sRNA landscape of the root of two Populus species (P. deltoides and P. trichocarpa) colonized by an EmF, L. bicolor, or an AmF, R. irregularis, with specific focus on the sORF and miRNA response. A striking observation was that these beneficial plant-fungal interactions triggered the expression of otherwise intergenic regions of the genome. Furthermore, despite being mutualistic symbiosis, only a small number of these transcripts were conserved within the species between the two fungi (7 in P. deltoides and 3 in P. trichocarpa; Figure 1B). This may be a reflection of the different classes of fungi, i.e., arbuscular mycorrhizal vs. ectomycorrhizal, suggesting that at least on the plant side and for the sRNAs, mutualistic associations do not necessarily have a common response/regulatory mechanism. Additionally, only 5 transcripts were shared between P. deltoides and P. trichocarpa (Figure 1D), leading to a hypothesis that the molecular differences (i.e., the low number of sRNAs shared between two Populus species under two different fungal treatments) in the expression of plant small RNAs in response to fungi might contribute to the phenotypic difference in colonization efficiency (see Table 1; Tagu et al., 2001). This hypothesis needs to be tested using experimental approaches such as gain-of-function via gene overexpression and loss-of-function via gene knockout in the future.

While we do not know for certain what the end products of each sRNA is, several lines of evidence exist for the functionality of sORFs translated from sRNAs (Andrews and Rothnagel, 2014; Ruiz-Orera et al., 2014). Given that non-coding RNA may be an important source of new peptides, our predicted coding potential of the sRNAs advocates for a functional role for several sORFs in the symbiotic interactions. This is further supported by the predicted subcellular targeting of these sORFs (Supplementary Table S7). Recently, Plett et al. (2017) showed that several Populus-encoded SSPs (≤ 250 AA), including a novel SSP that was not annotated in the P. trichocarpa genome, could enter the fungal hyphae, localize to the nucleus and affect growth of L. bicolor.

The differential expression of a large number of miRNAs identified in response to the fungal treatments (34 in P. deltoides and 130 P. trichocarpa) suggest that the plant miRNAs are an important regulator of this mutualistic symbiosis. In this respect, our differentially expressed miRNA list (Supplementary Tables S12, S13) contained previously identified miRNAs that are previously reported to be involved in mutualistic interactions, albeit in different plant species such as Medicago, tomato and rice, including miR156, miR160, miR170, miR167, miR393, and miR396 (Bazin et al., 2013; Etemadi et al., 2014; Wu et al., 2016). Unlike the sRNA landscape, several significantly differentially expressed miRNAs were shared between PTA and PTE treatments (Figure 2C), indicating that the miRNA landscape and hence, the molecular mechanism imparted by this level of regulation may be conserved to a certain extent within species. The latter is further substantiated by the same expression trend seen for 39 of the 44 miRNAs in common between PTA and PTE and all 4 miRNAs in common between PDA and PDE, where both miRNAs were either up- or down-regulated in the respective treatments (highlighted in blue in Supplementary Tables S12, S3). We however, also noted that the expression pattern of miRNAs belonging to the same family did also vary within a treatment. For example, miR156 family targeting members of the squamosa promoter binding and binding-like transcription factors were up- and down-regulated in PTA (Ptr_miRNA_245 [ptr-miR156b], Ptr_miRNA_329 [ptr-miR156d], Ptr_miRNA_127 [ptr-miR156e], Ptr_miRNA_131 [ptr-miR156g], and Ptr_miRNA_433 [ptr-miR156k]; Supplementary Table S13). This suggests a very intricate role for members of the same miRNA family in these plant-microbe interactions, and further, indicates the specificity to which the miRNAs functions in these tested relationships.

Given that the number of miRNAs identified in P. trichocarpa was approximately 4 times higher than in P. deltoides, and that a large number of miRNAs were in common between PTA and PTE, many with the same expression trend, we hypothesize that species-specific responsive miRNAs may play an important role in the preferential colonization of P. trichocarpa over P. deltoides, as previously reported (Tagu et al., 2001; Labbé et al., 2011). For example, consistent with previous findings in Solanum, Medicago, and Oryza, miR393, which negatively regulates arbuscule formation via an obstruction to the auxin perception, was also down-regulated in PTA in our study (Ptr_miRNA_289 [ptr-miR393b]) but was absent in PDA. The predicted target for Ptr_miRNA_289 [ptr-miR393b] was Potri.002G207800.1, annotated as Transport Inhibitor Response I (TIR/F-box/RNI-like superfamily protein; Supplementary Table S17). When auxin is perceived, TIR of the E3 ligase receptor complex targets the Aux/IAA protein for degradation via ubiquitination which releases the Auxin Response Factors (ARFs) to allow auxin induced gene expression of their targets (Sukumar et al., 2013; Etemadi et al., 2014). Gene components of this auxin signaling are associated with root development where auxin typically stimulates lateral root development and inhibits root elongation (Casimiro et al., 2003; Sukumar et al., 2013). Interestingly, in our study, the two ARFs targets identified, Potri.002G055000.2 and Potri.011G091900.4, should have increased expression since their corresponding miRNA, Ptr_miRNA_407 (ptr-miR167b) was down-regulated as well; thereby supporting the auxin- perception and signaling previously described (Sukumar et al., 2013).

It is also possible that AmF colonization of the P. deltoides roots, which still occurs albeit at a low frequency, proceeds via alternative means. Lauressergues et al. (2012) showed that there was an up-regulation of miR171h during the colonization of Medicago truncatula roots by R. irregularis. The gene target of miR171h is a GRAS transcription factor, NSP2, that is believed to be involved in the biosynthesis of a hormone, strigolactone which promotes germination and growth of AmF. Interestingly, the overexpression of miR171h is associated with a decrease in the target gene expression and an overall decrease in the fungal colonization; while the overexpression of a miR171h-resistant NSP2 gene showed over-colonization. These results suggest that that miR171h is a negative regulator of NSP2 thereby preventing over-colonization. Another study reported increased NSP2 gene levels in mycorrhizal roots despite the presence of miR171h; therefore, suggesting the presence of an addition component regulating miR171h-NSP2 (Hofferek et al., 2014). In our study, we noted an up-regulation of miR171 in PTA and PTE (Ptr_miRNA_339 [ptr-miR171i], common between the treatments; Supplementary Table S13) and down-regulation of miR171 in PDA and PDE (Pde_miRNA_228 [pde-MIR393f], also common between the treatments; Supplementary Table S12). According to Lauressergues et al. (2012), if miR171 is down-regulated, as in the case of PDA, then NSP expression will not be down-regulated thereby increasing mycorrhizal colonization. The alternative proposed by Hofferek et al. (2014) is also true where the NSP2 transcript was elevated during the progression of mycorrhizal colonization independent of miR171h levels, and in the case of PDA will once again promote increased colonization. This hypothesized route may assist in the symbiotic interaction in P. deltoides which has a lower propensity for colonization. A similar trend was noted with pathogenic fungi where a specific miRNAs were either repressed or induced in different plant species, suggesting functional differences in these alternant plants to common pathogens (Gupta et al., 2014). Moreover, as AmF symbiosis progresses miR171 levels should increase (Lauressergues et al., 2012; Hofferek et al., 2014). Ptr_miRNA_339 (ptr-miR171i) was upregulated in the PTA interactions. If miR171 is a negative regulator of NSP2 thereby preventing over-colonization (Lauressergues et al., 2012), the increase in expression in miR171 in PDA could serve as a checkpoint to prevent an imbalance in the mutualism between P. trichocarpa, with its high propensity to be colonized, and the fungus. Future experimental characterization and profiling of these responsive miRNAs and targets should shed some light on this method of colonization and control.

We were also able to identify putative gene targets to the significantly differentially expressed conserved and novel Populus miRNAs (Supplementary Table S17). Our gene ontology enrichment analysis showed consistency with a previous genome-wide study investigating the miRNAs responsive to arbuscular mycorrhiza in tomato, including transcription regulation, biological regulation and response to stimulus (Figure 3; Wu et al., 2016). It was interesting to note that the molecular functions of P. trichocarpa in mycorrhizal associations included transcription factor activity, nucleic acid binding and DNA binding. Plett et al. (2017) showed that SSPs could enter the hyphal nucleus of L. bicolor to affect the fungal physiology. Closer inspection of the miRNAs and their corresponding gene target revealed a large number of transcription factors associated with root formation including growth-regulating factors (GRFs), homeodomain-leucine zipper transcription factor, AP2 transcription factor and MYB transcription factors (Supplementary Table S17). We also found phytohormone related genes including ethylene-responsive transcription factor and auxin response factor, as noted in other studies.

In plants, sRNAs move to adjacent cells via the plasmodesmata and systemically via the vasculature system (Baulcombe, 2004). Recently, it was shown that cross-kingdom/organism RNA interference occurs via exosome-like extracellular vesicles that transfer the host sRNAs to the interacting partner at the infection site, which then consequently causes gene silencing (Wang et al., 2017; Cai et al., 2018; Shahid et al., 2018). Therefore, it is highly probable that plant-encoded sRNAs, during mycorrhizal symbiosis, extends beyond the self-regulation. In support of this hypothesis, we used the plant-encoded miRNAs to predict fungal gene targets in the respective genomes. We identified several fungal gene targets to the Populus-derived miRNAs, including transport proteins, transcription factors and several genes encoding proteins of unknown function in the two fungi analyzed (Supplementary Table S18).

Finally, our study adds to the growing number of reports suggesting the active role of plants in mutualistic associations with microbes (Akiyama et al., 2005; Scharf et al., 2016; Plett et al., 2017). While a number of studies have focused on the arbuscular mycorrhizal-responsive miRNAs (Devers et al., 2011; Wu et al., 2016), genome-wide studies on the endomycorrhizal-responsive miRNAs is in its infancy. Given the role of miRNAs in fungal interactions (Bazin et al., 2013; Etemadi et al., 2014) and the expanding evidence for sORFs in plant physiology and mutualistic symbiosis (Hanada et al., 2007, 2013; Andrews and Rothnagel, 2014; Plett et al., 2017), the current study provides a useful resource of potential regulators that are involved in arbuscular/endomycorrhizal symbiosis with Populus.

Data Availability

The Department of Energy (DOE) will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan). The RNA-Seq reads, along with the processed data (including the sequence and abundance measurements of differentially-expressed transcripts) are deposited in NCBI GEO (https://www.ncbi.nlm.nih.gov/geo/) with the accession code GSE117158.

Author Contributions

XY and JL conceived and supervised the project. RM and HY carried out bioinformatic analysis, interpreted the data, and wrote the manuscript. SJ, RH, and PV carried out the experiments. GT and FL supervised the study and interpreted the data. All authors read and commented on the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.