Reprogramming of Fundamental miRNA and Gene Expression during the Barley-Piriformospora indica Interaction

The interactions between plants and microorganisms, which are widely present in the microbial-dominated rhizosphere, have been studied. This association is highly beneficial to the organisms involved, as plants benefit soil microorganisms by providing them with metabolites, while microorganisms promote plant growth and development by promoting nutrient uptake and/or protecting the plant from biotic and abiotic stresses. Piriformospora indica, an endophytic fungus of Sebacinales, colonizes the roots of a wide range of host plants and establishes various benefits for the plants. In this work, an interaction between barley and the P. indica was established to elucidate microRNA (miRNA)-based regulatory changes in miRNA profiles and gene expression that occurred during the symbiosis. Growth promotion and vigorous root development were confirmed in barley colonized by P. indica. The genome-wide expression profile analysis of miRNAs in barley root showed that 7,798,928, 6,418,039 and 7,136,192 clean reads were obtained from the libraries of mock, 3 dai and 7 dai roots, respectively. Sequencing of the barley genome yielded in 81 novel miRNA and 450 differently expressed genes (DEGs). Additionally, 11, 24, 6 differentially expressed microRNAs (DEMs) in barley were found in the three comparison groups, including 3 dai vs. mock, 7 dai vs. mock and 7 dai vs. 3 dai, respectively. The predicted target genes of these miRNAs are mainly involved in transcription, cell division, auxin signal perception and transduction, photosynthesis and hormone stimulus. Transcriptome analysis of P. indica identified 667 and 594 differentially expressed genes (DEG) at 3 dai and 7 dai. Annotation and GO (Gene Ontology) analysis indicated that the DEGs with the greatest changes were concentrated in oxidoreductase activity, ion transmembrane transporter activity. It implies that reprogramming of fundamental miRNA and gene expression occurs both in barley and P. indica. Analysis of global changes in miRNA profiles of barley colonized with P. indica revealed that several putative endogenous barley miRNAs expressed upon colonization belonging to known micro RNA families involved in growth and developmental regulation.


Introduction
The interactions between plants and microorganisms, which are widely present in the microbial-dominated rhizosphere, have been well studied. This association is highly beneficial to the organisms involved, as plants benefit soil microorganisms by providing them with metabolites, while microorganisms promote plant growth and development by promoting nutrient uptake and/or protecting the plant from biotic and abiotic stresses [1,2]. The establishment and maintenance of mutualism requires genetic and epigenetic reprogramming and metabolomic regulation through the exchange of effector molecules between beneficial microorganisms and plants [3,4]. Beneficial microorganisms have a major role in crop production because of their impact on plant health and yield. Piriformospora indica (P. indica) is an endophytic fungus, belonging to the order Sebacinales that colonizes the roots of both monocotyledons and dicotyledons plants [5]. P. indica serves as an excellent model for beneficial microbes as it can form a mutually beneficial symbiosis with a series of crops such as Chinese cabbage, rice, wheat, cucumber, onion and banana, which can effectively promote their growth, nutrient absorption, accumulation of secondary metabolites, and resistance to disease damage [6][7][8][9][10][11][12]. P. indica has great potential in biological control and soil improvement and thus can play a positive role in agricultural production.
MicroRNAs (miRNAs) are a class of endogenous small noncoding RNAs (ncRNAs) which is evolutionary conserved and contains approximately 20-22 nucleotides [13,14]. They participate in the regulating gene expression and multiple physiological and biochemical processes by complementary functions with target gene mRNA. Many studies have shown that plant miRNAs play an important regulatory role in the interaction between plants and soil microbes, including promoting plant growth and development, stress response and hormone transduction [15]. MiRNAs recognize their mRNA target genes based on near-perfect complementarity and direct degradation or translational repression of homologous mRNA targets [16]. In addition, they tend to act as "early" regulators of signal transduction at the level of transcription factors (TFS) in various systems [17,18]. MiRNAs respond rapidly to infection by symbiotic bacteria. In soybean roots, a group of miRNAs which target a wide range of mRNAs were intensively up-or down-regulated by infection with the rhizobium bacterium Brodyrnia japonicum [19,20]. In the process of symbiosis, miRNAs were involved in the regulation of plant nutrient balance [21,22], hormone homeostasis and signal transduction [19], and spatial and temporal development of symbiosis nodules [20,23].
In addition, miRNAs also play a regulatory role in the process of abiotic stress in plants. Li et al. [24] found that after silencing BnmiR169n in rape seed, the drought tolerance of plants increased due to the increase of its target gene BnNFYA8. Under salt stress, the expression of miR169q in maize was inhibited and the expression of its target gene ZmNFYA8 was upregulated. ZmNFYA8 binds to the promoter of the antioxidant enzyme gene zmper1 and activates its expression, alleviating the toxic effect of ROS on plants and improving maize salt tolerance [25]. The research on tomato by Zhang et al. [26] also verified that miR394 was involved in the negative regulation of biological stress. The overexpression of miR394 inhibited the expression of its target gene LCR, and then inhibited JA synthesis-related genes, thus reducing the resistance of tomato to Phytophthora. It is found that MiR164 plays an important role in wheat leaf rust and poplar black spot defense [27,28].
The role of miRNAs as gene expression regulators in Sebacinalean symbiosis has been largely unexplored. Previous studies have shown that P. indica can induce root growth of Oncidium orchid which is closely related to microRNA [29]. Results indicated that the predicted miRNAs target genes are mainly participated in auxin signal perception and transduction, transcription, development, and plant defense. Several novel unique miRNAs were detected, for which a function could not yet be identified. Another research revealed fundamental sRNA and gene expression reprogramming at the onset of symbiosis between P. indica and the model grass species Brachypodium distachyon [30]. Their data suggests that a Sebacinalean symbiosis involves reciprocal sRNA targeting of genes during the interaction.
Based on comprehensive high-throughput sequencing and transcriptome analysis, we evaluated an interaction between barley and the beneficial fungal root endophyte P. indica to elucidate miRNA regulatory changes in gene expression and miRNA-mRNA interaction profiles. Additionally, we discuss the biological functions and potential regulatory mechanisms of miRNAs in barley growth and possible miRNA-based regulation that might be crucial for the establishment of the barely-P. indica symbiosis. This study will contribute to our understanding of the RNA-based growth promotion mechanism mediated by P. indica colonization.

