Phialocephala fortinii F5 was stored at −4°C on potato dextrose agar (PDA, containing 200 g/l⁻¹ potato extract, 20 g/l⁻¹ glucose, and 20 g/l⁻¹ agar, pH 7.0). The sequence of this strain’s rDNA internal transcribed space was submitted to the NCBI database under accession number KY910213.1. After 2 weeks of culture on PDA, the mycelia of F5 were collected using a sterile 5-mm diameter cork borer.

Raspberry seeds were sterilized in 75% alcohol solution for 30 s, washed three times with sterilized water, transferred to 2% sodium hypochlorite for 5 min, and then washed 5 times with sterilized water. The micropropagated seedlings of raspberry were transferred to a WPM medium. The seedlings were grown in the culture room at 28°C. After 3 months, aseptic raspberry seedlings were planted in a soil matrix (80% turfy soil and 20% vermiculite) and grown in a controlled growth room (25°C, 16/8 h, day/night, 1,600 lux). Twenty seedlings were inoculated with an F5 mycelium pellet obtained by liquid culture, and the remaining 20 were grown without F5 treatment as a control. The samples were labeled CK (control check, without inoculation), BI (presymbiotic phase, 3 weeks after inoculation), and UI (symbiotic phase, 5 weeks after inoculation). Sample roots were collected and separated into two parts. One part was used to check the F5 colonization rate, and the other part was flash frozen in liquid nitrogen and stored at −80°C for RNA extraction. To determine growth parameters, such as plant height, root length, fresh weight, and root activity (Peng et al., 2022), three biological replicates were performed to calculate the mean expression and standard deviation (SD).

Raspberry roots were fixed with formalin–acetic acid–alcohol (FAA) for 12 h. The roots were then treated with 10% KOH and heated for 1 h, followed by treatment with 1% HCl for 30 min at room temperature. Finally, the roots were stained with Trypan Blue (Liu et al., 2016). Dyed roots were detected after transparency with lactic acid and glycerin.

Total RNA was extracted from the collected root samples with TRIzol^® Reagent (Invitrogen, United States) according to the manufacturer’s instructions. Total RNA concentration and purity were measured by a NanoDrop 2000 (Thermo Scientific), and total RNA integration and quality were determined by agarose gel electrophoresis. RNA molecules that were less than 50 nt in length were enriched and ligated with proprietary adapters. The RNA samples ligated with adapters were reverse transcribed and amplified by PCR to produce sequencing libraries. The three sRNA libraries from raspberry roots were sequenced on an Illumina Hiseq 2,500 (Santiago, CA, United States) platform at Genedenovo Biotechnology Co., Ltd. (Guangzhou, China). The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive (Genomics, Proteomics & Bioinformatics 2021) in National Genomics Data Center (Nucleic Acids Res 2022), China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (GSA: CRA008006) that are publicly accessible at https://ngdc.cncb.ac.cn/gsa.

High-quality reads were generated by removing the adaptor reads, low-quality sequences, and fragments of less than 18 nt or more than 30 nt from the raw reads. To classify the clean reads, they were aligned with the small RNAs and annotated to Rfam databases (Rfam 11.0, rfam. janelia.org) to identify and filter out the ribosomal RNA (rRNA), transport RNA (tRNA), nuclear small RNA (snRNA), nucleolar small RNA (snoRNA), and repeat sequences. After the removal of known miRNAs, the remaining sequences were used to identify novel miRNAs according to their hairpin structures predicted by Mireap_V0.2 software. Novel miRNAs were predicted from non-annotated sRNA sequences using miRDeep2 software (Friedländer et al., 2012). The basic criteria were used to screen for new candidate miRNAs (Zhang et al., 2015).

To identify F5 symbiosis-responsive miRNAs, raw miRNA counts from each sample were used to normalize miRNA expression. The fold changes in miRNA levels in treated samples were calculated relative to those in the control. The miRNAs with a fold-change in expression log² (treatment/control) ≥ 2 and an adjusted value of p < 0.05 were considered differentially expressed and referred to as the F5 fungus-responsive miRNAs. Tool DEseq (version 1) was used to identify the differentially expressed miRNAs.

