Identification of MicroRNAs That Respond to Soybean Cyst Nematode Infection in Early Stages in Resistant and Susceptible Soybean Cultivars

Soybean cyst nematode (SCN) causes heavy losses to soybean yield. In order to investigate the roles of soybean miRNAs during the early stages of infection (1 and 5 dpi), 24 small RNA libraries were constructed from SCN resistant cultivar Huipizhi (HPZ) and the susceptible Williams 82 (W82) cultivar for high-throughput sequencing. By sequencing the small RNA libraries, a total of 634 known miRNAs were identified, and 252 novel miRNAs were predicted. Altogether, 14 known miRNAs belonging to 13 families, and 26 novel miRNAs were differentially expressed and may respond to SCN infection in HPZ and W82. Similar expression results were also confirmed by qRT-PCR. Further analysis of the biological processes that these potential target genes of differentially expressed miRNAs regulate found that they may be strongly related to plant–pathogen interactions. Overall, soybean miRNAs experience profound changes in early stages of SCN infection in both HPZ and W82. The findings of this study can provide insight into miRNAome changes in both HPZ and W82 at the early stages of infection, and may provide a stepping stone for future SCN management.


Introduction
Soybean [Glycine max (L.) Merr.] is rich in oils and proteins, which are in high demand all over the world. Soybean yield is negatively impacted by several soybean pathogens, such as soybean mosaic virus, Pseudomonas syringae, Phytophthora sojae, Phakopsora pachyrhizi, and Heterodera glycines [1]. Among these, Heterodera glycines, also known as soybean cyst nematode (SCN), causes heavy losses in soybean production of over 120 million USD in China [2] and 1.2 billion USD in the United States [3], making it the most damaging soybean pest worldwide. Therefore, how to restrict the soybean yield losses caused by SCN is a worldwide problem wherever soybean is grown.
To date, SCN is primarily managed by the use of resistant cultivars. In North America, the Peking and PI88788 cultivars have been the predominant resistant sources [4]. Although soybean cultivars conferring resistance to SCN can limit production losses, they can also cause virulence shifts in SCN populations, resulting in SCN with the ability to reproduce on resistant cultivars [5]. Thus, understanding the mechanisms behind soybean-SCN interactions is key in preventing soybean damage from SCN. Although great advances have been made in terms of resistant genes in the form of the changes in soybean miRNA responses to SCN in the early stages of infection need to be investigated, and it will be valuable to have a detailed understanding of soybean-SCN interactions at the miRNA level, due to the fact that soybean miRNAs play a key role in soybean disease defenses [27]. Here, we sequenced 24 miRNA libraries from an elite resistant soybean cultivar (Huipizhi) that confers resistance to diverse SCN races (race 1, 3, 4, 5, 7, and 14) [28] and from a susceptible cultivar (Williams 82) at 1 and 5 dpi of SCN infection.

HPZ Was Highly Resistant to SCN During the Early Stages of Infection/HPZ Represses the Development of SCN and Syncytium
The HPZ soybean cultivar has a broad resistant spectrum to several SCN races (race 1, 3, 4, 5, 7, and 14) [28]. To investigate the infection and development of SCN (race 3) in the roots of the HPZ and W82 cultivars during the early stages of infection (1 and 5 dpi), we inoculated each HPZ and W82 seedling with 2000 SCN J2s, and the infection rate and development state were recorded at 1 and 5 dpi. There were no significant differences in the total number of SCN juveniles in the roots of HPZ and W82 at 1 dpi ( Figure 1A-C). However, nearly 50% of SCN juveniles in W82 had developed to the swollen stage at 5 dpi, compared to only 10 % in HPZ ( Figure 1D-F). The fact that development of SCN was delayed in HPZ may be related to the establishment and development of syncytium. Therefore, we also measured the size of syncytium that SCN induced in the roots of HPZ and W82 at 5 dpi. The average syncytium size in W82 was 0.028 mm 2 , which was 1.5 times greater than the HPZ average (0.019 mm 2 ) at 5 dpi (Figure 2), and this difference was significant (p < 0.05). This result also suggests that resistant soybean cultivars do not impact the penetration of SCN, but that they instead degenerate and restrict the formation and development of syncytia using multiple mechanisms in the early stages of soybean-SCN interactions [29][30][31]. be investigated, and it will be valuable to have a detailed understanding of soybean-SCN interactions at the miRNA level, due to the fact that soybean miRNAs play a key role in soybean disease defenses [27]. Here, we sequenced 24 miRNA libraries from an elite resistant soybean cultivar (Huipizhi) that confers resistance to diverse SCN races (race 1, 3, 4, 5, 7, and 14) [28] and from a susceptible cultivar (Williams 82) at 1 and 5 dpi of SCN infection.