Novel miRNA Prediction
The characteristics of hairpin structure of miRNA precursor could be used to predict novel miRNA. We used miREvo [33] and mirdeep2 [34] to predict novel miRNAs and the minimum free energy of the small RNA tags unannotated in the former steps. At the same time, custom scripts were used to obtain the identified miRNA counts as well as base bias on the first position with certain length and on each position of all identified miRNAs, respectively.

Co-Expression Analysis of mRNA-miRNA
The functional annotation of identified miRNAs was performed using co-expression analysis [35]. Pearson's correlation coefficients between mRNAs and miRNAs were calculated based on the mRNAs FPKM values, and the putative target mRNA should have a value >0.99 or <−0.99. In addition, the TargetFinder [36] was used to predict the target mRNA of the miRNA. The mRNA-miRNA network was constructed using Cytoscape [37] software (Version 3.0.2) based on the correlations between mRNAs and miRNAs.

GO and KEGG Enrichment Analysis of Differentially Expressed Genes
Gene Ontology (GO) enrichment analysis of differentially expressed genes was implemented by the cluster Profiler R package, in which gene length bias was corrected. GO terms with corrected p-value less than 0.05 were considered significantly enriched by differential expressed genes. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism, and the ecosystem, from molecular-level information, especially large-scale molecular data sets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/ (accessed on 12 November 2022)). We used cluster Profiler R package and KOBAS software to test the statistical enrichment of differential expression genes. The term with a corrected p-value < 0.05 is considered to be significantly enriched in differentially expressed genes.

Quantitative Real Time Polymerase Chain Reaction (qRT-PCR) and Stem-Loop PCR for Validation of Sequencing Results
The RNA collected in 2.2 and M-MLV Reverse Transcriptase were used to make cDNA. After treatment with Dnase I (Sigma, Germany), the cDNA was used as a template for qRT-PCR to quantify selected miRNAs and mRNAs using the miRNA-specific primers and target mRNA specific primers [38]. The expression level of respective gene was determined by quantitative RT-PCR. Quantitative RT-PCR was measured by SYBR Green influorescence method as described previously. In brief, qPCR experiments were conducted on a Light Cycler96 Fast real-time PCR system (Roche). The reaction solution contains 2 × Ultra SYBR Mixture 10 µL, 100 ng cDNA template, 10 µM forward and reverse primers. HvUBIQUTIN was used as the control, and all experiments were conducted with at least three technical replications. The amplification program was applied as the following steps: the first initial activation step was performed at 95 • C, 5 min, then followed by 30 cycles (95 • C for 20 s, 56 • C for 35 s, 72 • C for 35 s, and 65 • C for 20 s). At the end of each cycle, melting curves were determined respectively to guarantee the amplification of the single-PCR product.
For the identification of miRNAs in barley, stem-loop RT-PCR was referenced [39]. Briefly, cDNA was synthesized from total RNA extracted from P. indica inoculated barley roots. Hairpin primer was designed and performed according to literature [40]. For each stem-loop reaction, the detail protocol was performed according to the manufacturer's instructions (Thermo Scientific, Waltham, MA, USA).
For primer annealing, the reaction was incubated at 16 • C for 30 min and then extended at 42 • C for 30 min. The universal stem-loop primer and specific miRNA primer (Table S5) were used in Endpoint PCR under the same conditions as described in target RNA amplification. PCR products were purified and cloned into the pGEM ® -T Vector Systems (Promega, Madison, WI, USA) following the manufacturer's instructions. Five colonies of each cloned miRNA were subjected to sequencing using an M13 forward primer (Ding guo, China).

