agronomy

: Phytophthora root rot, caused by Phytophthora sojae ( P. sojae ), is one of the most devastating diseases limiting soybean production worldwide. microRNAs (miRNAs) play major roles in regulating plant defense against pathogens. To understand the roles of soybean miRNAs during P. sojae infection, we analyzed four small RNA libraries from two soybean germplasms before and after P. sojae isolate JS08-12 infection. The cultivar Nannong 10-1 was resistant to JS08-12, whereas the 06-070583 line was susceptible to JS08-12. In total, 528 known and 555 putative novel miRNAs in soybean were identiﬁed from 97 million reads; 74 known miRNAs and 75 novel miRNAs that might be speciﬁcally related to Nannong10-1 responses to P. sojae ; and 55 known and 43 novel miRNAs expressed before and after infection in the susceptible line 06-070583. qRT-PCR provided similar miRNA expression patterns to those obtained by the small-RNA sequencing of the four libraries. Then, the potential target genes of these differentially expressed miRNA were predicted, which. encoded transcriptional factors, resistance proteins and transporters. Finally, we focused on the targets of the three legume-speciﬁc miRNAs (gma-miR1508, gma-miR1509, and gma-miR1510) and charted the miRNA–target interactions and networks based on the published degradome data.


Introduction
Soybean (Glycine max) is one of the most cultivated crops worldwide. Its seeds, which are used as a food source, contain approximately 40% protein and 20% oil. However, like any other major crop, the production of soybean is greatly affected by pathogens such as fungi, bacteria, viruses, and oomycetes. Oomycetes are eukaryotic microorganisms that can cause devastating diseases in plants and serious economic losses in some crops [1]. Phytophthora sojae, an oomycete, causes phytophthora root rot (PRR), which is responsible for overwhelming damage to soybean production throughout the world. microRNA(miRNA), a class of non-coding endogenous small RNA, have been proven to be the key players in regulating plant growth and development, immunity against pathogens, and responses to environmental stresses at the transcriptional and post-transcriptional level by degrading corresponding mRNA or inhibiting its translation in plants [2][3][4][5]. Highthroughput technology has predicted and discovered thousands of miRNAs in many different plant species [6,7]. Most of the predicted miRNAs are involved in regulating the developmental and productive processes of plants [8]. Over the last decade or so, several studies have confirmed that miRNA also responds to infection by the pathogens of oomycetes. In 2011, Guo et al. [9] studied the miRNA involved in soybean-P. sojae interactions using microarrays and identified the differentially expressed miRNAs under P. sojae infection. Cui et al. (2017) [10] reported that gma-miR1510a/b suppressed the expression of an NB-LRR domain gene and reduced resistance to P. sojae. The overexpression of Sp-miR396a-5p made tobacco more susceptible to P. nicotianae, indicating the importance of miRNA in regulating plant disease resistance [11]. Some biotrophic oomycetes, for example, H. arabidopsidis [12,13] and P. capsici [14], are controlled by the overexpression of miRNA; in these two instances, miR393 produces resistance in Arabidopsis. Wong et al. [15] also used Illumina sequencing technology to analyze and identify the miRNAs involved in the interaction between soybean and P. sojae, and by knocking down the level of mature miR393 enhanced the susceptibility of soybean to P. sojae and drastically reduced the expression of isoflavonoid biosynthetic genes.
As the original soybean-producing country, China holds numerous soybean germplasms, and a large number of PRR-resistance germplasms have been identified in previous studies [16,17]. The P. sojae isolate JS08-12 (virulent formula: 1a, 1b, 1c,1d, 1k, 2, 3a, 3b, 3c, 4, 5, 6, and 7) was used to inoculate Nannong 10-1 and 06-070583 by hypocotyl inoculations. The soybean line Nannong 10-1 is a resistant line with 100% living seedlings, and 06-070583 is a susceptible line with 100% dead seedlings [18]. The above two germplasms are native to China. A single dominant gene, RpsJS, which is located on soybean chromosome 18 (molecular linkage group G), is responsible for the PRR resistance in Nannong 10-1. With the objective of gaining deeper insights into the P. sojae resistance mechanisms in soybean as well as understanding the molecular interactions between soybean and pathogens, we constructed four libraries from the resistant cultivar Nannong 10-1 and the susceptible line 06-070583, uninfected (control) and infected by the P. sojae isolate JS08-12. By the high-throughput sequencing of small-RNA in these libraries and by analyzing the target genes of the miRNA, the results of the present study will help elucidate the regulatory and defense mechanisms of soybean microRNA and the potential targetsunder P. sojae attack.