Based on the sequences of known and novel miRNAs, the putative target genes were predicted by the software psRNATarget algorithm with the following parameters: (1) maximum expectation is less than 5.0; (2) length for complementary scoring (hspsize) is shorter than 19; (3) range of central mismatch leading to translational inhibition is 10–11 nt. The mature sequences of all miRNAs were submitted to predict target genes; and (4) the minimum free energy (MFE) of the miRNA/target duplex should be greater than or equal to 74% of its miRNA bound to a perfect complement.

To elucidate the regulatory function and pathways involved in symbiosis, the predicted target genes were further annotated and enriched against GO using Blast2GO (version 2.8, accessed online)¹ and against the KEGG pathway database using KAAS (KEGG Automatic Annotation Server, accessed online)², respectively.

Stem-loop primers were designed for qRT-PCR, as described by Chen et al. (2005). The primers are shown in Supplementary Table S1. To verify the relative expression of the miRNAs in the CK, BI, and UI groups, a reverse transcriptase reaction was performed using a HiFiScript cDNA Synthesis Kit (CW2569M, CWBIO, Beijing, China). The 20 ml mixture contained 0.5 mM each of dNTPs, stem-loop primers, 5 mg RNA template, 10 mM DTT, 200 U HiFiScript, and RNase-free water. Then, reverse transcription products were mixed with TB Green Premix Ex Taq (Tli RNaseH Plus, TaKaRa, Dalian, China) in a 96-well plate to start the real-time PCR reaction with a Roche LightCycler 480 II system using the following conditions: an initial denaturation step of 30 s at 95°C, 40 cycles of denaturation for 5 s at 95°C, 30 s at 60°C for annealing, and 30 s extension at 72°C. The relative expression of miRNAs was calculated based on the abundance of the reference gene U6 snRNA. The 2^–△△CT method was adopted to assess the relative miRNA expression from the qRT-PCR experiments (Livak and Schmittgen, 2001). Three biological replicates with three technical replicates were performed for each sample.

To evaluate the plant growth of DSE F5 on the growth of raspberry, we used F5 to inoculate R. idaeus seedlings. After 3 months, F5 inoculation had a dramatically positive effect on the fresh weight, plant height, root length, and root activity of raspberry. The F5-colonized raspberry seedlings exhibited more vigorous growth than the controls (Figures 1A,B). In particular, the root length of the F5-treated raspberry seedlings had significant advantages over that of the control. After 3 months of growth, F5-inoculated seedlings had a 43% greater fresh weight and 56% greater plant height than the control seedlings. Similarly, the root length and root activity of the F5-treated raspberry seedlings were 1.5 times greater than those of the controls (Figure 1G).

Microscopic observation confirmed that the fungal hyphae of F5 colonized the raspberry roots. Darkly pigmented hyphae were observed within the epidermal and cortical cells of F5-inoculated root hairs using trypan blue staining (Figures 1D–F). The features were identical to those previously described for DSEs. In contrast, the root cells of uninoculated seedlings were cleaner and showed no signs of infestation (Figure 1C).

This study aimed to identify miRNAs responsive to symbiosis between raspberry and DSEs. A total of nine sRNA libraries comprising three samples (CK: control check; BI: the period before infection; UI: the period under infection) were generated using the Illumina sequencing platform. Three biological replicates were performed per treatment. Sequencing results showed that there were 15,157,505, 13,487,132, and 10,542,855 raw reads for CK, BI, and UI, respectively, while the clean reads comprised over 73.20% of all reads in the appropriate sizes of 18–35 nt without low-quality or substandard reads (Table 1). By analyzing the lengths of clean read sequences in samples from different periods, the sRNA sequence lengths were found to be mainly distributed between 18 and 30 nt, and mainly concentrated in 21–24 nt, with the largest number of sRNAs at 21 nt. The sequence length decreased sharply in abundance after 25 nt (Figure 2). The distribution of sRNAs in this study was similar to that in other dicotyledons (Rajagopalan et al., 2006; Moxon et al., 2008). Different lengths of sRNAs have different biological functions. This difference in abundance suggests that the expression of miRNAs in samples at different stages can be significantly different.