P. indica Promote Root Growth and Plant Development
To investigate whether barley seed can develop a beneficial interaction with P. indica by applying seed soaking (SS) treatment with P. indica chlamydospore suspensions, seedling three days after germination in soil were subjected to colonization identification. The hyphae of P. indica were widely distributed over the root surface three days after germination ( Figure 1A) which indicating that the establishment of the beneficial symbiosis is successful. The biomass enhancement was observed in the SS treatment compared to that the non-inoculated barley, the stems and roots developed better than barley without P. indica co-cultivation at 3 and 7 days after inoculation (dai) ( Figure 1B,C). Furthermore, the branching of roots was also obvious ( Figure 1D). Comparison of shoot length in colonized vs. non-colonized plants grown in soil showed that P. indica increased the shoot length by 18.7% and 22.1% at 3 dai and 7 dai, respectively ( Figure 1E), and total grain weight/plant increased by 36.9% and 44.6% at 3 dai and 7 dai, respectively ( Figure 1F). We found that this method of soaking seeds was very effective in increasing barley biomass. Concordantly, shoot and weight analyses of barley seedlings revealed a significant increase in biomass ( Figure 1D) upon P. indica colonization.
the non-inoculated barley, the stems and roots developed better than barley wit indica co-cultivation at 3 and 7 days after inoculation (dai) ( Figure 1B,C). Furthermo branching of roots was also obvious ( Figure 1D). Comparison of shoot length in col vs. non-colonized plants grown in soil showed that P. indica increased the shoot len 18.7% and 22.1% at 3 dai and 7 dai, respectively ( Figure 1E), and total grain weigh increased by 36.9% and 44.6% at 3 dai and 7 dai, respectively ( Figure 1F). We fou this method of soaking seeds was very effective in increasing barley biomass. Co antly, shoot and weight analyses of barley seedlings revealed a significant increas omass ( Figure 1D) upon P. indica colonization.

Establishment of the Barley-P. indica Interaction Is Associated with Extensive Transcriptional Reprogramming
In order to evaluate how the mutualistic interaction affects, the miRNA profiles in the colonized barley root at 3 and 7 dai of P. indica and barley controls were subjected to miRNA sequencing. These time points were selected because the fungus-related growthpromoting effect was visible. After removal of the low-quality contaminant and adapter reads, 7,798,928, 6,418,039, and 7,136,192 clean read sequences were obtained from the libraries of mock, 3 dai, and 7 dai roots (Tables S1-S3), respectively. Sequencing data indicated that libraries are of high quality and can be used for further miRNA studies (Table S4).
Differentially expressed microRNAs in barley were found in the three comparison groups, including 3 dai vs. mock, 7 dai vs. mock, and 7 dai vs. 3 dai (Figure 2). There were 8 microRNAs upregulated and 3 micro RNAs down regulated in 3 dai vs. mock; there were 11 microRNAs upregulated and 13 micro RNAs down regulated in 7 dai vs. mock, and there were 3 microRNAs upregulated and 4 microRNAs down regulated in 7 dai vs. 3 dai. The clustering heat map showed that P. indica colonization affected the expression pattern of microRNAs ( Figure S1); alternatively, some miRNAs like novel-1, miR-444b, miR-1120 and miR397a were down regulated responding to the P. indica colonization; whereas miR171-3P and miR 6183 were upregulated at 3 dai and down-regulated at 7 dai; the other microRNAs including miR1681, novel-22, novel-3, miR168-3P, miR6198, miR6177, and novel 13 were up-regulated after P. indica colonization. To investigate interaction of barley by P. indica, differentially expressed genes (DEGs) in colonized P. indica in comparison to axenic mycelium samples were analyzed. Among the 667 DEGs of predicted unique P. indica genes, 499 were confirmed as up-regulated and 168 were down-regulated in 7 dai vs. control, and of 594 predicted unique P. indica genes, 517 were confirmed as upregulated and 168 were down-regulated at day 3 after colonization ( Figure 3). Gene ontology and annotation analysis indicated that the DEGs with the greatest changes were mainly concentrated in oxidoreductase activity, ion transmembrane transporter activity, phosphate transporter and tRNA THr modification in 3 dai vs. control, and other DEGs were concentrated in cellular respiration, ion transport, endonuclease activity, oxidoreductase activity, and phosphate transporter in 7 dai vs. control ( Figure 4). Top 20 P. indica DEGs during colonization at 3 dai and 7 dai vs. mock and 7 dai vs. 3 dai are listed in Tables 1-3. transporter activity, phosphate transporter and tRNA THr modification in 3 dai vs. control, and other DEGs were concentrated in cellular respiration, ion transport, endonuclease activity, oxidoreductase activity, and phosphate transporter in 7 dai vs. control ( Figure 4). Top 20 P. indica DEGs during colonization at 3 dai and 7 dai vs. mock and 7 dai vs. 3 dai are listed in Tables 1-3.