Plant Materials and P. sojae Isolate
The soybean germplasm Nannnong 10-1 and 06-070583 were used as plant materials and were obtained from the National Center for Soybean Improvement, China. P. sojae isolate JS08-12 was obtained from the Key Laboratory of the Monitoring and Management of Plant Diseases and Insects, Ministry of Agriculture, Nanjing Agricultural University, China.

Inoculation Method
Thirty seedlings of Nannnong 10-1 and 06-070583 were inoculated with P. sojae via a slant board assay, as described by Burnham et al. (2003) [19] with a few modifications. Thirty root-section samples were collected at 12 h after inoculation 5 mm below and above the lesion margin from each seedling with characteristic lesions. For mock inoculated plants, sections were taken from the same position as in the inoculated samples. The root samples were immediately treated with liquid nitrogen and stored at −80 • C in a refrigerator.

RNA Extraction, Library Construction, and Illumina Sequencing
For miRNA sequencing, the total RNA of 4 samples (S01-mock of Nannong 10-1; S02-JS08-12 treatment of Nannong 10-1; S03-mock of 06-070583; and S04-JS08-12 treatment of 06-070583) was extracted using TRIZol reagent (Invitrogen, Waltham, MA, USA) according to the supplier's manual. One-percent agrose gel electrophoresis was carried out to confirm the integrity and quality of total RNA. A Nanodrop 2000 spectrophotometer (IMPLEN, Westlake Village, CA, USA) was used to assess RNA purity and concentration. A Qubit 2.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) and Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) were used to measure the concentration and integrity of the RNA, respectively. The qualified RNA samples were selected for the construction of the following four small-RNA libraries using the Next Ultra Small RNA Sample Library Prep Kit for Illumina (NEB, Ipswich, NJ, USA) (Supplementary Figure S1). The Illumina HiSeq 2500 (Illumina, CA, USA) platform was employed to sequence the four small-RNA samples at Biomarker Technology Co. (Shanghai, China).