HPZ Was Highly Resistant to SCN During the Early Stages of Infection/HPZ Represses the Development of SCN and Syncytium
The HPZ soybean cultivar has a broad resistant spectrum to several SCN races (race 1, 3, 4, 5, 7, and 14) [28]. To investigate the infection and development of SCN (race 3) in the roots of the HPZ and W82 cultivars during the early stages of infection (1 and 5 dpi), we inoculated each HPZ and W82 seedling with 2000 SCN J2s, and the infection rate and development state were recorded at 1 and 5 dpi. There were no significant differences in the total number of SCN juveniles in the roots of HPZ and W82 at 1 dpi ( Figure 1A-C). However, nearly 50% of SCN juveniles in W82 had developed to the swollen stage at 5 dpi, compared to only 10 % in HPZ ( Figure 1D-F). The fact that development of SCN was delayed in HPZ may be related to the establishment and development of syncytium. Therefore, we also measured the size of syncytium that SCN induced in the roots of HPZ and W82 at 5 dpi. The average syncytium size in W82 was 0.028 mm 2 , which was 1.5 times greater than the HPZ average (0.019 mm 2 ) at 5 dpi (Figure 2), and this difference was significant (p < 0.05). This result also suggests that resistant soybean cultivars do not impact the penetration of SCN, but that they instead degenerate and restrict the formation and development of syncytia using multiple mechanisms in the early stages of soybean-SCN interactions [29][30][31]. . SCN (soybean cyst nematode) juveniles; infection rate and development state assay in HPZ and W82. Each soybean seedling was inoculated with 2000 second-stage SCN juveniles, and the infection rate was assayed at 1 dpi in W82 and HPZ (A-C), A: second-stage SCN juveniles infected the root of W82, bar = 500 μm, B: second-stage SCN juveniles infected the root of HPZ, bar = 500 μm, C: total number of second-stage SCN juveniles infected the root of W82 and HPZ, values are the mean of the number of SCN, bar, standard error. The development state was assayed at 5 dpi (D-F), D: SCN juveniles developed into the swollen stage in W82, bar =200 μm, (E) the development of SCN juveniles was delayed in HPZ, bar = 200 μm, (F) total number of SCN different state juveniles in W82 and HPZ at 5 dpi, values are the mean of the number of SCN, bar, standard error. Pictures of SCN juveniles in the soybean root were taken by an OLYMPUS DP80 light microscope (Olympus, Tokyo, Japan). Figure 1. SCN (soybean cyst nematode) juveniles; infection rate and development state assay in HPZ and W82. Each soybean seedling was inoculated with 2000 second-stage SCN juveniles, and the infection rate was assayed at 1 dpi in W82 and HPZ (A-C), A: second-stage SCN juveniles infected the root of W82, bar = 500 µm, B: second-stage SCN juveniles infected the root of HPZ, bar = 500 µm, C: total number of second-stage SCN juveniles infected the root of W82 and HPZ, values are the mean of the number of SCN, bar, standard error. The development state was assayed at 5 dpi (D-F), D: SCN juveniles developed into the swollen stage in W82, bar =200 µm, (E) the development of SCN juveniles was delayed in HPZ, bar = 200 µm, (F) total number of SCN different state juveniles in W82 and HPZ at 5 dpi, values are the mean of the number of SCN, bar, standard error. Pictures of SCN juveniles in the soybean root were taken by an OLYMPUS DP80 light microscope (Olympus, Tokyo, Japan). Nematode infection and development data were checked for normality with the Kolmogorov-Smirnov or Shapiro-Wilk test, then the t-test method was used to analyze these above data with SPSS software. Each soybean cultivar contained 10 replicates, ns: not significant, p < 0.05, p < 0.01 means a significant difference was found between each comparison.
Nematode infection and development data were checked for normality with the Kolmogorov-Smirnov or Shapiro-Wilk test, then the t-test method was used to analyze these above data with SPSS software. Each soybean cultivar contained 10 replicates, ns: not significant, p < 0.05, p < 0.01 means a significant difference was found between each comparison.

Figure 2.
Histological observation of syncytium size in HPZ and W82 at 5 dpi. (A) uninfected W82 root, bar = 100 μm, (B) SCN-infected W82 root, the syncytium induced by SCN was labeled with "S", bar = 50 μm, (C) uninfected HPZ root, bar = 100 μm, (D) SCN-infected HPZ root, the syncytium induced by SCN was labeled with "S", bar = 50 μm, (E) average size of syncytium in W82 and HPZ values are the mean syncytium size, bar, standard error. Pictures of syncytia were taken by an OLYMPUS DP80 light microscope (Olympus, Japan), and the size of the syncytia was measured with Cellsens standard software. Syncytium size data were checked for normality with the Kolmogorov-Smirnov or Shapiro-Wilk test, then the t-test method was used to analyze these above data with SPSS software. Size values of the five syncytia were used for t-test, p < 0.05 means a significant difference was found between W82 and HPZ.

Small RNA Sequenced from the 24 Libraries Constructed from HPZ and W82 Roots
Despite HPZ conferring strong resistance to SCN in the early stages of SCN parasitism, until now the changes and roles of miRNAs during this important stage remain to be elicited. Toward this goal, we constructed 24 small RNA libraries to sequence using an Illumina novaseq6000 platform, and a workflow diagram of library construction is shown in Figure 3. These 24 libraries provided a total of 1,020,975,683 sRNA raw reads, and the Q30 of all the libraries was ≥94%. After data quality control, a total of 643,468,338 sRNA clean reads were obtained, and each library contained a minimum of 19,622,055 clean reads (Table 1). Then, these clean reads were used for sRNA annotation and no less than 17.07% of the clean reads in each sample were using to detect known and novel miRNAs (Material S1).