After the sequencing data were filtered and quality-controlled, the clean reads were compared using the GenBank and Rfam databases to count those annotated to rRNA, tRNA, snRNA, and snoRNA. The proportion of rRNA was much higher than the proportion of other types, followed by tRNA and snRNA, and the lowest proportion was of snoRNA (Supplementary Table S2). The number of all annotated known and novel miRNAs distributed among the three comparisons of raspberry is shown in Figure 3. After removing small non-coding RNAs, such as rRNA, tRNA, snRNA, and snoRNA, the remaining valid sequences were compared with the mature sequences of plant miRNAs in the miRBase database.

A total of 361 known miRNAs were identified from three libraries. We describe the top 20 largest miRNA family reads among all miRNAs in Figure 4. Most miRNAs with high expression (> 100 in each group) were 21 nt in length. The hairpin structure of 50 miRNAs was successfully predicted using mireap_v0.2; the length of the miRNA hairpin ranged from 79 to 295 nt. The minimal folding free energy index (MFEI) of hairpin structures ranged from 0.51 to 1.18, and over 94% of the MFEIs of miRNAs were greater than 0.80. Detailed miRNA information is shown in Supplementary Table S3. These miRNA families were characterized by the presence of numerous miRNAs, as expected. Some miRNAs were only detected in a particular group, showing that certain miRNAs were group-specific (Supplementary Table S4). Most of these miRNA families were highly conserved in plants, such as miR159, miR166, and miR396. Among these, miR396 was the most conserved and existed in 46 plant species. This was followed by miR166 and miR159, which were present in 45 and 10 plant species, respectively. In general, the expression of known miRNAs in raspberry was not abundant; 85 known miRNAs were present in three comparisons at high abundance (> 100), and 276 known miRNAs were at low abundance (about 100). Most of these miRNAs have unknown functions in F5–raspberry symbiosis.

After removal of the known miRNAs and other noncoding RNAs, the remaining clean reads mapped to the raspberry transcriptome were used to predict novel miRNAs based on the secondary structure. Using Mireap_V0.2 software, 95 types of novel miRNAs were identified, and a list of novel miRNA sequences, hairpin MFEIs, precursor sequences, and structures are shown in Supplementary Table S5.

The expression levels of the novel miRNAs among the CK, BI, and UI libraries were compared (Figure 5). In CK, BI, and UI, 81, 85, and 59 novel miRNAs, respectively, were identified, showing that the number of novel miRNAs was much smaller than the number of known miRNAs. Based on Figure 5, we concluded that most of the novel miRNAs were equally expressed among the three libraries, while several of the differentially expressed novel miRNAs deserve further attention to discern their functions. For example, novel miRNAs, such as novel-m0090-5p, novel-m0091-3p, novel-m0094-3p, and novel-m0095-3p, were highly expressed only in the BI group; novel-m0072-5p, novel-m0073-5p, novel-m0076-5p, and novel-m0078-5p were highly expressed in the CK and BI groups but not in the UI group. Detailed information on the predicted novel miRNAs is listed in Supplementary Table S6.

To identify the differentially expressed miRNAs related to F5 in raspberry, we detected miRNA expression levels in nine sRNA libraries (including three biological duplicates). The relative expression of miRNAs from different libraries was calculated using log₂, with a fold change ≥ 2 and p-value < 0.05 as screening conditions. The results showed that 13, 22, and 21 miRNAs were differentially expressed in CK vs. BI (6 upregulated and 7 downregulated), BI vs. UI (15 upregulated and 7 downregulated), and CK vs. UI (16 upregulated and 5 downregulated), respectively (Supplementary Tables S7–S9).

Differential expression analysis revealed potential key regulatory miRNAs during F5 infestation in raspberry. A total of 25 known miRNAs were differentially expressed in the three comparisons. A total of 6 (4 upregulated and 2 downregulated) differentially expressed miRNAs were found in the CK vs. BI comparison group. In total, 15 and 16 known miRNAs were upregulated in the BI vs. UI and CK vs. UI comparison groups, respectively, of which 14 miRNAs were identical. The log2 of miR171-z, miR396, miR11532-z, miR11561-z, and miR1311-x exceeded 19.05, while miR156-y and miR7767-x were downregulated in both groups, indicating that the symbiosis-related miRNAs formed in vivo by F5 remained largely unchanged at the middle and late stages of raspberry infestation, acting together to promote symbiosis formation.