Identification of Known and Novel miRNAs
The identification of known miRNAs was used to align the clean reads to miRBase release 22 (http://www.mirbase.org, 1 May 2018) [26]. The identification of novel miRNA in soybean was performed according to the criteria for plant miRNA annotations [27]. First, the reads obtained by the sequencing were mapped onto the soybean reference genome by miRDeep2, and the secondary structure of miRNA precursors was predicted by the RNAfold program [28].

microRNA Expression Analysis
The expression of each miRNA was calculated and normalized using the TPM algorithm [29]. DESeq [30] was used to analyze the differentially expressed miRNAs. To detect the P. sojae-responsive miRNAs of each genotype, differential expression analysis was performed for the same genotype subjected to JS08-12. The criteria for differential expression were established as |log2 (fold change) | > 1 and FDR < 0.05 (Supplementary Figures S3-S5).

Validation of Differentially Expressed microRNA by Quantitative Real-Time PCR
The expression of randomly selected miRNAs was validated by stem-loop RT-PCR. Total RNA samples were treated with DNase I (Takara, Dalian, China) to remove residual genomic DNA. RNA was reverse-transcribed to cDNA, and qRT-PCR was carried out using Mir-X miRNA qRT-PCR SYBR ® Kits (Takara, Dalian, China), according to the supplier's instructions. The PCR cycling conditions were: 95 • C for 10 s, 40 cycles of 95 • C for 5 s, and 60 • C for 20 s. All reactions were repeated three times, and snoR1 served as the internal control. The relative expression level of miRNA was calculated using the 2 −∆∆CT method [31].

Functional Annotations of the Predicted Targets of the Differentially Expressed miRNA
Gene Ontology (GO) annotations, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, Clusters of Orthologous Groups of proteins (COG), Swiss-Prot, the Non-Redundant (NR) Protein Database, KOG, and the Pfam database were employed to investigate the putative biological functions of the target genes and biological processes possibly regulated by miRNAs in soybean.

High-Throughput Sequencing and Small-RNA Data Analysis
In the present study, four small RNA libraries were constructed from the roots of Nannong 10-1 and 06-070583 seedlings treated with P. sojae JS08-12 for 0 h and 12 h, respectively (S01: Nannong 10-1 mock; S02: Nannong 10-1 with P. sojae treatment; S03: 06-070583 mock; S04: 06-070583 with P. sojae treatment). Moreover, the four sRNA libraries were sequenced by Illumina technology and provided a total of 97,924,844 sRNA raw reads( Table 1). After low-quality and adapter sequences were removed, 69,325,641 reads ranging from 18 to 30 nucleotides (nt) were used for further analyses. Then, the known non-coding RNAs, including rRNA, snRNA, snoRNA, and tRNA, were annotated and removed. The numbers and proportions of the different types of small RNA are shown in Table 2.

High-Throughput Sequencing and Small-RNA Data Analysis
In the present study, four small RNA libraries were constructed from the roots of Nannong 10-1 and 06-070583 seedlings treated with P. sojae JS08-12 for 0 h and 12 h, respectively (S01: Nannong 10-1 mock; S02: Nannong 10-1 with P. sojae treatment; S03: 06-070583 mock; S04: 06-070583 with P. sojae treatment). Moreover, the four sRNA libraries were sequenced by Illumina technology and provided a total of 97,924,844 sRNA raw reads (Table 1). After low-quality and adapter sequences were removed, 69,325,641 reads ranging from 18 to 30 nucleotides (nt) were used for further analyses. Then, the known non-coding RNAs, including rRNA, snRNA, snoRNA, and tRNA, were annotated and removed. The numbers and proportions of the different types of small RNA are shown in Table 2. The ratio is equal to the separate reads divided by the total raw reads. S01 and S02 represent Nannong 10-1 seedlings treated with P. sojae for 0 h and 12 h, respectively. S03 and S04 represent 06-070583 seedlings treated with P. sojae for 0 h and 12 h, respectively.
These unannotated reads were mapped onto the G. max database using miRDeep2 computational software [25]. Figure 2 shows the locations of the mapped reads in different chromosomes of soybean. As reported in Table 3, 21, 22, and 24 nt were most abundant, which was consistent with previous studies involving small-RNA-sequencing analyses in soybean [36][37][38]. The length distribution between the two cultivars and P. sojae infection libraries was fairly similar. Small RNAs in the 21 nt and 22 nt classes showed higher abundance than those of the 24 nt class. S01 and S02 represent Nannong 10-1 seedlings treated with P. sojae for 0 h and 12 h, respectively. S03 and S04 represent 06-070583 seedlings treated with P. sojae for 0 h and 12 h, respectively. The ratio is equal to the separate reads divided by the total raw reads. S01 and S02 represent Nannong 10-1 seedlings treated with P. sojae for 0 h and 12 h, respectively. S03 and S04 represent 06-070583 seedlings treated with P. sojae for 0 h and 12 h, respectively.
These unannotated reads were mapped onto the G. max database using miRDeep2 computational software [25]. Figure 2 shows the locations of the mapped reads in different chromosomes of soybean. As reported in Table 3, 21, 22, and 24 nt were most abundant, which was consistent with previous studies involving small-RNA-sequencing analyses in soybean [36][37][38]. The length distribution between the two cultivars and P. sojae infection libraries was fairly similar. Small RNAs in the 21 nt and 22 nt classes showed higher abundance than those of the 24 nt class.

Identification of Conserved and Novel miRNA
Due to the specificity of dicer enzymes and dicer-like (DCL) enzymes, the final generation of mature miRNA is mainly concentrated in the 20 nt to 24 nt classes. In this study, 555 new miRNAs and 528 known miRNAs were predicted in the four samples. The length distributions of the known and new miRNAs identified are shown in Figure 3, vary from 19 to 25 nt, with the majority in the 20-21 nt. Since the AGO protein of the miRNA pathway preferentially associates 21 nt small RNAs starting with a 5 uridine, so a strong bias in the identification and cutting of precursor miRNA at the 5 end of the first base pair of U. The typical miRNA base proportion was obtained by analyzing the base preference of miRNA. For the first base preference of the 5 end and the base preference of the mature sequences of the total miRNAs, see Figure 3b,c. S01 and S02 represent Nannong 10-1 seedlings treated with P. sojae for 0 h and 12 h, respectively. S03 and S04 represent 06-070583 seedlings treated with P. sojae for 0 h and 12 h, respectively.

Identification of Conserved and Novel miRNA
Due to the specificity of dicer enzymes and dicer-like (DCL) enzymes, the final generation of mature miRNA is mainly concentrated in the 20 nt to 24 nt classes. In this study, 555 new miRNAs and 528 known miRNAs were predicted in the four samples. The length distributions of the known and new miRNAs identified are shown in Figure 3, vary from 19 to 25 nt, with the majority in the 20-21 nt. Since the AGO protein of the miRNA pathway preferentially associates 21 nt small RNAs starting with a 5' uridine, so a strong bias in the identification and cutting of precursor miRNA at the 5' end of the first base pair of U. The typical miRNA base proportion was obtained by analyzing the base preference of miRNA. For the first base preference of the 5' end and the base preference of the mature sequences of the total miRNAs, see Figure 3b,c.

Identification of Differentially Expressed miRNAs upon P. sojae Infection
To detect the effect of P. sojae on G. max miRNA expression, we performed a differential expression analysis between the libraries with and without P. sojae treatment. All miRNAs with more than one normalized read were analyzed by fold changes and FDR. Differential miRNA expression was considered to be indicted by an FDR lower than 0.01 and a fold change higher than two. A total of 74 conserved miRNAs (belonging to 34 families) and 75 novel miRNAs significantly changed in response to P. sojae in Nannong 10-1 ( Figure 4 and Table S1). Among them, there were 18 upregulated and 131 downregulated

Identification of Differentially Expressed miRNAs upon P. sojae Infection
To detect the effect of P. sojae on G. max miRNA expression, we performed a differential expression analysis between the libraries with and without P. sojae treatment. All miRNAs with more than one normalized read were analyzed by fold changes and FDR. Differential miRNA expression was considered to be indicted by an FDR lower than 0.01 and a fold change higher than two. A total of 74 conserved miRNAs (belonging to 34 families) and 75 novel miRNAs significantly changed in response to P. sojae in Nannong 10-1 ( Figure 4 and Table S1). Among them, there were 18 upregulated and 131 downregulated miRNAs. For the upregulated miRNAs, 12 conserved miRNAs (gma-miR1535a, gma-miR166m, gma-miR171c-3p, gma-miR171c-5p, gma-miR171n, gma-miR171o-3p, gma-miR171p, gma-miR171q, gma-miR319f, gma-miR4374a, gma-miR5037c, and gma-miR862a), belonging to seven families, and 6 novel miRNAs were identified. Sixty-two conserved miRNAs (belonging to 30 families) and sixty-nine novel miRNAs were downregulated in P. sojaetreated tissues of Nannong 10-1.
We also analyzed the differentially expressed miRNAs between 06-070583 with and without P. sojae treatment. Fifty-five conserved miRNAs and forty-three novel miRNAs were altered between the two samples (Table S2). Nine conserved miRNAs and thirty-nine novel miRNAs were upregulated in P. sojae-treated G. max roots; forty-six conserved miRNAs (belonging to 13 families) and four novel miRNAs were downregulated in P. sojae-treated G. max roots. regulated in P. sojae-treated tissues of Nannong 10-1.
We also analyzed the differentially expressed miRNAs between 06-070583 wi without P. sojae treatment. Fifty-five conserved miRNAs and forty-three novel m were altered between the two samples (Table S2). Nine conserved miRNAs and nine novel miRNAs were upregulated in P. sojae-treated G. max roots; forty-six con miRNAs (belonging to 13 families) and four novel miRNAs were downregulated sojae-treated G. max roots.

Validation of Sequencing Data by Stem-Loop qRT-PCR
To confirm the expression of the identified miRNAs and detect the dynam sponses to P. sojae infection, 11 significantly expressed miRNA families were selec stem-loop qRT-PCR to validate their expression patterns in resistant and susc plants, and the results coincided with the RNA-seq data, proving their validity. Te NAs were downregulated in the resistant cultivar Nannong10-1, and one microRN upregulated in the susceptible line 06-070583 post P. sojae inoculation (Table 4).

Validation of Sequencing Data by Stem-Loop qRT-PCR
To confirm the expression of the identified miRNAs and detect the dynamic responses to P. sojae infection, 11 significantly expressed miRNA families were selected for stem-loop qRT-PCR to validate their expression patterns in resistant and susceptible plants, and the results coincided with the RNA-seq data, proving their validity. Ten miRNAs were downregulated in the resistant cultivar Nannong10-1, and one microRNA was upregulated in the susceptible line 06-070583 post P. sojae inoculation (Table 4).

Target-Gene Prediction and Functional Annotations of the Predicted Targets of P. sojae-Responsive miRNAs in Resistant and Susceptible Soybean Germplasms
In total, 1484 target genes were predicted for 585 microRNAs (Table S3). The division of known and novel microRNAs resulted in 356 known miRNAs with 957 targeted genes and 229 novel miRNAs with 587 predicted target genes. Among the 1484 target genes, 1443 were annotated (Table S4). The numbers of annotated miRNA target genes in the  Table S5. The GO annotation showed that these target genes could participate in various cellular processes, such as "defense response", "the oxidation-reduction process", "protein phosphorylation", and "the regulation of transcription of DNA-template" (Figure 5). Based on the results of TOPGO, the two most basic "cellular component" categories were "nucleus" and "integral component of membrane"; in "biological process", the top two categories were "the regulation of transcription, DNAtemplate" and "defense response"; and in "molecular function", the top two categories were "DNA binding" and "ATP binding" in the resistant cultivar. Additionally TOPGO showed that the basic "biological processes" included "metabolic processes", "cellular processes", and "biological regulation", whereas the "molecular processes" included "binding" and "catalytic activity" and the cellular components included "cells", "cell parts", and "organelles" (Figure 6). Meanwhile, 21 KEGG pathways were classified in the resistant cultivar Nannong 10-1, and 10 KEGG pathways in the susceptible line 06-070583, and the most common pathways were "zeatin biosynthesis", "plant hormone signal transduction", "ribosome", "cyanoamino acid metabolism", "nucleotide excision repair", and "phenylpropanoid biosynthesis", which are shown below. We also found several pathways related particularly to disease response ( Figure 7) [39,40].

Discussion
Although many conserved and legume-specific miRNAs have been identified in soybean using high-throughput sequencing and bioinformatic analysis [36,[41][42][43], fewer studies have been published on microRNAs and their role in the interactions between soybean and oomycetes [9,10,44,45]. When we inoculated the resistant line Nannong 10-1 (containing the RpsJS gene) and the susceptible line 06-070583 with P. sojae JS08-12 (virulent formula; 1a, 1b, 1c, 1d, 1k, 2, 3a, 3b, 3c, 4, 5, 6 and 7) and analyzed them, we obtained 1083 known and putatively new microRNAs [18]. Upon further analysis, we obtained 74 known miRNAs and 75 novel miRNAs that were presumably involved in Nannong 10-1 responses to P. sojae, as well as 55 known and 43 novel miRNAs that were likely involved in the response of the susceptible line 06-070583 to P. sojae before and after infection. Subsequently, we predicted the target genes of these microRNAs. Finally, we focused on the legume-specific miRNAs gma-miR1508, gma-miR1509, and gma-miR1510 and used the published degradome data to analyze the resistance-related genes.
When investigating resistance-related genes, researchers focus on resistance (R) proteins, which lead to effector-triggered immunity (ETI) in plant cells upon the detection of pathogen effectors [46]. Most R proteins contain nucleotide-binding (NB) and proteinbinding leucine-rich repeat (LRR) domains. An examination of these NB-LRR proteins shows the presence of either a TIR (Toll interleukin 1 receptor) or CC (coiled coil) domain, which detect and recognize biotic effectors, i.e., fungi, oomycetes, bacteria, and viruses, to stimulate the plant defense system [46,47]. Remarkably, other 22 nt miRNAs targeting the NBS-LRR class of R genes have been found to be abundant and diverse in legumes and Solanum species [48,49]. For example, secondary siRNAs are produced via RDR6 and DCL4 by miR482, miR2109, and miR1507, which perform a silencing effect [50,51]. Hence, to target and regulate the rapidly evolving R genes, diversifying the secondary siRNAs is, in all probability, the most effective way to maximize their response [47]. Additional examples of these 22 nt miRNAs include: miR2118 (the passenger strand of soybean miR482), which targets the conserved P-loop motif of TIR-NBS-LRR; miR2109, which targets the TIR-1 motif of TIR-NBS-LRR; and miR1507, which targets the kinase-2 motif of CC-NBS-LRR [50]. In another experiment, two miRNAs (miR2109 and miR1507) were produced when soybean was infected with active P. sojae, but not when it was infected with heat-inactivated P. sojae, hinting at these miRNAs' possible role in ETI; the up-and downregulation of these miRNAs proved this hypothesis to be true [15]. To confirm this phenomenon, when pathogen effectors suppressed RNAi, the upregulation of miRNA targets (which, in this case, were the R genes) provided an abundance of R proteins to fight the pathogen [49]. Hence, our identification of the targeted genes gma-miR1508a ( Figure 8) and gma-miR1510b-3p ( Figure 9) through degradome analysis proved the involvement of said microRNAs in the defense against P. sojae. A previous study [15] showed that GmTNL16, encoding NBS-LRR-type proteins, was targeted by gma-miR1510; overexpressed GmTNL16 and silenced gma-miR1510 could improve the resistance to P. sojae in soybean hairy roots. Furthermore, JA-and SA-pathway-related genes, including JAZ, COI1, TGA, and PR genes, responded to P. sojae treatment.
microRNA are important regulators of plant hormone signals. ABA is a key plant hormone that plays an important role under multiple stress conditions by mediating the expression of stress-related genes and inducing stomatal closure, acting as the basic responder to environmental changes [52,53]. Several key genes in the microRNA biogenesis pathway, including HYL1, DCL1, HEN1, SE, and HASTY, were impaired with ABA-hypersensitive mutants during germination [54,55], suggesting that microRNAs are involved in ABA signaling. Moreover, miRNA159a was proven to be involved in ABA signal transduction [56]. The downregulation of gma-miR159d in cyanoamino acid metabolism was observed in the differential expression analysis of the resistant cultivar Nannong 10-1 during our data mining study. This suggested that gma-miR159d is suppressed when attacked by P. sojae, which in turn kick starts ABA production.
Many growth-related pathways, such as leaf and root architecture and vascular development, are modulated by a hormone named auxin, which is positively regulated by TIR (transport inhibitor response 1), as it promotes Aux/IAA proteins through ubiquitination. miR393 targets TIR, hence controlling auxin production. An increased level of miR393 would downregulate auxin signals and reduce plant growth. In the resistant cultivar Nannong 10-1, the downregulation of TIR1, which was targeted by gma-miR393 family members, strongly suggested the downregulation of ubiquitin mediated proteolysis, inhibiting cellular development and plant growth (Supplementary Figure S6a). The alternate responsive expression pattern suggested the instigation of a wide range of regulatory mechanisms by the miR393 family (Table S1).
gma-miR172k and gma-miR172f are inhibited in the resistant cultivar Nannong 10-1, and the target genes maybe regulate the ABF signals (Table S1). This suggested a positive correlation with stomatal closure upon P. sojae attack. When we investigated the susceptible line, a novel microRNA was found to be responsible for the downregulation of PP2C, suggesting an innate immunity signal in susceptible soybean (Supplementary Figure S6b).
We know that miR396a overexpression leads to a lowered stomatal density in Arabidopsis, proving to be an important factor in stress responses [52,57]. miR166 is crucial for cell development, as it regulates the class-3 homeodomain-leucine zipper transcription factors, which partially control lateral root development, auxiliary meristem initiation, and leaf polarity [11,[53][54][55]. Both of these microRNAs were downregulated in the resistant line upon P. sojae attack, suggesting innate stress response activation (Supplementary Table S1).
Cyanide, a nitrogen-rich compound that is effectively toxic to mammals, is naturally produced in plants, algae, fungi, and bacteria. The processing of cyanide is achieved through cyanogenesis, whereby cyanide is degraded and nitrogen is released, which is then used in growth. Cyanogenesis in a wide range of plants constitutes a chemical defense system against herbivores and pathogens [56,58]. The cyanoamino acid metabolism pathway, involving six novel microRNAs downregulated in the resistant cultivar and upregulated in the susceptible line, showed a negative correlation with pathogen-responsive genes in the Nannong 10-1 cultivar (Supplementary Figure S7a,b). These microRNAs showed the same pattern in starch and sucrose metabolism along with phenylpropanoid biosynthesis (Supplementary Figure S8a,b), suggesting the wide range of deactivated pathways controlled by these six novel microRNAs.
Furthermore, the downregulation of MKK1/2 by a novel microRNA was observed in the KEGG analysis, suggesting that through the downregulation of MKK1/2 as part of the plant-pathogen interaction, the signal was interrupted to deactivate pathogenresponsive genes and WRKY33, which in turn controls the oomycete defense-related genes (Supplementary Figure S9a,b).

Conclusions
Our results revealed the basic structure of microRNAs involved in plant metabolism under P. sojae attack and provided the functional annotation of eight novel microRNAs involved in soybean metabolism. Based on these findings, we could conclude that miRNAs are not just regulators but are key players in defense against P. sojae. This could be a breakthrough for understanding new aspects of plant defense mechanisms, as this draft could provide a framework for function identification of miRNAs and the development of new methods for improving P. sojae resistance in soybean.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/agronomy12122922/s1, Table S1: Differentially expressed miRNAs in Nannong 10-1 in response to P. sojae, Table S2: Differentially expressed miRNAs in 06-070583 in response to P. sojae, Table S3: miRNA target gene information, Table S4: Differential expression of miRNA target gene annotation, Table S5: Samples of different miRNA target gene annotation results, Figure S1: sRNA flowchart of sequencing information analysis, Figure S2: Homologue species of miRNA families found in different plants, Figure S3: RPKM and TPM distribution of four libraries along with their Pearson correlation, Figure S4: Correlation plot between the different samples, Figure S5: Differentially expressed miRNA volcano map, Figure S6a-b: Plant-hormone signaling pathway in S01 vs. S02 along with their DEGs, Figure S7: Cyanoamino acid metabolism showing the downregulated DEGs in S01 vs. S02 and S03 vs. S04, Figure S8: KEGG pathway of phenylpropanoid biosynthesis for differentially expressed miRNA target genes in S01 vs. S02 and S03 vs. S04, Figure S9: Plant-pathogen interaction pathway in S01 vs. S02 and S03 vs. S04.