Figure 2.
Histological observation of syncytium size in HPZ and W82 at 5 dpi. (A) uninfected W82 root, bar = 100 µm, (B) SCN-infected W82 root, the syncytium induced by SCN was labeled with "S", bar = 50 µm, (C) uninfected HPZ root, bar = 100 µm, (D) SCN-infected HPZ root, the syncytium induced by SCN was labeled with "S", bar = 50 µm, (E) average size of syncytium in W82 and HPZ values are the mean syncytium size, bar, standard error. Pictures of syncytia were taken by an OLYMPUS DP80 light microscope (Olympus, Japan), and the size of the syncytia was measured with Cellsens standard software. Syncytium size data were checked for normality with the Kolmogorov-Smirnov or Shapiro-Wilk test, then the t-test method was used to analyze these above data with SPSS software. Size values of the five syncytia were used for t-test, p < 0.05 means a significant difference was found between W82 and HPZ.

Small RNA Sequenced from the 24 Libraries Constructed from HPZ and W82 Roots
Despite HPZ conferring strong resistance to SCN in the early stages of SCN parasitism, until now the changes and roles of miRNAs during this important stage remain to be elicited. Toward this goal, we constructed 24 small RNA libraries to sequence using an Illumina novaseq6000 platform, and a workflow diagram of library construction is shown in Figure 3. These 24 libraries provided a total of 1,020,975,683 sRNA raw reads, and the Q30 of all the libraries was ≥94%. After data quality control, a total of 643,468,338 sRNA clean reads were obtained, and each library contained a minimum of 19,622,055 clean reads (Table 1). Then, these clean reads were used for sRNA annotation and no less than 17.07% of the clean reads in each sample were using to detect known and novel miRNAs (Material S1).  In total, 1 mL 0.1% water-agar mixture containing approximately 2000 SCN J2s (second stage juveniles) or water-agar mixture alone as the treatment or control, respectively, were added to the root system of soybean seedlings. Three soybean seedlings were pooled as one sample, each sample contained three biological replicates, and these samples were collected at 1 and 5 dpi, respectively. In total, 24 samples used for small RNA libraries construction were prepared, namely, 24 samples = 2 (treatments) × 2 (time-points) × 2 (cultivars) × 3 (biological replicates). In total, 1 mL 0.1% water-agar mixture containing approximately 2000 SCN J2s (second stage juveniles) or water-agar mixture alone as the treatment or control, respectively, were added to the root system of soybean seedlings. Three soybean seedlings were pooled as one sample, each sample contained three biological replicates, and these samples were collected at 1 and 5 dpi, respectively. In total, 24 samples used for small RNA libraries construction were prepared, namely, 24 samples = 2 (treatments) × 2 (time-points) × 2 (cultivars) × 3 (biological replicates).

Analysis of Differences in miRNA Expression Between HPZ and W82
To identify miRNAs responsive to SCN infection in both HPZ and W82 at the early stages of infection (1 and 5 dpi), miRNA differential expression (DE) profiling was completed between control groups (W1C, W5C, H1C, H5C) and SCN-treated groups (W1N, W5N, H1N, H5N) at different time points (1 and 5 dpi) in both HPZ and W82.
A total of 14 known miRNAs belonging to 13 families, and 26 novel miRNAs were differentially expressed in the 24 libraries (Material S3). Eleven miRNAs were differentially expressed between the susceptible cultivar W82 control library and the SCN-treated library at 1 dpi, while only three miRNAs (gma-miR3522, gma-miR408a-3p, and novel_miR_178) were upregulated at 5 dpi. Between the HPZ control library and the HPZ SCN-treated library, seven miRNAs were up-or downregulated at 1 dpi, and only three novel miRNAs (novel_miR_106, novel_miR_13, and novel_miR_34) were upregulated. We also conducted DE miRNAs analysis of the same cultivar infected by SCN at different time points in order to identify the soybean miRNAs that changed with SCN development. Eight miRNAs were downregulated in W5N compared to W1N, including conserved miRNAs, such as gma-miR408a, gam-miR171, gma-miR159, and gma-miR398d, while three miRNAs were upregulated, including gma-miR5037c and novel_miR_190. Six miRNAs were upregulated in H5N compared to H1N, including gma-miR5225, while three miRNAs were downregulated, such as miR171b-5p. To find the miRNAs that were differentially expressed between the resistant cultivar HPZ and the susceptible cultivar W82, we compared the SCN-infected libraries of the two cultivars. We found that three miRNAs were upregulated (novel_miR95, novel_miR190, and novel_miR136), and seven miRNAs were downregulated, including gma-miR408a, novel_miR120, and novel_miR214, in W1N compared to H1N. Nine miRNAs were downregulated and no miRNAs were upregulated in W5N compared to H5N (Table 3 and Material S3). After DE analysis of miRNAs in all DE sets, a hierarchical clustering analysis was performed on all DE miRNAs using TBtools software [34] based on the value of log10 (TPM1e-6). miRNAs with the same or similar expression behaviors were clustered. As shown in Figure 4, gma-408a seemed to respond strongly to SCN infection in the W82 susceptible cultivar at 1 and 5 dpi, whereas gma-4994-5p responded to SCN infection in the HPZ resistant cultivar at 1 and 5 dpi. We could deduce that gam-398d, gma-159a-5p, and novel_miR_34 responded to SCN infection and were highly expressed ( Figure 4).  To validate the sequencing data and our DE analysis of miRNAs, we conducted qRT-PCR using the same 24 RNA samples. In each DE set, three differentially expressed miRNAs found in silico were randomly selected for qRT-PCR experiments. The results of the qRT-PCR showed that the expressions of most of the miRNAs (19 of 24 selected miRNAs) agreed with the sequencing result. For instance, gam-miR408a-5p expression was 4.1 times higher in W1N compared to W1C according To validate the sequencing data and our DE analysis of miRNAs, we conducted qRT-PCR using the same 24 RNA samples. In each DE set, three differentially expressed miRNAs found in silico were randomly selected for qRT-PCR experiments. The results of the qRT-PCR showed that the expressions of most of the miRNAs (19 of 24 selected miRNAs) agreed with the sequencing result. For instance, gam-miR408a-5p expression was 4.1 times higher in W1N compared to W1C according to the NGS result (Material S4), and qRT-PCR found that its expression was 6.1 times higher ( Figure 5). The same consistent expression patterns could be found for other comparisons, like gma-miR3522 expression between H1C and H1N, and gma-miR398d expression between W5N and H5N (Material S4 and Figure 5).