Differentially Expressed miRNAs
Analysis of unique plant miRNAs in 7 dai-mock vs. 3 dai-mock revealed that 12 of the putative endogenous miRNAs were exclusively present in 7 dai-mock, 6 miRNAs were exclusively present in 3 dai-mock, and 12 miRNAs were present in both comparison group ( Figure 5). For the reads from the putative 3 dai-mock vs. 7 dai-3 dai, 13 miRNAs were exclusively present in 3 dai-mock, 4 miRNAs were exclusive present to 7 dai-3 dai, and 5 miRNAs were found in both comparison groups. Comparison between the unique miRNAs in 7 dai-mock vs. 7 dai-3 dai indicated that 19 of the putative endogenous miRNAs were exclusively present in 7 dai vs. mock, 4 miRNAs were found in 7 dai vs. 3 dai, and 5 miRNAs were present in both comparison group. Similarly, from the putative 3 dai-mock vs. 7 dai-3 dai, three miRNAs were exclusively present in 3 dai-mock, nine miRNAs were exclusively present in 7 dai-mock, one miRNA was exclusively present in 7 dai-3 dai, and two miRNAs were found in all three comparison groups.

Expression Patterns of miRNAs and Their Putative Targets in P. indica-Colonized Roots
To elucidate the regulatory function of miRNAs on their putative targets, real-time qPCR was performed to confirm the expression level by using specific primers for these miRNAs and their target genes (Table S5). The miRNAs novel_1, novel_22, novel_40, miR444b and their target genes involved in the growth-regulating factor, promoter-bindinglike protein and transcription factor were selected. QPCR results indicated that the selected miRNAs novel_1 was down-regulated at 3 dai compared with mock, and its target genes HORVU7Hr1G012380. 5 and HORVU5Hr1G055920. 2 were both down-regulated, too ( Figure 6A). In addition, the selected miRNAs miR444b was down-regulated at 3 dai compared to mock, and its target gene HORVU3Hr1G076030. 14 was up-regulated, but the other target gene HORVU1Hr1G006020. 8 was down-regulated ( Figure 6B). Moreover, the selected miRNAs novel_22 was up-regulated at 3 dai compared with mock, and its target gene HORVU4Hr1G082910. 18 was upregulated, but the target gene HORVU2Hr1G085210. 5 was downregulated ( Figure 6C). And the selected miRNAs novel_40 was down-regulated at 3 dai compared to mock, and its target gene HORVU5Hr1G054420. 4 was up-regulated, another target gene HORVU7Hr1G088630. 3 was down-regulated ( Figure 6D). The result showed that the expression of these miRNAs was consistent with the results of RNAseq data.
Functions of target mRNA corresponding to up-regulated and down-regulated miRNA in 3 dai vs. mock are listed in Tables 4 and 5. As shown, one miRNA corresponded to multiple target genes, and the regulatory function of miRNA on target genes may be either up-regulated or down-regulated. Data showed that four microRNAs corresponding to target RNAs had specific functional descriptions. And the functions of these target mRNA described as: Serine/threonine-protein kinase STY13, putative transcription factor RL9, leucine rich repeat family expressed, 60S ribosomal protein L35a-1, squamosa promoterbinding-like protein, squamosa promoter-binding-like protein 4/16, scarecrow-like protein 6/. And the functions of down-regulated miRNA related target mRNA were MADS-box transcription factor 57, 6-phosphogluconate dehydrogenase, general negative regulator of transcription subunit 3 isoform X4, and growth-regulating factor 1, 2, 4, 5-like isoform.

