Identification of Long Intergenic Noncoding RNAs in Rhizoctonia cerealis following Inoculation of Wheat

ABSTRACT Wheat sharp eyespot caused by Rhizoctonia cerealis is primarily a severe threat to worldwide wheat production. Currently, there are no resistant wheat cultivars, and the use of fungicides is the primary method for controlling this disease. Elucidating the mechanisms of R. cerealis pathogenicity can accelerate the pace of the control of this disease. Long intergenic noncoding RNAs (lincRNAs) that function in plant-pathogen interactions might provide a new perspective. We systematically analyzed lincRNAs and identified a total of 1,319 lincRNAs in R. cerealis. We found that lincRNAs are involved in various biological processes, as shown by differential expression analysis and weighted correlation network analysis (WGCNA). Next, one of nine hub lincRNAs in the blue module that was related to infection and growth processes, MSTRG.4380.1, was verified to reduce R. cerealis virulence on wheat by a host-induced gene silencing (HIGS) assay. Following that, RNA sequencing (RNA-Seq) analysis revealed that the significantly downregulated genes in the MSTRG.4380.1 knockdown lines were associated mainly with infection-related processes, including hydrolase, transmembrane transporter, and energy metabolism activities. Additionally, 23 novel microRNAs (miRNAs) were discovered during small RNA (sRNA) sequencing (sRNA-Seq) analysis of MSTRG.4380.1 knockdown, and target prediction of miRNAs suggested that MSTRG.4380.1 does not act as a competitive endogenous RNA (ceRNA). This study performed the first genome-wide identification of R. cerealis lincRNAs and miRNAs. It confirmed the involvement of a lincRNA in the infection process, providing new insights into the mechanism of R. cerealis infection and offering a new approach for protecting wheat from R. cerealis. IMPORTANCE Rhizoctonia cerealis, the primary causal agent of wheat sharp eyespot, has caused significant losses in worldwide wheat production. Since no resistant wheat cultivars exist, chemical control is the primary method. However, this approach is environmentally unfriendly and costly. RNA interference (RNAi)-mediated pathogenicity gene silencing has been proven to reduce the growth of Rhizoctonia and provides a new perspective for disease control. Recent studies have shown that lincRNAs are involved in various biological processes across species, such as biotic and abiotic stresses. Therefore, verifying the function of lincRNAs in R. cerealis is beneficial for understanding the infection mechanism. In this study, we reveal that lincRNAs could contribute to the virulence of R. cerealis, which provides new insights into controlling this pathogen.

threat to wheat production (5). For example, in China, wheat sharp eyespot caused more than 15 million dollars in economic damage from 2005 to 2008 (4), and in the decade since 2005, this disease affected 9 million hectares of wheat (6). In addition, wheat sharp eyespot disease has resulted in field losses in other countries such as Chile (7).
A reliable approach for controlling this disease is breeding resistant varieties (4). However, there are still no viable resistant cultivars for this disease (5, 6) despite many efforts to study resistance to R. cerealis at the wheat level. The use of fungicides is still the primary method for controlling this disease (5,8), negatively affecting the environment and increasing the cost to farmers. Revealing the molecular basis of R. cerealis development and infection may provide new insights into the control of this disease. Long noncoding RNAs (lncRNAs), noncoding RNAs (ncRNAs) with long transcripts (.200 nucleotides [nt]) and limited protein-coding potential (9)(10)(11), have served as essential regulators of transcription in various biological processes, including growth and development responses to abiotic and biotic stresses, with many studies on animals and plants (12,13). Based on genome location, lncRNAs are divided into three categories, long intergenic noncoding RNAs (lincRNAs), intronic noncoding RNAs (incRNAs), and long noncoding natural antisense transcripts (lncNATs) (14). Several recent studies have demonstrated that lncRNAs play essential roles in plant-pathogen interactions, most of which are lncRNAs in the host (15)(16)(17)(18)(19)(20). In contrast, there are just a few studies on lncRNAs in pathogens. An lncNAT, UvlncNATMFS, has been shown to regulate growth and stress responses in the phytopathogen Ustilaginoidea virens (13). Functional verification of lncRNA009491 and lncRNA012077 suggests that lncRNAs regulate virulence in Verticillium dahlia through positive and negative regulation (21). The roles of lncRNAs in phytopathogens could be used to alleviate plant diseases caused by pathogens. However, few lncRNAs have been functionally validated in fungi, mainly in yeasts (22,23). Therefore, exploring the roles of lncRNAs in R. cerealis may provide a new perspective to control the disease caused by R. cerealis and can also provide some information for better studying their roles in pathogenic filamentous fungi.
At present, there has been no report on lncRNAs in R. cerealis, and it is unknown whether lncRNAs play a role in R. cerealis development and infection. To address this question, we report the genome-wide identification of R. cerealis lincRNAs to identify lincRNAs potentially essential for development and infection processes. We predicted the functions of genes coexpressed with lincRNAs and identified a lincRNA that could affect the virulence of R. cerealis. Finally, we revealed that the lincRNA could regulate genes involved infection-related processes, not as a competitive endogenous RNA (ceRNA). Our findings could help elucidate how lincRNAs regulate R. cerealis infection of wheat and provide clues for disease control.