Target Prediction and Function Annotation
miRNAs exert their regulatory roles by interaction with transcripts of their target genes. In this way, miRNAs regulate diverse biological processes within the plant. To understand the biological changes that miRNAs regulate in soybean in response to SCN infection in both HPZ and W82, we used TargetFinder (v1.6) to predict the potential target genes of the identified DE miRNAs. A total of 8487 potential target genes were predicted (Material S4), before they were annotated by the Blast to Three miRNAs in each comparison were randomly selected for qRT-PCR validation, X axis stands for relative expression of miRNAs, NGS, the expression pattern of miRNAs obtained in high-throughput sequencing, qRT-PCR, and the miRNA expression pattern of miRNAs validated by qRT-PCR.

Target Prediction and Function Annotation
miRNAs exert their regulatory roles by interaction with transcripts of their target genes. In this way, miRNAs regulate diverse biological processes within the plant. To understand the biological changes that miRNAs regulate in soybean in response to SCN infection in both HPZ and W82, we used TargetFinder (v1.6) to predict the potential target genes of the identified DE miRNAs. A total of 8487 potential target genes were predicted (Material S4), before they were annotated by the Blast to the NR, Swiss-Prot, GO, KEGG, and Pfam databases. We noticed that the target genes of DE miRNAs were strongly linked to plant defense processes in both W82 and HPZ. For example, gma-miR408a regulated several genes that function as xylanase inhibitors, protein tyrosine kinases, and ubiquitin-conjugating enzymes, and gma-miR159a regulated genes with GAMYB domains or NB-ARC domains, which program cell death and pathogen responses. In addition, gma-miR398d regulated genes with copper/zinc superoxide dismutase (SODC) functions, gma-miR862b regulated genes encoding serine hydroxymethyltransferase, which is a key gene in Rhg loci, and gma-miR4415 was predicted to regulate several multicopper oxidase genes. Several novel miRNAs were predicted to regulate signal transduction-and phytohormone-related genes. For instance, novel_miR_106 regulated genes related to K+ potassium transporters, novel_miR_13 regulated genes with leucine-rich repeat domains, novel_miR_34 regulated genes encoding auxin response factors, and novel_miR_178 regulated genes with AP2 ethylene-responsive transcription factor domains (Material S4). These miRNAs may play essential roles in the interactions between SCN and both the HPZ resistant cultivar and the susceptible W82 cultivar (Table 4 and Material S4).
In order to gain insight into the biological processes that are influenced by the SCN infection in both HPZ and W82, we conducted KOG, GO, and KEGG pathway classification for the biological processes influenced by these potential target genes using a BMK cloud analysis platform. As expected, we found that signal transduction mechanisms, carbohydrate transport, metabolism, and defense mechanisms occurred very frequently in KOG classification ( Figure 6). In GO classification, immune system processes and responses to stimulus detoxification were found at high percentages in biological processes. Antioxidant activity, nutrient reservoir activity, and signal transducer activity occurred frequently in the molecular function (Figure 7). In KEGG pathway classification, 4.58% of genes were involved in plant-pathogen interaction in the organism systems, 5.30% of genes participated in amino acid biosynthesis, 4.03% were involved in starch and sucrose metabolism, and 3.54% played roles in phenylpropanoid biosynthesis (Figure 8).