Expression Patterns of miRNAs and Their Putative Targets in P. indica-Colonized Roo
To elucidate the regulatory function of miRNAs on their putative targets, real-t qPCR was performed to confirm the expression level by using specific primers for th miRNAs and their target genes (Table S5). The miRNAs novel_1, novel_22, novel_ miR444b and their target genes involved in the growth-regulating factor, promoter-bi ing-like protein and transcription factor were selected. QPCR results indicated that selected miRNAs novel_1 was down-regulated at 3 dai compared with mock, and its get genes HORVU7Hr1G012380. 5 and HORVU5Hr1G055920. 2 were both down-re lated, too ( Figure 6A). In addition, the selected miRNAs miR444b was down-regulate 3 dai compared to mock, and its target gene HORVU3Hr1G076030. 14 was up-regula but the other target gene HORVU1Hr1G006020. 8 was down-regulated ( Figure 6B). Mo over, the selected miRNAs novel_22 was up-regulated at 3 dai compared with mock, its target gene HORVU4Hr1G082910. 18 was upregulated, but the target g HORVU2Hr1G085210. 5 was downregulated ( Figure 6C). And the selected miRN novel_40 was down-regulated at 3 dai compared to mock, and its target g HORVU5Hr1G054420. 4 was up-regulated, another target gene HORVU7Hr1G08863 was down-regulated ( Figure 6D). The result showed that the expression of these miRN was consistent with the results of RNA-seq data.   Functions of target mRNA corresponding to up-regulated and down-regulated miRNA in 3 dai vs. mock are listed in Tables 4 and 5. As shown, one miRNA corresponded to multiple target genes, and the regulatory function of miRNA on target genes may be either up-regulated or down-regulated. Data showed that four microRNAs corresponding to target RNAs had specific functional descriptions. And the functions of these target mRNA described as: Serine/threonine-protein kinase STY13, putative transcription factor RL9, leucine rich repeat family expressed, 60S ribosomal protein L35a-1, squamosa promoter-binding-like protein, squamosa promoter-binding-like protein 4/16, scarecrow-like protein 6/. And the functions of down-regulated miRNA related target mRNA were MADS-box transcription factor 57, 6-phosphogluconate dehydrogenase, general negative regulator of transcription subunit 3 isoform X4, and growth-regulating factor 1, 2, 4, 5-like isoform.

Prediction of miRNA Target Genes
To understand the regulatory function of miRNAs during symbiosis, those miRNAs which were abundantly detected and significantly up-/down-regulated by P. indica were selected for further investigation. Eighty-one conserved miRNAs belonging to 15 families were selected for target gene prediction. Target-Finder software was used to predict miRNA target genes, 450 best fit target candidates were obtained (Table S6).
Subsequently, annotation and GO analysis were conducted by Goseq. A total of 450 targets were annotated and distributed in 38 categories. Clustering of microRNA target genes varied among different comparison groups. MiRNAs target genes were more concentrated in biological process signaling pathways than in cell components in 3 dai vs. mock ( Figure 7A): In the subcategory of biological process, most of the target genes of miRNA are concentrated in the cellular biosynthetic process, macromolecule biosynthetic process and nucleic acid metabolic process; In the subcategory of molecular function, transferase activity and copper ion binding are the main enrichment pathways of miRNA target genes. In the 7 dai vs. mock group, miRNA target genes were involved in biological processes, cell components, and molecular functional signaling pathways ( Figure 7B): most of the miRNA target genes are concentrated in the cellular biosynthetic process, macromolecule biosynthetic process, nucleic acid metabolic process and organic substance metabolic process in the subcategory of biological process, In the subcategory of cell components, most of the miRNA target genes are concentrated in intracellular organelle, membrane bounded organelle and nucleus; in the subcategory of molecular function, transferase activity, ribonucleotide binding, carbohydrate derivative binding, and anion binding are the main enrichment pathways of microRNA target genes. In the 7 dai vs. 3 dai group, miRNAs target genes were more concentrated in molecular function signaling pathways than in biological process ( Figure 7C). In the subcategory of biological process, the expression of miRNA target genes was mainly enriched in protein phosphorylation and phosphorous compound metabolism. Whereas, miRNA target genes, in the subcategory of molecular function were mainly concentrated in binding process of small molecules, ATP, ADP, copper ion, carbohydrate derivative, and purine nucleoside.