RESULTS
Identification, characterization, and expression profiles of lincRNAs in R. cerealis. To identify lincRNAs confidently in the R. cerealis genome, we performed the following transcript-filtering steps (Fig. 1A). First, transcripts of .200 nt long with class code "u" were chosen after annotating the assembled transcripts using Gffcompare, and a total of 2,680 transcripts were retained. Next, we kept 1,372 transcripts categorized as "noncoding" based on CPC2 (Coding Potential Calculator 2) and CNCI (Coding-Noncoding Index) and with no Pfam hit by Pfamscan (see Fig. S1 in the supplemental material). Subsequently, we removed 53 low-confidence lincRNAs (fragments per kilobase per million [FPKM] value of ,1 in all RNA sequencing [RNA-Seq] libraries). Finally, a total of 1,319 transcripts were identified as lincRNAs in R. cerealis (Data Set S1).
To characterize the genomic features of the lincRNAs in R. cerealis, we analyzed the exon distributions, lengths of the transcripts, gene expression levels, and GC contents of these lincRNAs and mRNAs. lincRNAs had an average of 3.7 exons, significantly fewer than the 6.4 exons found in mRNAs. Only 41.3% of the mRNAs had fewer than 4 exons, whereas 77.3% of the lincRNAs had fewer than 4 exons (Fig. 1B). Furthermore, twoexon lincRNA transcripts were the most abundant in this study. The mean and median lengths of lincRNAs were 1,132 bp and 783 bp, respectively, shorter than those of mRNAs, 1,490 bp and 1,238 bp, respectively. Overall, the lincRNA size distribution ranged from 200 to 7,842 bp, with 62% ranging from 201 to 1,000 bp, while most mRNAs (about 61.3%) were more than 1,000 bp long (Fig. 1C). In terms of FPKM values, the expression level of lincRNAs (mean, 2.8) was notably lower than that of mRNAs (mean, 15.89) (Fig. 1D). In addition, the GC content in lincRNAs (48.6%) was significantly lower than that in mRNAs (52%) (Fig. 1E). To determine whether R. cerealis lincRNAs could contribute to the infection process, we examined the expression of lincRNAs after inoculation (treatment) or noninoculation (control) of wheat at 7, 14, and 21 days postinoculation (dpi), respectively. We discovered that 654 lincRNAs were differentially expressed in the treatment group compared to the control group (false discovery rate [FDR] of ,0.05 and absolute log 2 fold change [jlog 2 fold changej] of $1). Of these 654 differentially expressed lincRNAs (DELs), 294, 312, and 442 lincRNAs were differentially expressed at 7, 14, and 21 dpi, respectively (Fig. 1F). Most of these DELs were downregulated, and the number of differentially expressed lincRNAs increased with increasing durations of infection. Totals of 96, 84, and 188 DELs were unique and differentially expressed at 7 dpi, 14 dpi, and 21 dpi, respectively (Fig. 1G). Especially at 21 dpi, about 42.5% of the DELs (442) were unique. One hundred eight of these DELs were shared by all three infection stages. These findings indicated that R. cerealis lincRNAs might play essential roles in infection. Functional analysis of infection and development associated with lincRNAs. To predict the putative function of lincRNAs in the infection process of R. cerealis, a coexpression network was constructed to associate 654 DELs and 10,520 differentially expressed genes (DEGs) derived from a previous study by Zeng et al. (24). A total of 11 modules (except for the gray module), which captured 624 DELs and 10,439 DEGs, were identified (Table 1). To verify whether lincRNAs are involved in the development-and infectionrelated processes of R. cerealis, the Gene Ontology (GO) and KEGG pathway databases were used to analyze the genes coexpressed with lincRNAs in each module (Data Set S2). Among these 11 modules, genes in the blue module were significantly enriched in the steroid biosynthesis (GO:0016126), biological process and steroid biosynthesis (Ko00100), mitogen-activated protein kinase (MAPK) signaling (Ko04011), and pathogenic Escherichia coli infection (Ko05130) pathways ( Fig. 2B and C). These findings suggested that the lincRNAs in the blue module likely regulated R. cerealis growth and infection. To further screen potential key lincRNAs in this module, we finally identified nine hub lincRNAs in the blue module with the criteria of an intramodular gene connectivity value of .0.4 and a degree of .50 (Fig. 2D).
Confirmation of the expression patterns of hub lincRNAs. We used quantitative real-time PCR (qRT-PCR) to examine the temporal expression patterns of these nine lincRNAs at 0, 2, 4, 6, 8, 10, 12, and 14 dpi (Fig. 3). Compared to the RcActin gene, linc2, linc6, linc7, and linc9 were significantly upregulated. At the same time, linc1, linc3, and linc5 were significantly downregulated. Generally, disease-related genes are highly expressed in the early stage of infection. Additionally, at 2 dpi, the linc6, linc7, and linc9 expression levels were 16, 23, and 16 times higher than those of RcActin, suggesting that these three lincRNAs may be essential for the process of R. cerealis infection. Therefore, these three lincRNAs were selected for further analysis.
MSTRG.4380.1 conferred pathogenicity to wheat. Due to the severe lag in functional genomics studies, there is no transformation system for R. cerealis. To further demonstrate whether these three hub lincRNAs were involved in the infection process, we knocked down their expression using the barley stripe mosaic virus host-induced gene silencing (BSMV-HIGS) system, which is widely used to verify the function of pathogenicity genes in fungi. At 14 dpi, the BSMV-infected wheat showed mild chlorotic mosaic symptoms (Fig. 4A), and the endogenous phytoene desaturase (PDS) gene was used as the positive control for BSMV infection symptoms. Subsequently, the leaves showing mild chlorotic mosaic symptoms were selected for R. cerealis infection. According to statistical analysis, only the lesion areas on the leaves of the inoculated linc9 HIGS plants showed macroscopic reductions compared to the negative controls (Fig. 4B). In contrast, the lesion areas of the plants with silenced linc6 and linc7 did not show a significant difference compared to the control. The biomass of R. cerealis decreased significantly in linc9-silenced plants (Fig. 4C). Furthermore, qRT-PCR analysis revealed that endogenous lincRNA was  (Fig. 4D). These results suggest that the silencing of linc9 could dramatically reduce the pathogenicity of R. cerealis, and this was selected for intensive study. Moreover, to elucidate the specific process that MSTRG.4380.1 participates in, we retrieved its coexpressed genes with the criterion of an intramodular gene connectivity value of .0.3. GO analysis indicated that these genes were enriched in categories related to membrane (GO:0016020), RNA binding (GO:0003723), calcium ion binding (GO:0005509), chitin synthase activity (GO:0004100), transferase activity, transferring hexosyl groups (GO:0016758), and lipid biosynthetic process (GO:0008610). KEGG pathway analysis showed that these genes were involved in the ribosome biogenesis (Ko03008), steroid biosynthesis (Ko00100), and sphingolipid metabolism (Ko00600) pathways ( Fig. 4E and F). Because lncRNA preferentially regulates the expression of nearby genes, we chose coding genes 10 kb upstream and 20 kb downstream of MSTRG.4380.1 as its cis target genes (Data Set S3A). Next, we intersected these genes with their coexpressed genes (Data Set S3B), but no candidate genes were retained, indicating that MSTRG.4380.1 regulates genes in trans.
RNA-Seq analysis of linc9 knockdown. Moreover, to identify genes regulated by linc9, we conducted RNA-Seq analysis with infected wheat leaves collected at 5 dpi lincRNAs Involved in Rhizoctonia cerealis Infection Microbiology Spectrum when the expression of linc9 was significantly knocked down compared with the control (Fig. 4D). A total of 953 differently expressed genes were identified in R. cerealis (jlog 2 fold changej of $1 and adjusted P [P adj ] value of #0.05), including 111 upregulated and 842 downregulated genes (Fig. 5A). Compared with the wheat inoculated with BSMV:: g , the expression of linc9 was significantly reduced in BSMV::linc9 plants (Fig. 5B). These results further demonstrated that HIGS could knock down the expression of linc9 and also indicated that the knockdown of linc9 could affect gene expression in R. cerealis. To determine the processes that linc9 is involved in, we implemented GO and KEGG pathway enrichment analyses with 842 downregulated and 111 upregulated genes. Interestingly, GO enrichment analysis of 842 downregulated genes showed that these genes were located mainly in the membrane and extracellular region and were associated with infection-related processes, including hydrolase activity and transmembrane transporter activity (Fig. 5C). KEGG enrichment analysis showed that these genes were associated with energy metabolism, such as starch and sucrose metabolism and amino sugar and nucleotide sugar metabolism (Fig. 5D). However, GO enrichment analysis showed that these upregulated genes were enriched in the glycolytic process and proton transmembrane transporter categories, and no KEGG item was enriched for these genes (Fig. S2). These results indicate that linc9 may positively regulate R. cerealis infection by regulating genes involved in hydrolase, transmembrane transporter, and energy metabolism activities. Identification and target prediction of miRNAs. lncRNAs can act as ceRNAs competing with target transcripts for microRNAs (miRNAs), which interferes with the function of miRNAs (25). To verify whether linc9 functions as a ceRNA, we conducted a small RNA (sRNA) sequencing (sRNA-Seq) analysis with infected wheat leaves collected at 5 dpi when the expression of linc9 was significantly knocked down compared with the control. After quality control, a total of 227,414,158 clean reads were retained, for which the length distribution of reads in each sample is shown in Fig. S3. After filtering rRNA, tRNA, snRNA, and snoRNA, the retained reads were used to identify novel miRNAs by miRDeep2, resulting in a total of 23 miRNAs being identified (Data Set S4). Next, we verified whether these miRNAs could target linc9 by using psRobot (26) and found that no miRNAs identified here targeted linc9, which indicates that linc9 does not regulate R. cerealis infection as a ceRNA.

DISCUSSION
R. cerealis is a devastating pathogen threatening worldwide wheat production (5,8). Although some progress has been made on the mechanism of resistance of wheat to R. cerealis (6), there are no wheat cultivars resistant to R. cerealis at the moment. Currently, the use of fungicides is the primary measure for the control of this disease, although fungicide use can cause environmental pollution and increase labor costs. Revealing the mechanism of R. cerealis infection and development may provide clues for disease control. lncRNAs play critical roles in diverse biological processes in fungi,

and lncRNAs may help gain new insights into the pathogen infection mechanism.
Here, we performed genome-wide identification and functional prediction of R. cerealis lincRNAs during infection stages and discovered a lincRNA, MSTRG.4380.1, that reduced pathogenicity after its expression was silenced by the regulation of genes associated with infection-related processes.
lncRNAs have been shown to play critical roles in diverse biological processes in a wide range of species. However, as for plant-pathogenic filamentous fungi, there are only a few studies on lncRNAs, including those of Ustilago maydis (27), Fusarium graminearum (28), Ustilaginoidea virens (13), Verticillium dahlia (21), and Magnaporthe oryzae (29,30). Naturally, the roles that lncRNAs can play in filamentous pathogenic fungi are less well known. Here, we performed genome-wide identification of R. cerealis lincRNAs, identifying a total of 1,319 lincRNAs. In addition, we found that R. cerealis lincRNAs may be involved extensively in multiple biological processes, for which half of the lincRNAs were differentially expressed during the different infection stages. These findings will provide the basis for further research on the role of R. cerealis lincRNAs and a reference for research on other pathogens.
It is challenging to confirm the role of lncRNA in RNA-mediated gene regulation (31). Although thousands of lncRNAs have been identified in phytopathogens (13,21,29,(32)(33)(34), only a few lncRNAs have been functionally characterized (13,21,33). Besides, functional genomics studies have lagged seriously, resulting from the lack of a genetic transformation system, which has hindered the functional verification of R. cerealis genes. Since the discovery of transkingdom RNA silencing, however, HIGS has been exploited as a powerful genetic tool to validate the function of pathogenicity genes during infections and has been applied to a variety of pathogenic fungi, including Blumeria graminis (35), Puccinia triticina (36), Fusarium graminearum (37), Puccinia striiformis f. sp. tritici (38), Verticillium dahlia (39), and Phakopsora pachyrhizi (40). Our latest work confirmed that HIGS can be applied to R. cerealis as a functional characterization tool (41). Here, A lincRNA, MSTRG.4380.1, was verified to affect the virulence of R. cerealis by HIGS, which reduced pathogenicity after its expression was silenced. Sprayinduced gene silencing (SIGS), which was also confirmed to be a practical tool for  (41), has recently served as an innovative method of crop protection (42)(43)(44). By combining HIGS and SIGS technologies, this discovery may provide new insight into the control of wheat sharp eyespot. Generally, lncRNAs may regulate the expression of neighboring genes in cis and distant genes in trans (45,46), while most of the lncNATs regulate neighboring genes in cis. For example, UvlncNATMFS might affect the growth of Ustilaginoidea virens by regulating its associated sense gene UvMFS (13). In contrast to lncNATs, the expression of lincRNAs is not significantly correlated with that of their neighboring genes (47), which is further supported by our results. The coexpressed coding genes of MSTRG.4380.1 were enriched in the membrane and involved in chitin synthase activity and steroid biosynthesis, indicating that these potential targets might contribute to infection-related processes and membrane formation. After knocking down MSTRG.4380.1, the downregulated genes were enriched mainly in hydrolase activity, transmembrane transporter activity, and energy metabolism activities, which further supports that MSTRG.4380.1 regulates genes involved in infectionrelated processes. Additionally, lncRNAs can also serve as ceRNAs to interfere with miRNA function by competing with target genes for miRNAs (25). However, in this study, our results suggested that for MSTRG.4380.1, this was not the case. Nonetheless, the 23 novel miRNAs identified in R. cerealis may provide a basis for the analysis of transkingdom miRNAs in the R. cerealis-wheat interaction.
In conclusion, this first identification of lincRNAs and miRNAs in R. cerealis will provide a basis for studying the roles of lincRNAs and miRNAs in R. cerealis-wheat interactions. It was verified that the lincRNA MSTRG.4380.1 could reduce the virulence of R. cerealis by regulating genes involved in infection-related processes. This may provide new insight into the control of wheat sharp eyespot by combining HIGS and SIGS technologies. However, much effort is still needed to understand the specific regulatory mechanism of MSTRG.4380.1.

MATERIALS AND METHODS
Identification of lincRNAs. Our most recent genome assembly and annotation of R. cerealis were used (24). Published dual-time-course RNA-Seq data were reanalyzed to identify lincRNAs involved in R. cerealis infection and growth (24). These data comprised 18 samples. Wheat leaves inoculated with R. cerealis were harvested at 7, 14, and 21 dpi, with three biological replicates for the treatment. Hyphae grown on potato dextrose agar (PDA) medium were collected simultaneously as the control. All raw data containing adaptors and low-quality reads were removed by using Fastp v0.20.1 (48). Next, the clean reads were aligned to the reference genome Rce_V1 using HISAT2 v2.2.1 (49). The alignments were assembled using StringTie v2. 1.4 (50) with the available R. cerealis annotation data as a reference. Next, the assembled transcripts were merged into a uniform set of transcripts for all samples by using StringTie. Subsequently, Gffcompare v0.10.1 (51) was used to annotate the assembled transcripts, and transcripts with class code "u" were selected. Furthermore, transcripts with lengths of ,200 nt were excluded. The coding potentials for the remaining transcripts were then predicted by using CPC2 (Coding Potential Calculator 2) (52), CNCI (Coding-Noncoding Index) (53), and Pfamscan (E value of ,1e25). Transcripts with no coding potential shared by the three software were retained. Finally, the remaining transcripts with a particular expression level (FPKM of .1 in at least one RNA-Seq library) were retained as candidate lincRNAs.
Analysis of differential expression patterns. The count matrix of lincRNAs was obtained using the prepDE.py program in StringTie. Differentially expressed lincRNAs were screened with the criteria of an adjusted P value of ,0.05 and a jlog 2 fold changej of $1 using the DESeq2 R package (54).
Weighted correlation network analysis. The expression data for DEGs derived from a previous study by Zeng et al. (24) and DELs were extracted to construct a coexpression network using the WGCNA (weighted correlation network analysis) R package (55) with a soft-threshold power of 13 and a minimum module size of 30. Other parameters were set to default values.
Functional enrichment analysis. The GO enrichment and KEGG pathway analyses were performed by using the clusterProfiler R package (56). The GO terms and KEGG pathways with P values of ,0.05 were significantly enriched.
Identification of hub lincRNAs. For the blue module, we calculated the connectivity of each gene based on its intramodular connectivity. lincRNAs with an intramodular connectivity value of .0.4 and a degree of .50 were defined as hub lincRNAs. Cytoscape v3.8.2 was used for network visualization (57).
RNA preparation and expression analysis. To define the temporal expression patterns of the lincRNAs during wheat-R. cerealis interactions, inoculated leaves were harvested at 2, 4, 6, 8, 10, 12, and 14 dpi. As a control, the hyphae of R. cerealis in vitro (PDA) were collected at the time of inoculation (0 dpi). Specific qRT-PCR primers (see Table S1 in the supplemental material) were designed to distinguish each lincRNA.
A spin column fungal total RNA purification kit (Sangon Biotech, China) was used to extract RNA.
lincRNAs Involved in Rhizoctonia cerealis Infection Microbiology Spectrum qRT-PCR was performed on the CFX Connect real-time system (Bio-Rad, USA) according to a procedure described previously by Yang et al. (58). The relative expression levels of target lincRNAs in R. cerealis were determined using the 2 2DDCT method (59); the expression levels of all of these lincRNAs were normalized to the expression level of RcActin. Each sample was analyzed in three biological replicates. Potential target gene prediction. Based on the genome location of MSTRG.4380.1 relative to the neighboring genes, the protein-coding genes transcribed within a 10-kb window upstream and a 20-kb window downstream were screened as potential cis-regulated target genes with the window function of BEDTools software (60).
BSMV-HIGS assays. The wheat variety Fielder was selected for BSMV-induced gene silencing in the target lincRNAs. A total of 30 plots of wheat seedlings were used in the experiment, with 12 seedlings in each plot. Among these 30 plots of wheat seedlings, 6, 6, 6, 6, 4, and 2 plots were used for linc6, linc7, lic9, g , mock, and TaPDS, respectively.
The capping and in vitro transcription of seven BSMV plasmids (BSMV::a, BSMV::b, BSMV:: g , BSMV:: linc6, BSMV::linc7, BSMV::linc9, and BSMV::TaPDS) were performed by using the Ribom7G cap analog and the RiboMAX large-scale RNA production system T7 (Promega, USA). The linearized in vitro transcripts of five fragments (BSMV:: g , BSMV::linc6, BSMV::linc7, BSMV::linc9, and BSMV::TaPDS) were individually mixed with BSMV:a and BSMV::b and then inoculated onto the second expanded leaves of the wheat seedlings (58). These virus-infected plants were maintained in a temperature-controlled chamber at 25°C. After 14 days, the fourth leaves of wheat with mild chlorotic mosaic symptoms were inoculated with R. cerealis isolate R0301 mycelium plugs and maintained at 23°C as described previously by Li et al. (41). The infected fourth leaves were harvested at 3 and 5 dpi for the qRT-PCR study. The fungal biomass was determined by qRT-PCR, as previously described by Huai et al. (61). After 7 days of inoculation, the phenotypes of the inoculated leaves were photographed.
RNA and sRNA library construction and sequencing. The infected fourth leaves treated with BSMV: linc9 and BSMV: g were harvested at 5 dpi for RNA and sRNA sequencing (three repeats for each treatment). The construction and sequencing of RNA and unique molecular identifiers small RNA were carried out by BGI (http://www.bgitechsolutions.com/).
RNA-Seq analysis. After sequencing, the clean reads were retained after the removal of adaptor sequences, contamination, and low-quality reads from raw reads by using SOAPnuke (62). Next, the quality of the data was examined by using fastQC. The clean reads were aligned to the R. cerealis reference genome using HISAT2 (49), and read counts and FPKM were determined by using featureCounts for gene expression. Significant differential expression was defined as a jlog2 fold changej of $1 and a P adj value of #0.05 using DESeq2. GO categories matched with significantly up-and downregulated genes were identified as described above.
Identification of miRNAs and target prediction of miRNAs. After sequencing, the raw reads were filtered by using SOAPnuke (62). Data filtering includes the removal of adaptor sequences, contamination, and low-quality reads from raw reads. Next, clean reads were annotated by using the cmscan program of the Infernal package (63,64), and reads with high similarities with rRNAs, tRNAs, snRNAs, and snoRNAs were removed. miRDeep2 was used to identify novel miRNAs with no related species in miRbase and quantify miRNAs (65). The potential targets of miRNAs were predicted by using psRobot with default parameters (26).