Discussion
Plant miRNAs play an active role in regulating diverse biological processes, including responses to biotic and abiotic stresses, and many recent studies have provided evidences of this. Jiang et al. Figure 8. KEGG pathway classification of target genes of differentially expressed miRNAs. The left Y axis is the name of the KEGG metabolic pathway, the right Y axis is the biological process that these genes participated, and the X axis is the number and the percentage of target genes.

Discussion
Plant miRNAs play an active role in regulating diverse biological processes, including responses to biotic and abiotic stresses, and many recent studies have provided evidences of this. Jiang et al. found that miR482b negatively impacted tomato resistance to Phytophthora infestans, with plants overexpressing miR482b displaying more serious disease symptoms [35]. Dai and coworkers proved that OsmiR396 negatively impacted rice resistance to brown planthoppers [36], and Ding et al. found that overexpression of miR160 could increase the sensitivity of cotton to high temperature stress [37]. Plant miRNAs also strongly respond to parasitism by plant parasitic nematodes. In Arabidopsis, Hewezi and coworkers studied the relationship between Arabidopsis miRNAs and parasitism by the beet cyst nematode (Heterodera schachtii) by sequencing the miRNA libraries constructed from the syncytia formed in Arabidopsis roots using LCM (laser capture microdissection) technology. Their results identified the miRNAs that were differentially expressed compared to control roots, and syncytia were predicted to regulate genes coding for transcription factors that regulate hormonal responses, cellular specialization, or cell expansion [38]. Several miRNAome sequencing experiments have since been conducted to investigate how plant miRNAs change when the plant is infected by parasitic nematodes, including studies on Arabidopsis infected with the root knot nematode M. javanica [39], Arabidopsis infected with M. incognita [22], tomato infected with the potato cyst nematode Globodera rostochiensis [40], and tomato infected with M. incognita [41]. Soybean cyst nematodes are one of the most devastating pathogens of soybean, casing heavy losses to soybean production worldwide [2,3]. To understand the function of miRNAs displayed during soybean SCN infection, three miRNAome changes experiments were conducted by SCN researchers. Li et al. investigated the miRNAs involved in infection by SCN race 3 using a resistant cultivar (HB) and a susceptible cultivar (L10). They found that 101 miRNAs were SCN responsive and that 20 miRNAs differed in their expression patterns between resistant and susceptible cultivars. Xu and coworkers sequenced miRNAome changes between soybean roots infected with SCN race 4 and control roots from resistant and susceptible cultivars. Their results showed that seven miRNAs were differentially expressed after SCN infection [25]. These previous studies provide us with new insight into how miRNAs respond to SCN infection; however, due to these studies lacking biological replicates, there may be some bias caused by different expression analysis methods. Recently, Tian et al. used two soybean cultivars (SCN susceptible KS4607 and SCN HG Type 7 resistant KS4313N) that were grown under SCN-infested and non-infested soil at two different time points (SCN feeding establishment and egg production). In their study, each treatment contained three replicates, and they identified 60 miRNAs that were specifically related to different resistant cultivars when infected by SCN [26]. To understand the miRNA changes at the early stages of infection, in our present study, we constructed 24 small libraries from W82 control plants, W82 plants infected by SCN, HPZ control plants, and HPZ plants infected by SCN at the early stages of SCN infection (1 and 5 dpi; namely W1C, W1N, H1C, H1N, W5C, W5N, H5C, and H5N, the workflow is shown in Figure 3). Each treatment consisted of three replicates in order to achieve high confidence in our high-throughput sequencing results because of the different expressions of miRNAs among each replicate, as shown in Supplementary Materials S2 and S3. Using the Illumina Novaseq6000, we obtained a total of 643,468,338 sRNA clean reads after removing the low-quality reads and adapters, which is greater than the 13,835,079 clean reads obtained in Xu et al. [25], the 44,515,792 reads reported in Li et al. [24], and the 58,628,603 reads in Tian et al. [26].
In this study, owing to the large number of reads obtained by Illumina Novaseq6000, a total of 634 known miRNAs were identified and 252 novel miRNAs were predicted. The most important aspect of high-throughput sequencing of miRNA changes is different expression analysis. We normalized the expression level of each miRNA using TPM, and used DEseq2 software to find miRNAs that were differentially expressed by setting the |log2(FC)| ≥ 1.00; FDR (False Discovery Rate) ≤ 0.01. We found 14 known miRNAs belonging to 13 families, and 26 novel miRNAs were differentially expressed in these 24 libraries (Supplementary Material S3). Our results differed greatly from those of other studies in terms of the number of DE miRNAs or families of miRNAs. This may be due to the specific soybean cultivars used for sequencing, the time points at which we collected the samples, and the software used for analysis [24][25][26]. In this study, there were more DE miRNAs in the susceptible soybean cultivar and the resistant cultivar at 1 dpi compared to 5 dpi. No miRNAs were downregulated, but only three miRNAs we found to be differentially expressed in comparison between W5C and W5N and between H5C and H5N, and the fold-changes of all DE miRNAs at 5 dpi were less than 1.5 (Table 3 and Material S4). This is an interesting phenomenon in our sequencing results, and our qRT-PCR confirmed the expressions of most of the DE miRNAs between W5C and W5N and between H5C and  H5N (Figure 7). The fact that more miRNAs were differentially expressed in both HPZ and W82 at 1 dpi compared to 5 dpi could be explained by extensive tuning of the host's response to nematode infection [42]. Koter et al. reported a sharp increase in host miRNA synthesis in the first days after nematode infection in tomato [40]. However, studies have shown that more miRNAs are upregulated in the late stages of nematode infection in soybean [26] and in Arabidopsis [38]. We considered that these differences may come from the type of plant, specific nematode, and the time points measured in these studies. Therefore, differences between sequencing experiments are a normal phenomenon. miRNAs exert control by influencing the expressions of target genes. In order to further understand which genes of soybean cultivar HPZ and W82 were regulated by miRNAs during soybean root SCN infection, we used TargetFinder (v1.6) to predict the potential target genes of the identified DE miRNAs. The legume-specific gma-miR3522, which was upregulated in W82 at 1 and 5 dpi and in HPZ at 1 dpi, was predicted to negatively control several genes encoding polyphenol oxidase (Material S4), which enhances the resistance of tomato to Pseudomonas syringae [43]. Of the other legume-specific miRNAs, gma-miR5037c was upregulated in W5N compared to W1N, and targets genes encoding GMC oxidoreductase, which is reportedly involved in the biocontrol activities of several plant pathogenic fungi [44]. gma-miR5225 was differentially expressed in the H1N vs. H5N comparison, and targets genes coding for protein tyrosine kinase, which plays essential roles in plant disease resistance, stress responses, and cell death in Arabidopsis [44]. Little is known about the function of these legume-specific miRNAs in soybean, and their roles remain to be elucidated, especially in terms of how they impact the interactions between soybean and SCN.
The function of conserved miRNA is relatively well understood, as they are highly conserved in function and structure. In our study, conserved miRNA gma-miR408a was found to respond to SCN infection in W82 at 1 and 5 dpi, and in HPZ at 1 dpi. On the other hand, Tian et al. found that gma-miR408a only responded to SCN infection in a susceptible cultivar [26]. The potential target genes of gma-miR408a encode multicopper oxidase, which contributes to plant resistance to several plant pathogens [45]. Furthermore, gma-miR159a and gma-miR398d were both differentially expressed in the comparison between W1N and W5N. gma-miR159 was reported to play a role in plant resistance to P. syringae and regulates the GAMYB transcript factor [46,47], while gma-miR398d targets genes coding for copper/zinc superoxide dismutase (SODC), which participates in plant responses to heat stress and barley powdery mildew fungus [48,49]. Downregulation of gma-miR398 induces the plant to produce superoxide dismutase 1 (CSD1) and superoxide dismutase 2 (CSD2), which improve tolerance to oxidative stress in Arabidopsis [50]. In our study, 252 novel miRNAs were predicted, and different expression analysis showed that several novel miRNAs were differentially expressed between W82 and HPZ, indicating that they may play important roles in responding to SCN infection. For example, novel_miR_178 was upregulated in W5N compared to W5C, which targets AP2/ EREBP (APETALA2/ethylene-responsive element binding protein) transcripts. It has been reported that overexpression for AP2/EREBP enhanced resistance to Phytophthora parasitica in tobacco [51]. Novel_miR_106 was upregulated in H5N compared to H5C, which targets several genes coding for leucine-rich repeat domains. This domain belongs to typical resistant genes (R gene), which play a positive role in inhibiting the invasion of plant pathogens [52]. The target genes of the DE miRNAs in every comparison mentioned above appear to be potential soybean miRNA candidates that are responsive to SCN infection in both HPZ and W82 in the early stages of infection.
In present study, we conducted KOG, GO, and KEGG pathway classification of these predicted target genes to gain insight into which soybean biological processes were altered by SCN infection. We found that the immune system, stimulus detoxification, plant-pathogen interaction, and phenylpropanoid biosynthesis occupied a high percentage in the KOG, GO, and KEGG pathway classifications. These terms were highly related to plant defense, for example, the phenylpropanoid biosynthesis term was reported to be involved in plant resistance to disease [53]. Further, nutrient reservoir activity, signal transducer activity, biosynthesis of amino acids, and starch and sucrose metabolism also constituted high percentages in these classifications, indicating that the soybean metabolism was dramatically altered in response to SCN invasion. As Siddique and Grundler reviewed, in nematode feeding sites, nutrients like amino acids experience profound changes, and sucrose and starch are increased during nematode development [54]. These metabolic changes to the target genes of the DE miRNAs may increase our understanding of the response of both HPZ and W82 to SCN.
As has been mentioned in this paper, several experiments were conducted to understand host plant miRNAome changes during the plant parasitism by nematodes. These studies focused on the miRNAome changes caused by CN (cyst nematode, Heterodera spp.) [24][25][26]40] and RKN (root knot nematode, Meloidogyne spp.) [22,38,39,41]. To date, the function of miR159, miR172, and miR319 in the interaction between plants and plant parasitic nematodes are the few miRNA functions that have been deeply elucidated [22,55,56], and many questions remain surrounding the roles of miRNAs in plant parasitism by nematodes. For instance, cross-kingdom RNA silencing probably also occurs during interactions between plants and plant parasitic nematodes [57]. In other plant species, recent studied have indicated that plant miRNAs can silence the transcripts of plant pathogens that are important for virulence [21], and that miRNAs from the plant pathogen also can silence genes of their hosts that are important for immunity [58]. The question remains as to whether soybean miRNAs silence SCN genes that are important for parasitism. Combining high-throughput sequencing technologies and novel technologies like short tandem target mimic (STTM), which has already been applied to investigate the functions of plant endogenous miRNAs in Arabidopsis, soybean, and tomato [35,59,60], and CRISPR/Cas9, a powerful tool for creating soybean mutants [61], with the large datasets obtained in miRNAome sequencing studies will allow the concrete roles of specific miRNAs in the interaction between SCN and soybean to be understood. Then, the regulatory roles of miRNAs can be used respond to SCN as a powerful, highly efficient strategy to improve the yield and resistance of soybean to SCN.