Identification of miRNAs Related to Transcription Factor and Other Key Pathway Regulation in Barley Roots Colonization by P. indica
To further study the function of miRNAs in the growth of barley roots in response to P. indica colonization, the GO terms of different expression miRNAs targets were annotated based on http://geneontology.org/ (accessed on 12 November 2022) and http://www. uniprot.org/uniprot/ (accessed on 12 November 2022). The results showed that seven miRNAs might be involved in the regulation of gene transcription, because their target genes have transcription factor activity ( Table 6). In the 3 dai vs. mock group, miR6189 was up-regulated and miR 6214 and miR444b were down-regulated; in the 7 dai vs. mock group, miR6190 was up-regulated and miR6214, miR397a were down-regulated; in the 7 dai vs. 3 dai group, miR6180 was down-regulated and miR6214 was up-regulated; all their target genes have transcription factor activity.
Besides the transcription factor activity regulation, we found that miRNAs including miR6214, miR6180, miR6189, miR444b, miR6190, and miR397a might also be involved in other key signal pathways because their target genes participated in these pathways including cell division, auxin stimulus, photosynthesis, hormone stimulus, and chlorophyllide oxygenase activity, which contribute to the plant growth promotion. It was observed that there were four miRNAs response to auxin stimulus in the 3 dai vs. mock and 7 dai vs. mock comparison groups. Those miRNAs were all up-regulated compared with mock while they were diversely (positively or negatively) correlated to their target genes (Table 7). The chart shows that miRNA_novel 1 was downregulated and its target genes HORVU7Hr1G012380 and HORVU5Hr1G055920 were down-regulated. (B) The chart shows that miRNA_444b was downregulated and its target genes HORVU3Hr1G076030 was upregulated but HORVU5Hr1G006020.8 were down-regulated. (C) The chart shows that miRNA_novel 22 was upregulated and its target Figure 7. The relative expression of miRNAs and their target genes involved in the growth-regulating factor, promoter-binding-like protein, and transcription factor. (A) The chart shows that miRNA_novel 1 was downregulated and its target genes HORVU7Hr1G012380 and HORVU5Hr1G055920 were down-regulated. (B) The chart shows that miRNA_444b was downregulated and its target genes HORVU3Hr1G076030 was upregulated but HORVU5Hr1G006020.8 were down-regulated. (C) The chart shows that miRNA_novel 22 was upregulated and its target genes HORVU4Hr1G082910 was upregulated, too. The other target gene HORVU2Hr1G085210 was downregulated.

Discussion
To cope with rapidly changing environments, plants employ a large number of mechanisms that provide phenotypic plasticity and allow fine-tuning of stress response actions. Advances in molecular biology have made great strides in identifying genomic regions and underlying mechanisms that influence transcriptional and post-transcriptional biotic and abiotic stress pathways regulation. In plants, miRNAs evolve and contribute to the complexity of these molecules through a series of pathways, and play an important role in the regulation of stress response activity. It has been shown that one way in which plants respond to environmental stress is through the activity of miRNAs. MiRNAs, as important regulatory molecules of plant biotic and abiotic stress response, are the driver molecules of RNA interference (RNAi), and ensure the up-regulation and down-regulation of target genes, and participate in important biological processes [41]. MiRNAs are small, 20-22 nt noncoding RNAs that regulate gene expression by post-transcriptional gene silencing in most eukaryotes [42]. RNA interference (RNAi) regulates gene expression by inducing degradation of messenger RNA (mRNA) or inhibiting its translation. MiRNAs play crucial roles in plant development, including the formation of embryo, meristem, leaf, and flower [43] as well as the responses to biotic and abiotic stresses [44]. P. indica is well known to be able to establish beneficial interactions with many different hosts, including monocotyledons such as barley, wheat, rice, corn, and dicotyledons such as Arabidopsis and tobacco [8], even on Brassicaceae family that cannot be colonized by mycorrhizal fungi [11]. Colonization of the roots by P. indica results in enhanced biomass production as well as increased resistance against biotic and abiotic stresses. We established and studied the interaction between barley-a major cereal crop-and P. indica-a beneficial endophyte with an exceptionally large host range. The functions and regulatory mechanisms of miRNAs in barley growth regulation during the symbiosis with P. indica were explored in our work. We showed that P. indica colonizes barley, resulting growth promotion in shoot, alterations in root architecture, and improved grain development.