Additionally, we analyzed novel miRNA expression patterns. Fourteen novel miRNAs were differentially expressed, with a total of 7 (2 upregulated and 5 downregulated) differentially expressed miRNAs in CK vs. BI. A total of 4 miRNAs were downregulated in BI vs. UI, and 4 (1 upregulated and 3 downregulated) differentially expressed miRNAs were identified in CK vs. UI. Among them, novel-m0002-3p and novel-m0021-3p were downregulated in both the BI vs. UI and CK vs. UI comparison groups. Overall, the upregulated miRNAs were more abundant than the downregulated miRNAs in plants during the F5 infestation of raspberry. These results suggest that both known and novel miRNAs may play specific and important roles in the response of raspberry to DSE infestation.

miRNAs perform their functions by regulating the expression of their target genes or repressing their translation. Therefore, the prediction and identification of target genes is an important part of understanding the various biological functions of miRNAs. We performed target gene prediction using psRNATarget, and all miRNAs were predicted to have corresponding target genes. Among the identified miRNA-target pairs, target gene prediction was generally not a simple one-to-one relationship but rather a mismatched relationship. A single miRNA may regulate multiple target genes. In this article, 24,112 target genes were predicted for known miRNAs, while 1,367 target genes were obtained for novel miRNAs (Supplementary Tables S10, S11).

There was significant differential miRNA expression between CK, BI and UI (Figure 6). To further investigate the role of responsive symbiotic miRNAs in F5–raspberry symbiosis, we selected miRNAs differentially expressed in F5-infected raspberries for target gene prediction in the middle and late stages. A total of 561 target genes targeted by 37 differentially expressed miRNAs were identified. The annotations of some of the target genes responding to the co-expressed miRNAs are shown in Table 2. A large proportion of these target genes had specific or presumed functions, and these target genes were involved in the regulation of diverse metabolic processes.

To further analyze the possible regulatory pathways involved in miRNA target genes, GO and KEGG functional enrichment analyses were performed for all target genes in this study. GO enrichment analysis showed that many target genes were associated with biological process (GO:0008150), metabolic process (GO:0008152), cellular process (GO:0009987), organic substance metabolic process (GO:0071704), and cellular metabolic process (GO:0044237) in the biological process. For cellular components, the main terms were cell (GO:0005623), cell part (GO:0044464), and membrane (GO:0016020), and important terms for molecular functions were binding (GO:0005488), catalytic activity (GO:0003824), and transferase activity (GO:0016740; Figure 7).

Kyoto Encyclopedia of Genes and Genomes functional enrichment showed that the pathways with the highest enrichment of target genes were Biosynthesis of secondary metabolites (ko01110), Plant–pathogen interaction (ko04626), Glycerophospholipid metabolism (ko00564), Oxidative phosphorylation (ko00190), and Biosynthesis of amino acids (ko01230; Figure 8). Notably, five genes were annotated to the plant–pathogen interaction pathway, namely Unigene0062736 (CPK2), Unigene0014710 (MEKK1), Unigene0047782 (EFR), Unigene0017630 (CML38), and Unigene0020119 (At4g27190). All five genes are involved in plant immune responses, which play an important role in symbiosis formation between plants and rhizosphere microorganisms and are usually specifically repressed from rhizosphere microorganisms, thus promoting the invasion and proliferation of rhizosphere microorganisms (Cao et al., 2017). In this study, five genes were targeted by miR11011-z, miR2109-z, miR8007-y, and miR390-x. All four miRNAs were upregulated in symbiosis; therefore, it was hypothesized that miRNAs may affect F5–raspberry symbiosis by regulating plant immune responses.

Based on the identification of symbiotic miRNAs, quantitative real-time (qRT)-PCR analysis was used to verify the expression patterns of five representative miRNAs. Although the expression levels of each miRNA differed between the high-throughput sequencing and qRT-PCR analyses, the overall trend was consistent, indicating that the relative expression level identified by the sequencing data was reliable and suitable for further analyses (Figure 9).