Soybean, Nematode Population Culture
Glycine max, "Huipizhi" (HPZ) and Williams 82 (W82) were selected as resistant and susceptible cultivars to SCN race 3, respectively. Soybean seeds were surface sterilized with 1% NaClO for 10 min and were washed at least five times with sterilized water. Then, the seeds were geminated in PVC tubes (height × width = 10 × 3 cm) containing equal ratios of sterilized sand and soil in a climatic chamber (light/dark = 16/8, 23-26 • C, 50% relative humidity). Hoagland's nutrient solution was added to soybean seedlings once every three days. The Heterodera glycines (soybean cyst nematode, SCN) race 3 population was propagated on W82 in infested soil in a greenhouse in the Nematode Institute of Northern China (Shenyang, China). After two months, SCN cysts were extracted from the infested soil on a 60-mesh sieve (250 µm), and harvested cysts were then crushed on an 80-mesh sieve (180 µm). Eggs were collected on a 500-mesh sieve (25 µm), and were further purified using 35% (w/v) sucrose solution. Eggs were sterilized with 0.1% NaClO for 10 min before being rinsed with sterilized water several times to remove any traces of NaClO. The sterilized eggs were transferred to a modified Baermann pan with 3 mM ZnSO 4 at 25 • C in the dark for 5 days to allow them to hatch, and the freshly hatched pre-parasitic second-stage juveniles (J2s) of SCN were then harvested [62,63].