Transcriptional Changes Detected during the Barley-P. indica Interaction
To investigate the interaction of barley with P. indica, we analyzed DEGs in colonized P. indica in comparison to axenic mycelium samples. Gene ontology analysis indicated DEGs with the greatest changes were mainly concentrated in oxidoreductase activity, ion transmembrane transporter activity, phosphate transporter, and tRNA THr modification in 3 dai vs. control, and other DEGs were concentrated in cellular respiration, ion transport, endonuclease activity, oxidoreductase activity and phosphate transporter in 7 dai vs. control. We noticed that the phosphate transporter gene which could promote the uptake of phosphorus by plant roots and provide essential nutrients for plants [45] was upregulated in both 3 dai vs. control and 7 dai vs. control. These results indicate that the interaction between barley and P. indica initiates the expression of proteins related to phosphorus transport. Other up-regulated genes encode enzymes involved in ion transport and ion transmembrane transporter activity. These proteins may modulate metal ions transport from soil into the plants. In fact, there are many transporters in P. indica which are responsible for the transfer of phosphorus [45] and sulfur elements [46] and the transport of metal elements. It is believed that these transporters play an important role in assisting the plant to absorb large and trace elements.
Roots of barley plants also displayed substantial transcriptional reprogramming following P. indica colonization. GO analysis indicated DEGs in barley root enrichment in nucleic acid metabolic and macromolecule biosynthetic, cellular marcomolecule biosynthetic and transferase activity associated processes in 3 dai vs. mock. Barley root DEGs exhibiting the greatest changes in expression between 7 dai and noncolonized plants are related to the primary metabolic process, cellular process, organic substance, transferase activity, and phosphotransferase activity. Previous research indicated that DEGs in Brachypodium distachyon after P. indica colonization involved in catalytic and oxidoreduction associated processes [30]. This also indicates that different plants have different transcriptome responses upon P. indica colonization. In 3 dai vs. mock, of the downregulated barley genes, several encode proteins commonly associated with stress responses, including a peroxidase and a putative protein kinase. Several downregulated DEGs also encode transcription factors, including growth-regulating factor 1/2/4/5/6-like (Table S7). One report also demonstrated that several DEGs encode transcription factors, including MYB-related, GRAS, and bZIP were downregulated in Brachypodium distachyon after P. indica colonization. Of the upregulated barley genes, several encode leucine rich repeat family and squamosa promoter-binding-like protein (Table S7). In 7 dai vs. mock, several upregulated DEGs also encode squamosa promoter-binding-like protein 16, serine/threonine-protein kinase SIS8, and some down regulated DEGs encode MADS-box transcription factor 57, pyruvate kinase, cytosolic, isozyme, and transcription initiation factor TFIID subunit 6 (Table S8). In 7 dai vs. 3 dai, the downregulated target genes were mainly encoded scarecrow-like protein and growth-regulating factor (Table S9).