Nematode Inoculation, Penetration, Development Evaluation, and Histological Observation
Ten days after soybean seedlings geminated, a hole was dug around the soybean seedlings. One milliliter of 0.1% water-agar mixture containing approximately 2000 SCN J2s or water-agar mixture alone as a treatment or control, respectively, were added to the root systems of soybean seedlings, and these seedlings were grown under the conditions described above. Infections were synchronized by washing the infected roots at 24 h post-inoculation.
To analyze the penetration and development of SCN in the roots of HPZ and W82, 10 roots of each cultivar were harvested at 1 and 5 dpi, and were stained using the method reported in Bybd et al. [64]. Then, SCN juveniles were counted using a Nikon SMZ800 stereomicroscope (Nikon, Japan). To observe the histological changes in HPZ and W82 at 1 and 5 dpi, SCN-infected and non-infected roots of HPZ and W82 were cut into 1-cm cubes before being put into Carnoy's fix solution, dehydrated with gradient ethanol, and embedded in Technovit 7100 (Heraeus Kulzer, Hanau, Germany) according to the manufacturer's instructions. Embedded roots were sectioned (3 µm) using a Leica CM1860 UV microtome (Leica, Wetzlar, Germany) before being stained by 1% toluidine blue (Sigma-Aldrich, St. Louis, MO, USA). Pictures of the syncytia were taken using an OLYMPUS DP80 light microscope (Olympus, Tokyo, Japan), and the size of the syncytia was measured using Cellsens standard software.

Sequencing Sample Preparation
At 1 and 5 dpi, i.e., the early stages of soybean-SCN interactions for both resistant and susceptible cultivars, soybean roots were washed with running tape water and cleaned with sterilized water. Roots were then flash frozen in liquid nitrogen and stored in a −80 • C freezer. Three independent biological replicates were collected, and each replicate contained the roots of three soybean seedlings. The sampled roots were as follows: Non-inoculated W82 at 1 dpi (W1C), SCN-inoculated W82 at 1 dpi (W1N), non-inoculated HPZ at 1 dpi (H1C), and SCN-inoculated HPZ at 1 dpi (H1N). Likewise, the samples of W5C, W5N, H5C, and, H5N were taken at 5 dpi, giving a total of 24 samples.