Barley miRNAs Detected in the Barley-P. indica Interaction
The role of miRNAs as gene expression regulators in P. indica symbiosis has been largely unexplored in barley. The report demonstrated that P. indica promot plant growth associated miRNAs in Oncidium orchid roots [29]. Another study indicated that sRNAs reprogrammed after P. indica colonization in Brachypodium distachyon [30]. A high-throughput sequencing and comparative expression analysis were conducted. In total, 7,798,928, 6,418,039 and 7,136,192 clean reads were obtained from the libraries in control and P. indicacolonized root libraries. Differentially expressed microRNAs were found in the three comparison groups, including 3 dai vs. mock, 7 dai vs. mock and 7 dai vs. 3 dai (Figure 2). Analysis of putative endogenous barley miRNAs expressed during P. indica colonization identified 42 miRNAs. Some of them have unknown targets, whereas some of them have more than one target. For example, the miRNA hvu-miR1120 targets mRNAs (HORVU1Hr1G080480) encoding 6-phosphogluconate dehydrogenase and other targets mRNAs (HORVU4Hr1G002170) encoding general negative regulators of transcription.
Four key miRNAs, including hvu-miR6181, novel_22, novel_44, and hvu-miR6198, were up-regulated both in 3 dai and 7 dai vs. mock in our work. Their targets gene functions were mainly involved in vacuolar protein-sorting-associated protein and squamosa promoterbinding-like protein. Other important miRNAs, such as hvu-miR1120, hvu-miR444b, novel_1, and hvu-miR397a, were down-regulated both in 3 dai and 7 dai vs. mock. Among them, hvu-miR6189, hvu-miR6214, hvu-miR444b, hvu-miR6190, hvu-miR397a, their targets mRNAs actually possess transcription factor activity (Table 4). In Arabidopsis, repression of transcription factors by the miR165/166 family modulates root growth, maintenance of the shoot apicalmeristem, and the development of leaf polarity; miR156-mediated downregulation of SPLs modulates developmental timing, lateral root development, branching, and leaf morphology [30]. This suggests that miRNAs, which regulate transcription factors, play an important role in plant growth and development. MiRNA including mir156, mir166, and mir169, which target mRNAs for the transcription factor genes SPL, PHV/PHB, and NTF, were also abundantly detected in Oncidium orchid roots after P. indica colonization [29]. This further indicates that the colonization of P. indica could reprogram miRNAs which participate in transcription factors regulation.
These miRNAs identified in barley include hvu-miR6180, hvu-miR6189, hvu-miR6214, hvu-miR444b, and hvu-miR6190, predicted to target genes involved in hormone activity, cell division, and photosynthesis pathway (Table 5). Because these miRNAs predict a variety of targets that are associated with plant growth and development, this group of miRNAs may play an important role in reprogramming plant cells during P. indica symbiosis establishment.

MiRNA and Target mRNA
MiRNAs regulate gene expression by mediating target gene silencing at transcriptional (TGS) and post-transcriptional (PTGS) levels, including DNA methylation, histone modification, translational repression, and RNA silencing [47,48]. They play important roles in plant development, differentiation, and response to biotic and abiotic stresses [48][49][50][51][52][53][54][55]. In the analysis of the interaction between miRNA and target genes, it was found that miRNA can be positively correlated with target genes or negatively correlated with target genes at 3 dai and 7 dai of P. indica. The predicted target genes of these miRNAs are mainly involved in transcription, cell division, auxin signal perception and transduction, photosynthesis, and hormone stimulus (Tables 4 and 5). In fact, miRNAs can bind not only to mRNAs but also to long non-coding RNAs (lncRNA). LncRNAs acting as potential competing endogenous RNA-harboring microRNA response elements (MREs), thereby competing with mRNAs for shared miRNA and thus regulate miRNA-mediated gene silencing [44]. In our work, the lncRNA-mRNA-miRNA network was constructed, too ( Figures S2 and S3). For instance, several key lncRNAs were found to have correlations with the miRNA including novel_1, novel_19, novel_22, novel_44hvu-miR6181, hvu-miR6198, and hvu-miR444b in 7 dai vs. mock; several key lncRNAs were also found to have correlations with the miRNAs, including novel_1, novel_2, novel_22, hvu-miR6181, hvu-miR6198, and hvu-miR444b in 3 dai vs. mock. These data suggest that these key miRNAs excavated play an important role in the regulation of plant growth in response to P. indica colonization. Specific regulatory functions of miRNAs on target genes and lncRNAs is the main research content in the future work.

Conclusions
We reported the miRNA profiling of barley after colonization by P. indica. The miR-NAs and their target genes illustrated that the physiological metabolism of barley is reprogrammed in response to the symbiotic interaction. Genes participating in transcription, cell division, auxin signal perception and transduction, photosynthesis, and hormone stimulus are major targets of the P. indica-induced miRNAs in barley roots. Therefore, we propose that P. indica alters the miRNA pattern to establish an intricate network for growth promotion and developmental reprogramming and enhances resistance in barley roots. Several novel unique miRNAs were detected, for which a function could not yet be identified. Further investigations on the molecular mechanism of miRNAs in symbiotic interactions are of huge significance.
Author Contributions: L.L.: conceptualization, methodology, investigation, formal analysis, visualization, writing-original draft preparation, reviewing, and editing. N.G.: investigation, formal analysis, visualization. Y.Z. and Z.Y.: software, data curation, formal analysis. S.L. and A.L.: investigation, formal analysis, writing-reviewing and editing, Z.W. reviewing and funding the project. All authors have read and agreed to the published version of the manuscript. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: All data generated or analyzed during this study are available in this article and its supplementary information files. All miRNA-sequencing data related to the present study have been deposited in the National Centre for Biotechnology Information (NCBI) under the Bioproject accession number PRJNA898289. (https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA898 289 (accessed on 12 November 2022)).