RNA Isolation, Construction of Small RNA Libraries, and Sequencing
MicroRNA was extracted from all 24 samples with an EASYspin plant microRNA extraction kit (Aidlab biotechnologies Co., Ltd, Beijing, China), according to the manufacturer's instructions. The purity, concentration, and integrity of the microRNA samples were tested using an Agilent Bioanalyzer 2100 (Agilent, Waldbronn, Germany). MicroRNA libraries were constructed using the NEB Next Multiplex Small RNA Library Prep Kit for Illumina. Then, these microRNA libraries were sequenced on an Illumina novaseq6000 at Biomarker company (Beijing, China), and 50 single-end reads were generated. All the raw data were uploaded to SRA (https://dataview.ncbi.nlm.nih.gov/ object/PRJNA579002?reviewer=c9ctcl8rfej71udnds5o2jcn0g).

Analysis of Sequencing Data, Differential miRNA Expression Profiling, and Target Prediction
For data quality control, raw data (raw reads) in a fastq format were firstly processed through FastQC. In this step, clean reads were obtained by removing reads containing adapters, reads containing ploy-N, and low-quality reads. Further, reads were also trimmed and cleaned by removing the sequences smaller than 18 nt or longer than 30 nt. At the same time, Q30, GC-content, and sequence duplication levels of the clean data were calculated. All downstream analyses were based on high quality clean data. To comparatively analyze the 18 to 30 nt unannotated clean reads, Bowtie [65] was used to filter the clean reads using Silva database, GtRNAdb database, Rfam database, and Repbase database sequence alignment, and to filter out ribosomal RNAs (rRNA), transfer RNAs (tRNA), small nuclear RNAs (snRNA), small nucleolar RNAs (snoRNA), and other ncRNAs and repeats. The remaining reads were used to detect known miRNAs and novel miRNAs, which were predicted through comparisons with the soybean reference genome (Glycine max.Wm82.a2.v1) [32] and with known miRNAs from miRbase (v22) [33]. miRDeep2 [66] was used for novel miRNA prediction.
The expression levels of the miRNAs in each sample were calculated and normalized using the TPM (transcripts per million) algorithm [67], and the raw counts of miRNAs were used for the analysis of differential miRNA expression that was conducted by using DESeq2(1.10.1) [68]. miRNAs that were found to have |log2(FC)| ≥ 1.00; FDR (False Discovery Rate) ≤ 0.01 by DESeq2 were assigned as differentially expressed. The potential target genes of the miRNAs were predicted using TargetFinder (v1.6), and these predicted target genes were then annotated by Blast to the NR, Swiss-Prot, GO, KEGG, and Pfam databases. KOG, GO, and KEGG pathway classification was implemented on the BMK cloud analysis platform.

Quantitative Real-Time PCR (qRT-PCR) Analysis to Validate miRNA Expression
To validate the expressions of the miRNAs identified in the sequencing analysis, three miRNAs with different expression patterns in different treatments were randomly selected for qRT-PCR assays. miRNA First Strand cDNA Synthesis (Sangon Biotech, Shanghai, China) was used for polyadenylation and reverse transcription of all miRNAS, according to the manufacturer's instructions. The MicroRNAs qPCR Kit (Sangon Biotech, Shanghai) was used for miRNA qPCR with the following steps: 95 • C for 30 s, 40 cycles at 95 • C for 5 s, and 60 • C for 30 s. All qRT-PCR experiments were performed on the CFX Connect Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). Soybean U6 snRNA was used as the internal reference gene for the miRNAs, and the data were quantified using the 2 −∆∆Ct method, all the primers used for qRT-PCR were shown in Supplementary Material S5. [69].

Statistical Analysis
Nematode infection, development, and syncytium size data were checked for normality and homogeneity of variance using the Kolmogorov-Smirnov or Shapiro-Wilk test. T-tests were then used to analyze the above data using SPSS software (SPSS Inc., Chicago, IL, USA).

Conclusions
SCN is the most damaging pathogen of soybean in terms of yield losses, and control strategies mainly rely on the few resistant cultivars available. The development of new technologies to study biological questions, such as high-throughput sequencing, CRISPR/Cas9 gene editing, and other novel technologies, will provide us with new methods to understanding the mechanisms behind how soybean responds to SCN infection, which we can use to manage SCN. In this study, we constructed 24 small RNA libraries for high-throughput sequencing, in order to investigate the response of soybean miRNAs to SCN in both a susceptible cultivar (W82) and a resistant cultivar (HPZ) at the early stages of infection (1 and 5 dpi). In total, 634 known miRNAs were identified, and 252 novel miRNAs were predicted. Among these miRNAs, a total of 14 known miRNAs belonging to 13 families and 26 novel miRNAs may be involved in the response to SCN infection, as they were differentially expressed in these 24 libraries constructed from HPZ and W82. These miRNAs included legume-specific miRNAs, conserved miRNAs, and novel predicted miRNAs (for instance, gma-miR3522, gma-miR408a, and novel_miR_106). Further analysis of the potential target genes of these DE miRNAs found that the target genes were highly related to plant defense, and that signal transduction, immune system processes, and nutrient metabolism were responsive to SCN invasion in KOG, GO, and KEGG pathway analysis. In summary, our findings provide new knowledge of the mechanisms behind soybean-SCN interaction at the miRNA level, and may act as a stepping stone for controlling the damage caused by SCN.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

SCN
soybean cyst nematode miRNA microRNA qRT-PCR quantitative real-time PCR