Non-coding RNAs identification and regulatory networks in pathogen-host interaction in the microsporidia congenital infection
BMC Genomics volume 24, Article number: 420 (2023)
The interaction networks between coding and non-coding RNAs (ncRNAs) including long non-coding RNA (lncRNA), covalently closed circular RNA (circRNA) and miRNA are significant to elucidate molecular processes of biological activities and interactions between host and pathogen. Congenital infection caused by vertical transmission of microsporidia N. bombycis can result in severe economic losses in the silkworm-feeding industry. However, little is known about ncRNAs that take place in the microsporidia congenital infection. Here we conducted whole-transcriptome RNA-Seq analyses to identify ncRNAs and regulatory networks for both N. bombycis and host including silkworm embryos and larvae during the microsporidia congenital infection.
A total of 4,171 mRNAs, 403 lncRNA, 62 circRNAs, and 284 miRNAs encoded by N. bombycis were identified, among which some differentially expressed genes formed cross-talk and are involved in N. bombycis proliferation and infection. For instance, a lncRNA/circRNA competing endogenous RNA (ceRNA) network including 18 lncRNAs, one circRNA, and 20 miRNAs was constructed to describe 14 key parasites genes regulation, such as polar tube protein 3 (PTP3), ricin-B-lectin, spore wall protein 4 (SWP4), and heat shock protein 90 (HSP90). Regarding host silkworm upon N. bombycis congenital infection, a total of 14,889 mRNAs, 3,038 lncRNAs, 19,039 circRNAs, and 3,413 miRNAs were predicted based on silkworm genome with many differentially expressed coding and non-coding genes during distinct developmental stages. Different species of RNAs form interacting network to modulate silkworm biological processes, such as growth, metamorphosis and immune responses. Furthermore, a lncRNA/circRNA ceRNA network consisting of 140 lncRNAs, five circRNA, and seven miRNAs are constructed hypothetically to describe eight key host genes regulation, such as Toll-6, Serpin-6, inducible nitric oxide synthase (iNOS) and Caspase-8. Notably, cross-species analyses indicate that parasite and host miRNAs play a vital role in pathogen-host interaction in the microsporidia congenital infection.
This is the first comprehensive pan-transcriptome study inclusive of both N. bombycis and its host silkworm with a specific focus on the microsporidia congenital infection, and show that ncRNA-mediated regulation plays a vital role in the microsporidia congenital infection, which provides a new insight into understanding the basic biology of microsporidia and pathogen-host interaction.
Microsporidia are a group of fungi-related obligate intracellular parasites that infect most animal groups including humans [1, 2]. As the first reported microsporidia, Nosema bombycis was described by Nageli in 1857. As of now, more than 1,500 species belonging to 197 genera have been reported , and 700 species are parasitic on insects alone . Microsporidia are capable of vertical transmission, N. bombycis transovarial transmission can cause a highly fatal silkworm disease (pébrine). Therefore, N. bombycis is the only mandatory inspection target for silkworm egg production by silkworm feeding countries.
It has been known that ncRNAs, including miRNAs, lncRNAs, and circRNAs play an important role in many biological processes via regulating gene expression at both transcriptional and post-transcriptional levels [5,6,7]. MiRNAs, usually referred to as small ncRNAs ranging in size from 19 to 22 nucleotides, are well conserved in eukaryotic organisms. They are thought to be a class of critical post-transcriptional regulators to play a vital role in host–pathogen interaction. In the nucleus, primary miRNA transcripts (pri-miRNAs) are processed by the microprocessor complex, consisting of an RNA binding protein DiGeorge Syndrome Critical Region (DGCR8) and a ribonuclease III enzyme, Drosha, and generate precursor miRNAs (pre-miRNAs) . The pre-miRNAs are subsequently cleaved by Dicer to produce an ∼22-nt miRNA/miRNA* duplex . Although both strands of the duplex are necessarily produced in equal amounts by transcription, their accumulation is asymmetric at the steady state. By convention, the most abundant duplex strand is defined as the mature miRNA strand, whereas the less abundant strand is known as the miRNA* strand . For example, Mycobacteria tuberculosis H37Rv induced an increase in miR-33 and miR-33* expression in peritoneal macrophages to inhibit autophagy and enhance host lipid metabolism to facilitate their survival . In addition, B. bassiana deploys bba-milR1 to modulate mosquito immunity by targeting the mosquito genes Spätzle4 (Spz4) and CLIP-domain serine protease (CLIPB9) . RNA interference (RNAi) and related RNA silencing pathways are widespread in animals, plants, fungi, and protozoans, and control gene silencing in all living cells . RNAi relies on small noncoding RNAs, such as miRNAs and siRNAs. Two proteins, Dicer and Argonaute (Ago), are responsible for small noncoding RNA biogenesis. Among microsporidia, Huang et al. showed that 11 species including N. apis, N. ceranae and N. bombycis maintain both Dicer and Ago orthologs , suggesting that those microsporidia retained the canonical RNAi signaling pathway. A few studies have revealed that miRNAs played an essential role in N. bombycis-silkworm interactions [15, 16]. However, the identification of miRNAs in N. bombycis-silkworm interactions, especially in the congenital infection of N. bombycis in silkworm including embryo and larvae stages, still remain elusive.
LncRNAs are a novel class of ncRNAs with nucleotides from 200 bp to several kilobases in size. Structurally, lncRNAs are similar to mRNAs, and are subject to post-transcriptional modifications such as capping, polyadenylation, and splicing. Unlike miRNAs, lncRNAs are generally less conserved in eukaryotic organisms . LncRNAs can compete with endogenous RNAs, acting as miRNA sponges to reduce their regulatory effect on target mRNAs involved in various of biological processes [18, 19]. In microsporidia, Guo et al. have reported that 83 novel lncRNAs including 59 antisense lncRNAs, 21 long intergenic noncoding RNAs, and three sense lncRNAs produced by N. ceranae . In addition, Wu et al. have described that 1,180 lncRNAs were identified in the silkworm, including 6,250 intergenic lncRNAs, 474 intronic lncRNAs and 5,086 natural antisense lncRNAs . Chen et al. have demonstrated that lncRNAs and their regulatory networks were involved in the host honeybee response to N. ceranae infection . Furthermore, one lncRNA lncRNA4.2 can act as a decoy and titrate away dimerization of Toll, resulting in dysregulation of the Toll signaling pathway in silkworm after virus infection . Up to date, the function of N. bombycis and silkworm lncRNAs in host–pathogen interactions, especially in the congenital infection of N. bombycis in silkworm is still obscure.
CircRNAs are another novel class of ncRNAs with a special structure, which is characterized with a covalent closed-loop without 5’ caps or 3’ poly (A) tails. Those circRNAs are produced from back-splicing of precursor mRNA (pre-mRNA) in eukaryotes, and are generally expressed at low levels . Many circRNAs are more abundant than their linear counterparts , and their expression levels are generally tissue-specific and cell-type-specific . There is increasing evidence that circRNAs function as sponges for miRNAs [27, 28], regulating transcription, splicing , or polypeptide production [30, 31]. In microsporidia, 204 circRNAs produced by N. ceranae have been identified by using deep sequencing technology, of which include 174 exonic circRNAs and 30 intergenic circRNAs . Furthermore, ceRNA regulatory networks were constructed, and 15 N. ceranae circRNAs were found to act as sponges of the corresponding three miRNAs (Nce-miR9769, Nce-miR20213, and Nce-miR34537). However, current research on microsporidia and host circRNAs mainly focuses on circRNAs identification, and the studies on the function of circRNAs in host-microsporidia interactions are scant.
To obtain a better understanding of various species of ncRNAs in N. bombycis and its host B. mori, we have performed small RNA and strand-specific library construction for high-throughput RNA sequencing and genome-wide analysis of the coding and noncoding RNAs in the N. bombycis congenitally infected and uninfected silkworm embryos of 5-day, and larvae of 1, 5, 10-day, respectively. The expression and cross-talk of mRNA, miRNA, lncRNA, and circRNA were analyzed. Based on ceRNA hypothesis, the lncRNA/circRNA-miRNA- mRNA networks were constructed to understand the functional effects of lncRNAs/circRNAs, through sponging miRNAs. Furthermore, we predicted that parasite and host genes were potentially targeted by parasite and/or silkworm miRNAs. This is the first pan-transcriptome analysis of miRNA, lncRNA, circRNA and mRNA expression profiles and their cross-talks in the embryos and larvae of insects following the congenital microsporidia infection, which provides a new insight into understanding the basic biology of microsporidia and pathogen-host interaction.
Materials and methods
Silkworm embryos and larva congenitally infected by N. bombycis
The silkworm 305 strain was reared at 26℃ with a 70% relative humidity and 12-h light/dark cycle. The N. bombycis isolate CQ1 (Chongqing, China) was conserved in the China Veterinary Culture Collection Center (CVCC No. 102059). The silkworms were orally challenged by mature spores with approximately 106 spores per larva at the fifth instar larvae and were then reared under the above feeding conditions until cocoon. Surviving pupae hatched into adults and fertilized female adult laid progeny eggs. Lastly, the progeny eggs were collected and hatched, and hatched larvae were reared under the above-mentioned feeding conditions until the 10th day of larvae.
Embryos and larvae collection
To confirm the congenital N. bombycis infection in silkworm, eggs produced by mothers infected with N. bombycis were hatched. The hatched eggs of 5-day and larvae of 1, 5, 10-day were collected to detect the presence of spores by hematoxylin & eosin (H&E) staining and transmission electron microscopy (TEM). In addition, to sequence the whole transcriptome, a 0.6 g hatched eggs of 5-day and larvae of 1, 5, 10-day from either the N. bombycis congenitally infected or uninfected control groups were collected.
The collected eggs were immersed in 10% potassium hydroxide solution (w/v) at room temperature for 7 min to soften cell walls followed by fixing with Smith fixative (0.5 g potassium dichromate, 2.5 mL glacial acetic acid, 10 mL methanol, 87.5 mL distilled water) for 12 h. The silkworm tissues were dehydrated in a series of ethanol and then embedded in paraffin and cross-sectioned into 5–7 μm slices. After deparaffination, the sections were treated with xylene and washed with distilled water. Lastly, the sections were stained with H&E straining according to the manufacturer’s protocol, and the sections were observed under a light microscope.
Transmission electron microscopy
The collected silkworm embryos and larvae were fixed with the 2.5% glutaraldehyde 0.1 M phosphate buffer saline (PBS) solution. The sections were cut using an ultramicrotome (Lecia, Wetzlar, Germany) and placed on nickel grids. The sections were stained in 3% uranyl acetate, followed by lead citrate. Lastly, the sections were washed three times in ddH2O, and photographed with a JEM-1400 Plus TEM at an accelerating voltage of 80 kV.
Non-infected (normal) eggs were placed in an incubator at -26℃. On the fifth day, a certain number of eggs were taken out and divided into three samples, designated NI-E5-1 (non-infected embryos day 5 replicate 1), NI-E5-2 (non-infected embryos day 5 replicate 2), and NI-E5-3 (non-infected embryos day 5 replicate 3). The remaining eggs were placed and hatched into larvae in an incubator. On the first day of larvae, a certain number of larvae were taken out and divided into three samples, designated NI-L1-1 (non-infected larvae day 1 replicate 1), NI-L1-2 (non-infected larvae day 1 replicate 2), and NI-L1-3 (non-infected larvae day 1 replicate 3). On the fifth day of larvae, a certain number of larvae were taken out and divided into three samples, designated NI-L5-1 (non-infected larvae day 5 replicate 1), NI-L5-2 (non-infected larvae day 5 replicate 2), and NI-L5-3 (non-infected larvae day 5 replicate 3). On the tenth day of larvae, a certain number of larvae were taken out and divided into three samples, designated NI-L10-1 (non-infected larvae day 10 replicate 1), NI-L10-2 (non-infected larvae day 10 replicate 2), and NI-L10-3 (non-infected larvae day 10 replicate 3). Similarly, N. bombycis-infected eggs were also placed in an incubator at -26℃。Infected-samples were taken at the same way as non-infected samples, named I-E5-1 (infected embryos day 5 replicate 1), I-E5-2 (infected embryos day 5 replicate 2), I-E5-3 (infected embryos day 5 replicate 3), I-L1-1 (infected larvae day 1 replicate 1), I-L1-2 (infected larvae day 1 replicate 2), I-L1-3 (infected larvae day 1 replicate 3), I-L5-1 (infected larvae day 5 replicate 1), I-L5-2 (infected larvae day 5 replicate 2), I-L5-3 (infected larvae day5replicate 3), I-L10-1 (infected larvae day 10 replicate 1), I-L10-2 (infected larvae day 10 replicate 2), and I-L10-3 (infected larvae day 10 replicate 3), respectively. A total of 24 samples were sequenced and analyzed.
RNA extraction, library construction and sequencing
Three replicates were included at each time point, and each sample included 0.6 g eggs or larvae. RNA extraction, library construction, and RNA-sequencing were carried out by the Gene Denovo Biotechnology Co. (Guangzhou, China). Total RNA was extracted using TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer’s protocols. The integrity and quality of RNA were assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) and checked using RNase free agarose gel electrophoresis. For strand-specific cDNA libraries, rRNAs were removed to retain mRNAs and ncRNAs using Ribo-Zero rRNA removal kit (H/M/R) (Illumina, San Diego, USA) according to the manufacturer’s protocols. The enriched mRNAs and ncRNAs were fragmented into short fragments by using fragmentation buffer and reverse transcribed into cDNA with random primers. Second-strand cDNA were synthesized by DNA polymerase I, RNase H, dNTP (dUTP instead of dTTP) and buffer. Next, the cDNA fragments were purified with QIAquick polymerase chain reaction (PCR) extraction kit (Qiagen, Venlo, The Netherlands), end repaired, A base added, and ligated to Illumina sequencing adapters. Then UNG (Uracil-N-Glycosylase) were used to digest the second-strand cDNA. The digested products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced to perform pair-end sequencing using Illumina HiSeq™ 4000 by Gene Denovo Biotechnology Co. (Guangzhou, China). For small RNA cDNA library, after total RNA was extracted by Trizol reagent kit (Invitrogen, CA, USA), the RNA molecules in a size range of 18–30 nt were purified from total RNA on a urea polyacrylamide gel electrophoresis (PAGE). The purified small RNA was ligated with 3’ RNA adaptor to add a 3’ hydroxyl group. The 5’ adapters were then ligated to the RNAs as well. Reverse transcription was done to convert the RNA to cDNA, which was then selectively enriched by 12 cycles of PCR. The PCR products were purified on PAGE in a size range of 140 to 160 bp. Libraries were sequenced to perform single-end sequencing using Illumina HiSeq TM 2500 by Gene Denovo Biotechnology Co. (Guangzhou, China).
For mRNAs and lncRNAs, reads obtained from the sequencing machines included raw reads containing adapters or low-quality bases which would affect the following assembly and analysis. Thus, to get high quality clean reads, reads were further filtered by fastp (version 0.18.0) . The parameters were as follows: 1) removing reads containing adapters; 2) removing reads containing more than 10% of unknown nucleotides (N); 3) removing low quality reads containing more than 50% of low quality (Q-value ≤ 20) bases. High-quality reads were aligned with reference N. bombycis (http://silkpathdb.swu.edu.cn/datasets) or silkworm (https://silkdb.bioinfotoolkits.net /main/species-info/-1) genome using HISAT2 (version 2.1.0) . Then, guided by the Ensembl gene annotation file, the StringTie software program  was used to the reconstruction of transcripts. Differential gene expression and transcript abundance (expressed as fragment per kilobase of transcript per million mapped reads (FPKM) values) were calculated using the RESM program . Accordingly, absolute a fold change ≥ 2 and a P value < 0.05 were calculated based on FPKM, and differentially expressed (DE) mRNAs and lncRNAs were identified. All DE mRNAs were mapped to GO terms (http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway [37,38,39]. LncRNA target genes were predicted by locations in relation to nearby genes, and KEGG pathway analyses were carried out on these target genes [37,38,39]. In addition, we performed lncRNA antisense and trans-regulation analysis.
For circRNAs, 20 mers from both ends of the unmapped reads were extracted, and aligned to the reference genomes to identify unique anchor positions within the splice site using the methods described by Memczak et al. , respectively. Then anchor reads that aligned in the reverse orientation indicated circRNA splicing was subjected to Find_circ software  to identify circRNAs. To identify differentially expressed circRNAs across samples or groups, the edgeR package (http://www.rproject.org/) was used. We identified circRNAs with a fold change ≥ 2 and a P value < 0.05 in a comparison between samples or groups as significant differentially expressed circRNAs. KEGG pathway analyses were performed for the DE circRNA-associated genes [37,38,39].
For miRNAs, reads obtained from the sequencing machines included dirty reads containing adapters or low-quality bases which would affect the following assembly and analysis. Low-quality sRNA reads were first filtered form the raw sequences. Several types of impurities such as 3’ adaptor null reads, insert null reads, and 5’ adaptor contaminants, were removed from the remaining sRNA reads. Subsequently, reads smaller than 18 nt or having Poly (A) tails were removed from the high-quality sRNA reads. Lastly, all of the clean tags were aligned with small RNAs in GenBank database (Release 209.0)  to identify and remove rRNA, scRNA, snoRNA, snRNA and tRNA. Meanwhile all of the clean tags were aligned with small RNAs in Rfam database (Release 11.0)  to identify and remove rRNA, scRNA, sonRNA, snRNA and tRNA. All of the clean tags were also aligned with reference genomes. Those mapped to exons or introns might be fragments from mRNA degradation, so these tags were removed. The tags mapped to repeat sequences were also removed. All of the clean tags were then searched against miRbase database (Release 22)  to identify exist miRNAs or known miRNAs. Additionally, all of the unannotated data were aligned with N. bombycis or silkworm genome. According to their genome positions and hairpin structures predicted by software MiReap_v0.2, the novel miRNA candidates were identified. We identified miRNAs with a fold change ≥ 2 and a P value < 0.05 in a comparison as DE miRNAs. Based on the sequences of the miRNAs, the candidate target genes were predicted by RNAhybrid , Miranda  and TargetScan  with default settings. KEGG pathway analyses were performed for the DE miRNA-targeted genes [37,38,39].
Analysis of the ceRNA (lncRNA/circRNA) regulatory network
CeRNA regulatory networks were constructed based on the overlapping miRNA target between the mRNA and lncRNA/circRNA. The nodes in the networks include lncRNAs/circRNA, miRNAs, and mRNAs. The lncRNAs/circRNA-miRNA-mRNA regulatory networks were visualized with Cytoscape (v3.6.0) ( http://www.cytoscape.org/).
Statistical significance was assessed by one-way ANOVA for comparisons among multiple groups. Statistical significance was defined as P value < 0.05.
Validation of the congenital N. bombycis infection in silkworm embryos and larvae
To confirm the congenital N. bombycis infection in silkworm, embryos of 5-day from the parent silkworms infected with N. bombycis were used to detect the microsporidia spores by H&E staining. A few N. bombycis spores appeared at the region of the yolk granules (Fig. S1A). Histological examination of larvae of silkworm reared in isolation on 1, 5 and 10-day revealed that a number of spores were present in various cells such as muscle cells, adipose cells and epidermis cells (Fig. S1B-D). The presence of N. bombycis spores in silkworms derived from congenital infection was further confirmed by TEM. Consistent with H&E staining, microsporidia spores were mainly around the yolk granules in embryos of 5-day as revealed by TEM (Fig. S1E). On larvae of 1-day, spores were predominantly present in the midgut of silkworms (Fig. S1F) and spread to fat body and skin tissues on 5 and 10-day larvae (Fig. S1G, H). In contrast, no spores were observed in the 5-day embryos and 1, 5, 10-day larvae from parent silkworms uninfected with N. bombycis.
The summary of the whole transcriptome sequencing
To investigate the RNA transcripts of host and pathogen at different developmental stages, the specimen including 5-day infected and uninfected embryos, 1, 5, 10-day infected and uninfected larvae were collected and processed for RNA-seq analysis using N. bombycis and B. mori genome as references. For strand-specific cDNA libraries, each of 24 samples tested in the platform produced an average of 27.43 Gb of data. In the infected eggs and larvae, with 1.55% to 6.91% of these reads mapped to the N. bombycis genome, and in the uninfected eggs and larvae, there was negligible reads that can be mapped to the N. bombycis genome (Table 1). For small RNA libraries, the sequencing analysis also found negligible tags that can be mapped to the N. bombycis genome in the uninfected eggs and larvae, and 23.02% to 56.73% of tags can be mapped to the N. bombycis genome in the infected eggs and larvae (Table 2). These data showed that the experimental materials and whole transcriptome sequencing have good reliability.
Transcriptome-wide identification of N. bombycis mRNAs, lncRNAs, circRNAs and miRNAs in congenitally infected silkworm embryos and larvae
As shown in Fig. S2, we performed principal component analysis (PCA) based on the original mRNA, lncRNA, circRNA, and miRNA data. For mRNA data, PCA demonstrated a clear separation between infected groups and uninfected groups at the four different stages, and the distance between infected and uninfected groups at the fifth day of embryos of silkworms was relatively short (Fig. S2A). Similarly, based on lncRNA data, the infected groups and uninfected samples at different stages were distinguished by PCA (Fig. S2B). For circRNA data, PCA could not show an unambiguous separation between infected groups and uninfected groups at the fifth day of embryos and first day of larvae of silkworms, and a clear separation between infected groups and uninfected groups was present at the fifth and tenth day of larvae (Fig. S2C). For miRNA data, PCA showed that the overall distance between different groups was relatively long, indicating that there were major differences among samples (Fig. S2D).
We subsequently analyzed four RNA molecules (mRNA, lncRNA, circRNA and miRNAs) of N. bombycis. For N. bombycis mRNAs, a total of 4,051 mRNAs were aligned to reference genes and 104 mRNAs were defined as novel genes (Table S1). Additionally, the total numbers of the identified N. bombycis mRNA genes at different developmental stages were comparable with 3,782, 3,881, 4,024, and 3,960 mRNA genes in 5-day embryos (T1 stage), 1-day (T2 stage), 5-day (T3-stage), 10-day (T4 stage) larvae, respectively (Table 3), and 3,640 genes were identified in all four above mentioned stages and only 19, 22, 65, and 41 genes were stage-specifically expressed genes based on Venn diagram analysis (Fig. 1A and Table S2), indicating that N. bombycis mRNAs expression was not substantially affected by host growth and development.
With regard to the N. bombycis lncRNAs, 403 putative novel lncRNAs were identified that can be classified into five categories according to their genomic locations relative to protein coding genes: 8 sense lncRNAs, 166 antisense lncRNAs, 17 bidirectional lncRNAs, 168 intergenic lncRNAs and 44 other lncRNAs. There were 362, 380, 392, and 388 lncRNAs that were predicted in T1 (5-day embryos), T2 (1-day larvae), T3 (5-day larvae), and T4 (10-day larvae) stages, respectively (Table 3). When comparing lncRNAs profiles at different stages, we found that only two, one, six and three lncRNAs were uniquely expressed at T1, T2, T3 and T4 stage, respectively, and 352 lncRNAs were shared by all four stages, which account for 87.6% of all predicted N. bombycis lncRNAs (Fig. 1B).
A total of 62 novel circRNAs were predicted in the N. bombycis specimens collected in our study, of which 25, 31, 46, and 43 circRNAs were detected in T1 (5-day embryos), T2 (1-day larvae), T3 (5-day larvae), and T4 (10-day larvae) stages, respectively (Table 3). For circRNA type, 62 circRNAs can be classified into five categories: 9 antisense circRNAs, 16 exon_intron circRNAs, 12 intergenic circRNAs and 25 one_exon circRNAs. 62 circRNAs remained with sizes ranging from 129 to 37,391 bp. Interestingly, only 17 N. bombycis circRNAs were identified in all four stages, accounting for 27.4% of all predicted N. bombycis circRNAs identified in this study. Zero, four, ten, and eight circRNAs were uniquely expressed in T1, T2, T3, and T4 stages, respectively (Fig. 1C), suggesting that the expression pattern of circRNAs was different from that of mRNAs and lncRNAs, and N. bombycis circRNAs expression was more likely to be affected by host development pressure.
Regarding miRNA, a total of 284 miRNAs were identified based on the N. bombycis genome, and there were 277, 244, 273, and 277 miRNAs that were detected in T1 (5-day embryos), T2 (1-day larvae), T3 (5-day larvae), and T4 (10-day larvae) stages, respectively (Table 3), with 241 miRNAs being identified in all four stages (Fig. 1D). Consistent with expression profile of mRNAs and lncRNAs, only a small fraction of miRNA is stage specific. One, zero, three, and two miRNAs were uniquely expressed in T1, T2, T3, and T4 stages, respectively. In addition, all 284 miRNAs can be classified into intergenic miRNAs.
Expression profiles analysis of N. bombycis mRNAs, lncRNAs, circRNAs and miRNAs in congenitally infected silkworm embryos and larvae
To identify expression pattern of N. bombycis coding and non-coding RNAs in four different stages of congenitally infected silkworm, differentially expressed (DE) RNA between four stages were characterized based on the combined criteria of absolute fold change ≥ 2 and a P value < 0.05. As shown in Fig. 2, different subsets of RNAs were up or down- regulated in N. bombycis at different stages. Many differentially expressed mRNAs, lncRNAs and miRNAs were found, however rare differentially expressed circRNAs were shown at the different developmental stages.
To further visualize the expression of DE mRNAs, we generated heat maps with down-regulated RNAs showing in blue and up-regulated RNAs showing in red. As shown in Fig. 3A, there were 574 DE genes that were fallen into five clusters according to the different gene expression patterns (Table S3). KEGG pathways were enriched by up-and down-regulated DE genes. The expression of genes in cluster I rapidly increased from the fifth day of embryos to the first day of larvae, and gradually decreased. Genes were associated with ribosome biogenesis in eukaryotes, protein processing in endoplasmic reticulum, aminoacyl-tRNA biosynthesis, RNA polymerase, and pyrimidine metabolism (Fig. S3A). Genes in cluster II continued to increase from the fifth day of embryos to the fifth day of larvae, and remained temporally stable at the tenth day of larvae. These genes were mainly involved in purine metabolism, RNA polymerase, pyrimidine metabolism, non-homologous end-joining, and SNARE interactions in vesicular transport (Fig. S3B). The expression of genes in cluster III was highest at the fifth day of embryos, and gradually decreased in larvae. KEGG pathways analysis showed that genes were mainly involved in starch and sucrose metabolism, biosynthesis of antibiotics, metabolic pathways, carbon metabolism, and microbial metabolism in diverse environments (Fig. S3C). The expression of genes in cluster IV decreased from the fifth day of embryos to the first day of larvae, and occurred a same trend from fifth day of larvae to tenth day of larvae. As shown in Fig. S3D, these genes were mainly involved in inositol phosphate metabolism, fructose and mannose metabolism, glycolysis/gluconeogenesis, biosynthesis of amino acids, and protein export. Genes in cluster V continued to increase their highest expression at the tenth day of larvae, and these genes were associated with biosynthesis of antibiotics, metabolic pathways, MAPK signaling pathway-yeast, sphingolipid metabolism, and amino sugar and nucleotide sugar metabolism (Fig. S3E).
Similarly, we generated a heat map based on 55 DE lncRNAs of N. bombycis (Fig. 3B, Table S3). We found that most lncRNAs were highly expressed in 10-day larvae, and those lncRNAs can be classified into four clusters. The expression of lncRNAs in cluster I was highly expression at the fifth day of embryos and the tenth day of larvae, and remained low expression at the other stages. LncRNAs in cluster II continued to increase their highest expression at the tenth day of larvae. Sever lncRNAs (MSTRG10794.1, MSTRG5153.1, MSTRG5153.2, MSTRG5153.3, MSTRG12418.1, MSTRG2172.1, MSTRG2172.2, and MSTRG2172.3) in cluster III was highest at the fifth day of embryos, and gradually decreased in larvae. The expression of lncRNAs in cluster IV increased from the fifth day of embryos to the first day of larvae, and occurred a same trend from fifth day of larvae to tenth day of larvae. Furthermore, lncRNA-mRNA association analysis revealed target genes of lncRNAs in antisense, cis, and trans-regulation (Table S4).
The heat map showed that the expression levels of most DE miRNAs in larvae were higher than those in embryos, and 125 DE miRNAs (Table S3) target genes were divided into 4 clusters (Fig. 3C, Table S5). The expression of miRNAs in cluster I increased from the fifth day of embryos to the fifth day of larvae, and showed downregulation at the tenth day of larvae. These miRNAs were mainly involved in mRNA surveillance pathway, cell cycle-yeast, aminoacyl-tRNA biosynthesis, ether lipid metabolism, and glycerophospholipid metabolism (Fig. S4A). The expression of miRNAs in cluster II was highest at the tenth day of larvae, and showed low expression at the other stages. In second cluster miRNAs targeted genes were associated with protein processing in endoplasmic reticulum, glycerophospholipid metabolism, RNA polymerase, glycine, serine and threonine metabolism, and methane metabolism (Fig. S4B). The expression of miRNAs in cluster III continued to increase from the fifth day of embryos to the fifth day of larvae, and remained temporally stable at the tenth day of larvae. As shown in Fig. S4C, miRNAs targeted genes in cluster III were associated with aminoacyl-tRNA biosynthesis, protein processing in endoplasmic reticulum, starch and sucrose metabolism, proteasome, and glycolysis /gluconeogenesis. The expression of miRNAs in cluster IV continued to increase their highest expression at the tenth day of larvae. As shown in Fig. S4D, miRNAs targeted genes in cluster IV were involved into protein processing in endoplasmic reticulum, cell cycle-yeast, aminoacyl-tRNA biosynthesis, glycolysis/gluconeogenesis, and arachidonic acid metabolism.
Intraspecies ceRNA networks of N. bombycis lncRNAs and circRNAs in congenitally infected silkworm embryos and larvae
To discover the function of lncRNAs and circRNAs of N. bombycis ceRNA networks for DE lncRNAs and circRNAs were constructed. We first apply the principle of base complementary pairing to match the miRNAs with mRNAs, lncRNAs and circRNAs, resulting in values of the correlation between miRNA-mRNA, miRNA-lncRNA or miRNA-circRNA node. Then, only the negatively correlated pairs of miRNA-mRNA, miRNA-lncRNA or miRNA-circRNA were retained, and the final ceRNA networks included lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA networks were constructed (Table S6, 7), and this ceRNA networks include 18 lncRNAs, one circRNA, 28 miRNA, and 113 mRNA were visualized by Cytoscape software (Fig. 4A). Meanwhile, Sankey diagrams showed the ceRNA networks involving 14 key parasites genes such as PTP3, ricin-B-lectin, SWP4, HSP90, 18 lncRNAs, one circRNA, and 20 miRNAs (Fig. 4B). Furthermore, we used topology method to calculate the connection degree of each gene to illustrate its importance in the ceRNA network (Table S8). Three miRNAs (novel-m0224-3p, novel-m0211-3p and novel-m0065-3p), three lncRNAs (MSTRG.5153.1, MSTRG.5153.2 and MSTRG.5153.3), and one circRNA (novel-circ-000006) had high connection degrees.
Transcriptome-wide identification and characteristics of B. mori mRNAs, lncRNAs, circRNAs and miRNAs in the congenital N. bombycis-infected silkworm embryos and larvae
For B. mori mRNA, 13,592 genes were aligned to reference genes and 937 genes were defined as novel genes (Table S9). Furthermore, the total numbers of the predicted B. mori genes at different stages were comparable with 13,571, 13,431, 13,232, and 13,213 genes in 5-day embryos (CK1 and T1), 1-day (CK2 and T2), 5-day (CK3 and T3), 10-day (CK4 and T4) larvae, respectively (Table 4), and common 11,722 genes were detected in four stages, accounting for 78.7% of all predicted B. mori genes and only 440, 228, 145, and 199 genes were stage-specifically expressed mRNA genes based on Venn diagram analysis (Fig. 5A).
Regarding lncRNA, a total of 3,038 B. mori lncRNAs were identified in the uninfected and infected embryos and larvae, and there were 2,917, 2,988, 2,991, and 2,987 lncRNAs that were detected in 5-day embryos (CK1 and T1), 1-day (CK2 and T2), 5-day (CK3 and T3), 10-day (CK4 and T4) larvae, respectively (Table 4). Additionally, there were 2,836 lncRNAs that were identified in all four stages, and only 11, 1, 4, and 1 lncRNAs were stage-specifically expressed lncRNA genes (Fig. 5B).
With regard to the B. mori circRNAs, a total of 19,039 lncRNAs were identified in the uninfected and infected embryos and larvae, and there were 11,512, 13,615, 12,522, and 12,623 circRNAs that were identified in CK1 and T1, CK2 and T2, CK3 and T3, and CK4 and T4 stages, respectively (Table 4) with 6,527 circRNAs being identified in all four stages and 1,275, 1,439, 1,094, and 1,155 circRNAs being identified in CK1 and T1, CK2 and T2, CK3 and T3, and CK4 and T4 stages, respectively (Fig. 5C).
Similarly, as shown in Table 4, a total of 3,403 B. mori miRNAs were identified in all four above mentioned stages. There were 2,890, 2,326, 1,736, and 1,534 miRNAs that were detected in CK1 and T1, CK2 and T2, CK3 and T3, and CK4 and T4 stages, respectively (Table 4). Additionally, there were 1,114 miRNAs were shared in four stages, accounting for 32.6% of all predicted B. mori miRNAs, and 697, 199, 103, and 107 miRNAs that were stage-specifically expressed miRNA genes based on Venn diagram analysis (Fig. 5D).
Analysis of differentially expressed mRNAs, lncRNAs, circRNAs and miRNAs of B. mori in the congenital N. bombycis-infected silkworm embryos and larvae
To understand the role of various forms of host RNAs in N. bombycis congenital infection, we next analyzed the differentially expressed mRNAs, lncRNAs, circRNAs, and miRNAs by comparing uninfected and congenitally infected specimen. We found that there are a large number of differentially expressed mRNAs, lncRNAs, circRNAs and miRNAs in the embryos and larvae congenitally infected by N. bombycis. As shown in Fig. 6A, 460 (205 up-regulated and 255 down-regulated), 377 (182 up-regulated and 195 down-regulated), 628 (359 up-regulated and 269 down-regulated), and 864 mRNAs (538 up-regulated and 326 down-regulated) were DE genes during the congenital N. bombycis infection in the embryos on day 5, larvae on day 1, 5, and 10, respectively. For lncRNAs, there were 77 (31 up-regulated and 46 down-regulated), 87 (39 up-regulated and 48 down-regulated), 140 (74 up-regulated and 66 down-regulated), and 218 DE lncRNAs (144 up-regulated and 74 down-regulated) after the congenital N. bombycis infection in the embryos on day 5, larvae on 1, 5 and 10 days, respectively (Fig. 6B). Regarding circRNAs, there were 171 (88 up-regulated and 83 down-regulated), 237 (105 up-regulated and 132 down-regulated), 234 (134 up-regulated and 100 down-regulated) and 307 (195 up-regulated and 112 down-regulated) DE circRNAs in the silkworm 5-day embryos, 1, 5 and 10-day larvae on of the congenital infected group versus the control group, respectively (Fig. 6C). As shown in Fig. 6D, in the silkworm embryos on 5 days, 358 DE miRNAs (162 up-regulated and 196 down-regulated) were observed in the silkworm embryos following the congenital N. bombycis infection; there were 351 DE miRNAs (108 up-regulated and 243 down-regulated), 166 DE miRNAs (32 up-regulated and 134 down-regulated), 139 DE miRNAs (53 up-regulated and 86 down-regulated) in the larvae on 1, 5 and 10 days of the infected group versus the control group, respectively.
Moreover, we analyzed exclusive and shared differentially expressed mRNAs, lncRNAs, circRNAs and miRNAs over all the stages and stage-specific differentially expressed mRNAs, lncRNAs, circRNAs and miRNAs (Fig. 7A-D). Figure 7 represented that there were three mRNAs (protein lethal (2) essential for life-like, nuclear receptor GRF and transferrin precursor), one lncRNA (MSTRG.23380.1), no circRNA and six miRNAs (bmo-miR-2780a-5p, bmo-miR-2780b, bmo-miR-3001, bmo-miR-306b, miR-8517-x and novel-m0636-3p) that were common differentially expressed genes over all the stages.
Functional analysis of congenital N. bombycis infection-induced DE mRNAs, lncRNAs, circRNAs and miRNA in the silkworm embryos and larvae
To reveal the potential functions of DE mRNAs, we analyzed KEGG pathways of those up-regulated and down-regulated mRNAs (Fig. S5A-D). KEGG enrichment analysis of DE mRNAs showed that those genes involved in cell cycle, oocyte meiosis, Toll and Imd signaling pathway, longevity regulating pathway-multiple species, and antigen processing and presentation. In order to investigate the immune status of the host, B. mori immune function related DE genes were constructed to form four heat maps (Fig. 8A-D). We found that immune recognition receptors (beta-1,3-glucan recognition protein 2 precursor (βGRP2), fibrinogen alpha chain, and c-type lectin 4), immune regulatory genes (ovarian serine protease and serpin-15), immune signaling pathway gene (Spz3) and immunoeffector genes (pro-phenol oxidase and iNOS) were down-regulated in 5-day embryos exposed to the congenital N. bombycis challenge, while those genes and other immune-related genes such as peptidoglycan recognition protein, scavenger receptor, serine protease 7, dorsal and most antimicrobial peptide genes were up-regulated in larvae, especially in 5-day larvae exposed to the congenital N. bombycis challenge. The above results suggest that immune system involved in Toll and Imd signaling pathway and hemolymph melanization is in immunosuppression in the embryos exposed to the congenital N. bombycis infection, and host immune system involved in Toll and Imd signaling pathway, hemolymph melanization, and phagocytosis is activated by N. bombycis infection in the larvae, especially in the early and middle larvae.
One of the functions of lncRNAs is cis-regulation of their neighboring genes on the same allele. LncRNAs in less than 100 Kb up/downstream of a gene were likely to be cis-regulators. We identified lncRNA-mRNA pairs including 99 lncRNAs and 105 mRNAs (Table S10). The cis-target genes of lncRNAs were then subjected to enrichment analysis of KEGG pathways. KEGG enrichment analysis showed that systemic lupus erythematosus, Toll and Imd signaling pathway, drug metabolism-other enzymes, longevity regulating pathway-multiple species, caffeine metabolism, and alcoholism were enriched under the congenital N. bombycis stress (Fig. S6). We specifically focused on lncRNAs regulating immune-related gene expression in cis (Fig. 9). A string of tandem genes encoded cecropin and enbocin, which were lncRNAs (MSTRG18915.1, MSTRG18915.2 and MSTRG18915.3) nearly genes (Fig. 9A). In addition, one gene (moricin) was predicted to be three lncRNAs (MSTRG15336.1, MSTRG15337.1 and MSTRG15338.1) nearly genes that were likely to be regulated by lncRNAs in cis (Fig. 9B).
We performed the functional enrichment analysis of DE circRNA source genes to study the main functions of these genes (Fig. S7A-D). KEGG enrichment analysis of DE circRNA source genes showed that those genes involved in relaxin signaling pathway, PI3K-Akt signaling pathway, focal adhesion, secretion and action, microRNAs and thyroid hormone synthesis.
MiRNAs recognize and bind to partially complementary sites in the 3′ untranslated regions of target genes and regulate the expression of the target transcripts. To study the functions of the targeted genes, we carried out the functional enrichment analysis of DE miRNA target genes (Fig. S8A-D). KEGG enrichment analysis of DE miRNA target genes showed many important biological processes including protein processing in endoplasmic reticulum, apoptosis-multiple species, endocytosis, VEGF signaling pathway and cell cycle were enriched. In addition, the basic metabolisms of life such as carbohydrate, lipid and amino acid metabolisms were also targets of miRNAs.
Intraspecies ceRNA networks of DE lncRNAs and circRNAs in B. mori in the congenitally infected by N. bombycis
It was known that lncRNAs and circRNAs can be considered as two types of ceRNAs, which play significant role in response to various pathogens infection [23, 46]. To access the function of lncRNAs and circRNAs in the congenital N. bombycis infection in the silkworm embryos and larvae, dynamical ceRNA networks were analyzed for four different stages during the congenital N. bombycis infection. In the 5-day silkworm embryos, lncRNA and circRNA ceRNA networks were constructed following the congenital N. bombycis infection (Table S11, 12). The lncRNA ceRNA network contained 59 lncRNAs and 103 miRNAs, and the circRNA ceRNA network consisted of 64 circRNAs and 101 miRNAs. Similarly, in the silkworm 1, 5 and 10-day larvae, the lncRNA/circRNA ceRNA networks were constructed following the congenital N. bombycis infection, respectively (Table S13-18). According to previous reports, ceRNA network has two characteristically distinct types: one is common to all stages, which is defined as the common ceRNA networks, and other is unique for each stage, which is defined the unique ceRNA networks . We found that lncRNA/circRNA ceRNA networks are almost the unique ceRNA networks.
The innate immune response of silkworms plays crucial roles in the resistance to N. bombycis infection ; therefore, we next focused on the differential expression of immune-related genes in the silkworm embryos and larvae during the congenital N. bombycis infection. There were 8 immune-related DE mRNAs including Toll-6, fibrinogen alpha chain, Serpin-6, Caspase-8, Caspase-1, iNOS, hemolin-interacting protein (Hemolin-IP), Pellino, and TBC1 domain family member 14 (TBC1d14). Based on 8 DE mRNAs, we further constructed and illustrated the ceRNA networks. The lncRNA-miRNA-mRNA pathways were constructed, including 140 lncRNAs and 6 miRNAs (Fig. 10A). And we used topology method to calculate the connection degree of each gene to illustrate its importance in the lncRNA ceRNA network (Table S19). Similarly, the circRNA-miRNA-mRNA pathways were constructed, including 5 circRNAs and 4 miRNAs (Fig. 10B), of which, lncRNA ceRNA networks shares the 4 mRNAs and 3 miRNAs with circRNA ceRNA networks. Furthermore, we used topology method to calculate the connection degree of each gene to illustrate its importance in the circRNA ceRNA network (Table S20).
Intraspecies and cross-species miRNAs targets analysis
It is known that miRNAs act not only in the cells where they are biologically formed but they also function in cells from different species . We screened 125 parasite and 1,376 silkworms differentially expressed miRNAs, and their intraspecies and cross-species miRNAs targets were predicted in N. bombycis and silkworm interaction (Fig. 11). For 125 parasite miRNAs, intraspecies and cross-species miRNA targets analysis showed that 4,378 parasite genes and 6,951 host genes were predicted to be those parasite miRNAs targets, respectively. For 1,376 host miRNAs, 5,283 silkworm genes and 595 parasite genes were predicted to be those host miRNAs targets, respectively. Moreover, 5,176 host genes and 454 parasite genes may be regulated not only by pathogen miRNAs, but also by host miRNAs. KEGG pathways analysis showed that 5,176 host genes were mainly involved in cellular processes, environmental information processing, genetic information processing, organismal systems, metabolism, and human disease, and 454 parasite genes were involved in cellular processes, metabolism, environmental information processing and genetic information processing. In this analysis, miRNA-mediated cross-species regulation is worth further investigation. For example, two N. bombycis miRNAs (novel-m0083-3p and novel-m0044-3p) are highly expressed, and their targets are silkworm Pelle and Spz3, two key genes in Toll signaling pathway, which are downregulated in embryos after N. bombycis congenital infection. We assume that microsporidia deploy their miRNAs to attenuate silkworm immune responses by suppressing the expression of the related immune genes. In addition, one N. bombycis miRNA (novel-m0172-3p) is possible to target ATP binding cassette subfamily D member 3 (ABCD3) to inhibit β-oxidation of fatty acid to reduce lipid consumption for parasite multiplication. It seems that cross-species miRNAs play an important role in host-microsporidia interactions.
Microsporidia, naturally occurring obligate intracellular fungal-related parasites, have been at the forefront of efforts to develop biocontrol agents for agriculture and forestry pests due to their host specificity and transovarial transmission [50,51,52]. As a typical species of microsporia, N. bombycis is an important pathogen of silkworm B. mori causing the epidemic disease pébrine through horizontal and vertical transmission. N. bombycis may be a model of congenital infection for studying microsporidia pathogenicity and pathogen-host interactions. In microsporidia-host interplays, insects recognize pathogens by its innate immune system, and parasites attenuated host immune responses through various molecular mechanisms . However, these underlying mechanisms remain elusive. Here, we showed the pan-transcriptome including the coding RNA and ncRNAs for N. bombycis and host B. mori with a specific focus on microsporidia congenital infection, and suggest that ncRNA-mediated regulation plays a vital role in the microsporidia congenital infection, which provides solid foundation to study fungal pathogenicity and pathogen-host interactions in embryos and larvae.
Non-coding RNAs comprise various RNA species, including miRNA, lncRNA and circRNA, which are recognized as key regulators of gene expression to regulate various physiological processes, including growth, development, metabolism and immune responses [54,55,56]. Predictably, we obtained 4,046 reference genes and 125 novel genes compared with the 4,460 annotated genes in N. bombycis (http://silkpathdb.swu.edu.cn). Furthermore, in our study we identified that 403 lncRNAs and 284 miRNAs in N. bombycis, which are higher abundance compared to the corresponding species reported for N. ceranae [20, 32, 57]. Recent studies suggest that N. bombycis modulates small RNA (sRNA)-mediated regulation that increases fungal pathogenicity .
We found that a large number of mRNAs were identified in N. bombycis at various stages of infection, only a handful of genes were stage-specifically expressed genes for each stage. Similarly, there were 351 lncRNAs (87.34% of all predicted N. bombycis lncRNAs) and 241 miRNAs (84.86% of all predicted N. bombycis miRNAs) that were identified in all four above mentioned stages. However, regarding circRNAs, it is not the same case, only 17 of 62 circRNAs were detected in all four stages of infection, suggesting that the major part of circRNAs encoded from N. bombycis exhibit specific development or stress-inducible expression patterns. CircRNAs usually show specific cell-type, tissue, and developmental expression patterns in animals and plants [58, 59]. CircRNAs also exhibit abiotic and biotic stress-inducible expression patterns [60, 61].
As expected, we found that a variety of DE mRNAs, lncRNAs, and miRNAs encoded by N. bombycis in the process of N. bombycis infecting silkworm embryos and larvae, suggesting that those mRNAs, lncRNAs, and miRNAs may play a critical role in the biological processes of N. bombycis in silkworms. In contrast, only three N. bombycis circRNAs were differentially expressed during the process. Interestingly, PTP3 was significantly downregulated in 1, 5, 10-day larvae infected with N. bombycis compared with 5-day infected embryos. Previous studies have shown that PTP3 is involved in the sporoblast-to-spore polar tube biogenesis and plays a significant role in the control of the polar tube extrusion in E. cuniculi . Most importantly, some parasites lncRNAs (MSTRG.7750.2, MSTRG.10054.3, MSTRG.11899.2 and MSTRG.11875.2) were predicted to regulate PTP3 expression by sponging miRNAs. Additionally, two parasite lncRNAs (MSTRG.585.1 and MSTRG.1988.2) were predicted to regulate their neighbor gene PTP3 expression in cis, because lncRNAs in less than 100 Kb up/downstream of a gene were likely to be cis-regulators. Interestingly, 32 down-regulated miRNAs of N. bombycis were detected in 1-day larvae, compared to 5-day embryos infected with N. bombycis. Those highly expressed miRNAs of N. bombycis in 5-day embryos may regulate host genes expression to evade the host innate immune response, resulting in contribute to the congenital infection of N. bombycis in silkworm embryos. For example, one parasite miRNA (novel-m0044-3p) is highly expressed, and its target Toll ligand Spz3 expression is decreased after N. bombycis infection in the embryos. Recent studies showed that Spz3 is mainly involved in the melanization process in the stripe pattern formation in caterpillars . Another melanization factor pro-phenol oxidase plays an important role in defense against intruding parasitoids , and is downregulated after N. bombycis infection in embryos. These data suggest that microsporidia have evolved to inject parasites miRNAs into the host cells to suppress hemolymph melanization in embryos for the successful development of their progeny. This is an example of the evolution of a cross kingdom miRNA-mediated defense mechanism in microsporidia.
Salmena et al. first proposed that ceRNA cross-talk between RNAs, both coding and non-coding, through miRNA response elements (MREs), forming large-scale regulating network across the transcriptome . CeRNA networks link the function of protein-coding mRNAs with that of non-coding RNAs such as miRNA, lncRNA, and circRNA. In this study, based on known functional genes of N. bombycis, we constructed N. bombycis lncRNA/ circRNA-miRNA-mRNA ceRNA networks by bioinformatics analysis. We found that one miRNA can bind to different lncRNAs to regulate the same gene. For example, one miRNA (novel-m0174-5p) is predicted to bind to three lncRNAs (MSTRG.2013.1, MSTRG.2923.1 and MSTRG.4851.1) to regulate transposase expression. One lncRNA was predicted to bind to multiple miRNAs to regulate the same expression. For instance, one lncRNA (MSTRG.11875.2) can bind to two miRNAs (novel-m0206-5p and novel-m0224-3p) to regulate PTP3 expression. To our knowledge, this is the comprehensive analysis of microsporidia ncRNAs available, and it will be interest in the future to investigate the function of key lncRNAs, circRNAs, and miRNAs and lncRNA/circRNA-related ceRNA networks in N. bombycis.
Regarding silkworms, we found that a total of 14,889 mRNAs, 3,038 lncRNAs, 19,039 circRNAs, and 3,403 miRNAs were identified by the B. mori reference genome in the silkworm embryos and larvae with or without the congenital N. bombycis infection. Interestingly, for lncRNA, we identified 1,492 lncRNAs with more than 90% homologous sequences reported previously in B. mori . And for miRNA, we identified 427 miRNAs with 100% homologous sequences reported previously in B. mori . In addition, 19,039 circRNAs and many stage-specific circRNAs were identified, suggesting that circRNAs may play an important role in silkworm development, metabolism, and immune response system. Most importantly, we found that the number of miRNAs was significantly reduced by 50% after N. bombycis congenital infection. We also found that the number of small RNAs was twice as much as that of infected samples in small RNA libraries. The data suggest that microsporidia congenital infection greatly interferes with host miRNA biosynthesis.
In the present study, we detected many DE mRNAs during the congenital N. bombycis infection. This is consistent with our previous findings that silkworm mRNAs alter during the congenital N. bombycis infection in embryos and larvae. KEGG enrichment analysis of silkworm DE mRNAs showed that there were fewer immune-related genes in the embryos on day 5, and more immune-related genes involved in 1, 5, 10-day larvae exposed to the congenital N. bombycis challenge. Moreover, the expression of immune-related genes, such as βGRP2, Spz3 and pro-phenol oxidase, decreased in 5-day embryos after infected with N. bombycis compared with uninfected samples. While most immune genes, such as peptidoglycan recognition protein like (PGRP-L), Toll-like receptor 3 and antimicrobial peptide genes, were upregulated in 1, 5, 10-day larvae exposed to the congenital N. bombycis challenge compared with uninfected samples. The results further demonstrated that parasites miRNA-mediated regulation leads to host embryo immunosuppression to contribute to vertical transmission. Of course, the underlying mechanisms of immunosuppression or immune activation in the embryos and larvae exposed to the congenital N. bombycis infection remain to be determined. Furthermore, we predicted that there are many ncRNAs (lncRNAs, circRNAs and miRNAs) of silkworms during the congenital N. bombycis infection, and plays an important role in the fight against microsporidia infection. For example, we predicted that there were three silkworm genes (alpha-crystallin A chain, HSP25.4, and zinc transporter 2-like) that are likely to be regulated by one lncRNA (MSTRG.23380.1). HSP25.4, a member of the heat shock protein family, is highly expressed following NPV infection in Antheraea Pernyi and functions as an immune booster . We assume that lncRNA (MSTRG.23380.1) is possible to participate in the immune response to N. bombycis in silkworm by cis regulating the expression of HSP25.4. It also has been reported that N. ceranae-responsive lncRNAs may regulate the expression of neighboring genes by acting in cis fashion . For B. mori miRNAs, we found that there was significantly decreased in the number of miRNAs in the silkworm embryos and larvae exposed to the congenital N. bombycis challenge, compared with uninfected embryos and larvae. Both miRNA-33 and miRNA-122, which can target lipid metabolism genes to reprogram host lipid metabolism [68,69,70], are down-regulated during the congenital N. bombycis infection. miR-276 is associated with amino acid metabolism , which is down-regulated during the congenital N. bombycis infection. In addition, miR-306, a conserved miRNA, is worth mentioning as it showed that it is all differentially expressed genes in the silkworm exposed to the congenital N. bombycis infection in four different stages, compared with uninfected silkworm embryos and larvae. Some studies have been reported that miR-306 regulates host immune response in the insect during pathogens infection [72, 73]. In this study, we did not investigate translation capacity of these silkworm circRNAs, and regarding circRNA-protein interactions is elusive in microsporidia infection. Therefore, these ncRNAs may play an important role in the host–pathogen interaction, and the underlying regulatory mechanism requires further investigation.
Recent evidence has shown that lncRNA and circRNA, acting as ceRNA, have an important role in regulating miRNA target genes in host–pathogen interactions in insect. In the study, the lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA pathways were established, respectively. Moreover, based on eight immune-related genes, which was differentially expressed in silkworm exposed to the congenital N. bombycis challenge, we predicted lncRNA-miRNA-mRNA pathways and circRNA-miRNA-mRNA pathways, respectively. Two LncRNAs (MSTRG.20503.1 and MSTRG.9320.1) was predicted to modulate the TBC1d14 gene by acting as a decoy of miR-8523 during the silkworm immune response to the congenital N. bombycis stress. A circRNA (circ_002617) may function as a ceRNA to regulate the Toll-6 gene by acting as a decoy of bmo-miR-283-5p in immune response to the congenital N. bombycis infection. Therefore, the regulatory networks of the lncRNAs and circRNAs need to be proven in the future.
In summary, we analyzed the pan-transcriptome including the coding RNA (mRNA) and ncRNAs (lncRNA, circRNA and miRNA) for both parasites and their host silkworms at different development stages. This is the first pan-transcriptome and cross-talk of the coding and non-coding RNAs in the microsporidia: N. bombycis congenital infection of silkworm embryos and larvae, and shows that ncRNA-mediated regulation plays a vital role in the microsporidia congenital infection, which provides new insights into understanding the basic biology of microsporidia and their host–pathogen interaction.
Availability of data and materials
The raw data generated was submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/sra) with the BioProject accession number (PRJNA953616).
Long non-coding RNA
Covalently closed circular RNA
Polar tube protein 3
Spore wall protein 4
Heat shock protein 90
Inducible nitric oxide synthase
CLIP-domain serine protease
Competing endogenous RNA
Hematoxylin & eosin
Transmission electron microscopy
Fragment per kilobase of transcript per million mapped reads
Phosphate buffer saline
Polymerase chain reaction
DNA-directed RNA polymerases I and III subunit
- ABC transporter M2:
ABC transporter B family member 2
- ABC transporter M24:
ABC transporter B family member 24
Kyoto Encyclopedia of Genes and Genomes
Beta-1,3-glucan recognition protein 2
Peptidoglycan recognition protein like
TBC1 domain family member 14
ATP binding cassette subfamily D member 3
National Center for Biotechnology Information
Sequence Read Archive
Didier ES. Microsporidiosis: an emerging and opportunistic infection in humans and animals. Acta Trop. 2005;94(1):61–76.
Stentiford GD, Becnel JJ, Weiss LM, Keeling PJ, Didier ES, Williams BAP, Bjornson S, Kent ML, Freeman MA, Brown MJF, Troemel ER, Roesel K, Sokolova Y, Snowden KF, Solter. Microsporidia-emergent pathogens in the global food chain (Trends in Parasitology 32, 336–348; April 2, 2016). Trends Parasitol. 2016;32(8):657.
Vavra J, Lukes J. Microsporidia and “the art of living together.” Adv Parasitol. 2013;82:253–319.
Pan G, Bao J, Ma Z, Song Y, Han B, Ran M, Li C, Zhou Z. Invertebrate host responses to microsporidia infections. Dev Comp Immunol. 2018;83:104–13.
Perry RB, Ulitsky I. The functions of long noncoding RNAs in development and stem cells. Development. 2016;143(21):3882–94.
Yang B, Jia L, Ren H, Jin C, Ren Q, Zhang H, Hu D, Zhang H, Hu L, Xie T. LncRNA DLX6-AS1 increases the expression of HIF-1alpha and promotes the malignant phenotypes of nasopharyngeal carcinoma cells via targeting MiR-199a-5p. Mol Genet Genomic Med. 2020;8(1): e1017.
Yang C, Wu D, Gao L, Liu X, Jin Y, Wang D, Wang T, Li X. Competing endogenous RNA networks in human cancer: hypothesis, validation, and perspectives. Oncotarget. 2016;7(12):13479–90.
Gregory RI, Yan KP, Amuthan G, Chendrimada T, Doratotaj B, Cooch N, Shiekhattar R. The Microprocessor complex mediates the genesis of microRNAs. Nature. 2004;432(7014):235–40.
Kim VN. MicroRNA biogenesis: coordinated cropping and dicing. Nat Rev Mol Cell Biol. 2005;6(5):376–85.
Mah SM, Buske C, Humphries RK, Kuchenbauer F. miRNA*: a passenger stranded in RNA-induced silencing complex? Crit Rev Eukaryot Gene Expr. 2010;20(2):141–8.
Ouimet M, Koster S, Sakowski E, Ramkhelawon B, van Solingen C, Oldebeken S, Karunakaran D, Portal-Celhay C, Sheedy FJ, Ray TD, Cecchini K, Zamore PD, Rayner KJ, Marcel YL, Philips JA, Moore KJ. Mycobacterium tuberculosis induces the miR-33 locus to reprogram autophagy and host lipid metabolism. Nat Immunol. 2016;17(6):677–86.
Cui C, Wang Y, Liu J, Zhao J, Sun P, Wang S. A fungal pathogen deploys a small silencing RNA that attenuates mosquito immunity and facilitates infection. Nat Commun. 2019;10(1):4298.
Agrawal N, Dasaradhi PV, Mohmmed A, Malhotra P, Bhatnagar RK, Mukherjee SK. RNA interference: biology, mechanism, and applications. Microbiol Mol Biol Rev. 2003;67(4):657–85.
Huang Q. Evolution of Dicer and Argonaute orthologs in microsporidian parasites. Infect Genet Evol. 2018;65:329–32.
Dong Z, Zheng N, Hu C, Deng B, Fang W, Wu Q, Chen P, Huang X, Gao N, Lu C, Pan M. Nosema bombycis microRNA-like RNA 8 (Nb-milR8) increases fungal pathogenicity by modulating BmPEX16 gene expression in its host, Bombyx mori. Microbiol Spectrum. 2021;9(2): e0104821.
Hu C, Dong Z, Deng B, Wu Q, Chen P, Lu C, Pan M. MicroRNA-6498–5p Inhibits Nosema bombycis proliferation by downregulating BmPLPP2 in Bombyx mori. J Fungi. 2021;7(12):1051.
Kour S, Rath PC. Long noncoding RNAs in aging and age-related diseases. Ageing Res Rev. 2016;26:1–21.
Chen L, Zhou Y, Li H. LncRNA, miRNA and lncRNA-miRNA interaction in viral infection. Virus Res. 2018;257:25–32.
Jalali S, Bhartiya D, Lalwani MK, Sivasubbu S, Scaria V. Systematic transcriptome wide analysis of lncRNA-miRNA interactions. PLoS ONE. 2013;8(2): e53823.
Guo R, Chen D, Xiong C, Hou C, Zheng Y, Fu Z, Liang Q, Diao Q, Zhang L, Wang H, Hou Z, Kumar D. First identification of long non-coding RNAs in fungal parasite Nosema ceranae. Apidologie. 2018;49(5):660.
Wu Y, Cheng T, Liu C, Liu D, Zhang Q, Long R, Zhao P, Xia Q. Systematic identification and characterization of long non-coding RNAs in the silkworm, Bombyx mori. PLoS ONE. 2016;11(1): e0147147.
Chen D, Chen H, Du Y, Zhou D, Geng S, Wang H, Wan J, Xiong C, Zheng Y, Guo R. Genome-wide identification of long non-coding RNAs and their regulatory networks involved in Apis mellifera ligustica response to Nosema ceranae infection. Insects. 2019;10(8):245.
Satyavathi V, Ghosh R, Subramanian S. Long non-coding RNAs regulating immunity in insects. Non-Coding RNA. 2017;3(1):14.
Salzman J, Gawad C, Wang PL, Lacayo N, Brown PO. Circular RNAs are the predominant transcript isoform from hundreds of human genes in diverse cell types. PLoS ONE. 2012;7(2): e30733.
Jeck WR, Sorrentino JA, Wang K, Slevin MK, Burd CE, Liu J, Marzluff WF, Sharpless NE. Circular RNAs are abundant, conserved, and associated with ALU repeats. RNA. 2013;19(2):141–57.
Salzman J, Chen RE, Olsen MN, Wang PL, Brown PO. Cell-type specific features of circular RNA expression. PLoS Genet. 2013;9(9): e1003777.
Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, Kjems J. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495(7441):384–8.
Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, Loewer A, Ziebold U, Landthaler M, Kpcks C, le Noble F, Rajewsky N. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495(7441):333–8.
Ashwal-Fluss R, Meyer M, Pamudurti NR, Ivanov A, Bartok O, Hanan M, Evantal N, Memczak S, Rajewsky N, Kadener S. circRNA biogenesis competes with pre-mRNA splicing. Mol Cell. 2014;56(1):55–66.
Legnini I, Di Timoteo G, Rossi F, Morlando M, Briganti F, Sthandier O, Fatica A, Santini T, Andronache A, Wade M, Laneve P, Rajewsky N, Bozzoni I. Circ-ZNF609 is a circular RNA that can be translated and functions in Myogenesis. Mol Cell. 2017;66(1):22-37 e29.
Yang Y, Fan X, Mao M, Song X, Wu P, Zhang Y, Jin Y, Yang Y, Chen LL, Wang Y, Wong CC, Xiao X, Wang Z. Extensive translation of circular RNAs driven by N(6)-methyladenosine. Cell Res. 2017;27(5):626–41.
Guo R, Chen D, Chen H, Xiong C, Zheng Y, Hou C, Du Y, Geng S, Wang H, Dingding Z, Yilong G. Genome-wide identification of circular RNAs in fungal parasite Nosema ceranae. Curr Microbiol. 2018;75(12):1655–60.
Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587–92.
Benson DA, Karsch-Mizrachi I, Lipman DJ, Ostell J, Wheeler DL. GenBank. Nucleic Acids Res. 2008;36(Database issue):D25-30.
Griffiths-Jones S, Bateman A, Marshall M, Khanna A, Eddy SR. Rfam: an RNA family database. Nucleic Acids Res. 2003;31(1):439–41.
Griffiths-Jones S, Grocock RJ, van Dongen S, Bateman A, Enright AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34(Database issue):D140-144.
Kruger J, Rehmsmeier M. RNAhybrid: microRNA target prediction easy, fast and flexible. Nucleic Acids Res. 2006;34(Web Server issue):W451-454.
Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. 2003;5(1):R1.
Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell. 2005;120(1):15–20.
Chen X, Yang T, Wang W, Xi W, Zhang T, Li Q, Yang A, Wang T. Circular RNAs in immune responses and immune diseases. Theranostics. 2019;9(2):588–607.
Kong C, Yao YX, Bing ZT, Guo BH, Huang L, Huang ZG, Lai YC. Dynamical network analysis reveals key microRNAs in progressive stages of lung cancer. PLoS Comput Biol. 2020;16(5): e1007793.
Ma Z, Li C, Pan G, Li Z, Han B, Xu J, Lan X, Chen J, Yang D, Chen Q, Sang Q, Ji X, Li T, Long M, Zhou Z. Genome-wide transcriptional response of silkworm (Bombyx mori) to infection by the microsporidian Nosema bombycis. PLoS ONE. 2013;8(12): e84137.
Strzyz P. microRNA communication in plants. Nat Rev Mol Cell Biol. 2021;22(12):775.
Sweeney AW, Becnel JJ. Potential of microsporidia for the biological control of mosquitoes. Parasitol Today. 1991;7(8):217–20.
Dakhel WH, Latchininsky AV, Jaronski ST. Efficacy of two Entomopathogenic Fungi, Metarhizium brunneum, strain F52 alone and combined with Paranosema locustae against the migratory grasshopper, Melanoplus sanguinipes, under laboratory and greenhouse conditions. Insects. 2019;10(4):94.
Pu X, Feng M. Efficacy of emulsifiable formulations of two entomopathogenic fungi against small green leafhoppers on tea plant. Ying Yong Sheng Tai Xue Bao. 2004;15(4):619–22.
Szumowski SC, Troemel ER. Microsporidia-host interactions. Curr Opin Microbiol. 2015;26:10–6.
Fernandes JCR, Acuna SM, Aoki JI, Floeter-Winter LM, Muxel SM. Long non-coding RNAs in the regulation of gene expression: physiology and disease. Non-Coding RNA. 2019;5(1):17.
Oberbauer V, Schaefer MR. tRNA-derived small RNAs: biogenesis, modification, function and potential impact on human disease development. Genes. 2018;9(12):607.
Li X, Liu CX, Xue W, Zhang Y, Jiang S, Yin QF, Wei J, Yao RW, Yang L, Chen LL. Coordinated circRNA biogenesis and function with NF90/NF110 in viral infection. Mol Cell. 2017;67(2):214-227 e217.
Huang Q, Evans JD. Identification of microRNA-like small RNAs from fungal parasite Nosema ceranae. J Invertebr Pathol. 2016;133:107–9.
Song Q, Ando A, Jiang N, Ikeda Y, Chen ZJ. Single-cell RNA-seq analysis reveals ploidy-dependent and cell-specific transcriptome changes in Arabidopsis female gametophytes. Genome Biol. 2020;21(1):178.
Darbani B, Noeparvar S, Borg S. Identification of circular RNAs from the parental genes involved in multiple aspects of cellular metabolism in Barley. Front Plant Sci. 2016;7:776.
Ye CY, Chen L, Liu C, Zhu QH, Fan L. Widespread noncoding circular RNAs in plants. New Phytol. 2015;208(1):88–95.
Wang Z, Liu Y, Li D, Li L, Zhang Q, Wang S, Huang H. Identification of circular RNAs in Kiwifruit and their species-specific response to bacterial canker pathogen invasion. Front Plant Sci. 2017;8:413.
Peuvel I, Peyret P, Metenier G, Vivares CP, Delbac F. The microsporidian polar tube: evidence for a third polar tube protein (PTP3) in Encephalitozoon cuniculi. Mol Biochem Parasitol. 2002;122(1):69–80.
KonDo Y, Yoda S, Mizoguchi T, Ando T, Yamaguchi J, Yamamoto K, Banno Y, Fujiwara H. Toll ligand Spatzle3 controls melanization in the stripe pattern formation in caterpillars. Proc Natl Acad Sci USA. 2017;114(31):8336–41.
Cooper EL. Comparative immunology. Integr Zool. 2006;1(1):32–43.
Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta stone of a hidden RNA language? Cell. 2011;146(3):353–8.
Zhou QZ, Zhang B, Yu QY, Zhang Z. BmncRNAdb: a comprehensive database of non-coding RNAs in the silkworm, Bombyx mori. BMC Bioinformatics. 2016;17(1):370.
Zhang C, Zhu B, Dai LS, Liu C, Luo X. Nucleopolyhedroviruses (NPV) induce the expression of small heat shock protein 25.4 in Antheraea pernyi. Gene. 2016;591(2):327–32.
Clerbaux LA, Schultz H, Roman-Holba S, Ruan DF, Yu R, Lamb AM, Bommer GT, Kennell JA. The microRNA miR-33 is a pleiotropic regulator of metabolic and developmental processes in Drosophila melanogaster. Dev Dyn. 2021;250(11):1634–50.
Ghosh J, Bose M, Roy S, Bhattacharyya SN. Leishmania donovani targets Dicer1 to downregulate miR-122, lower serum cholesterol, and facilitate murine liver infection. Cell Host Microbe. 2013;13(3):277–88.
Masoudi F, Sharifi MR, Pourfarzam M. Investigation of the relationship between miR-33a, miR-122, erythrocyte membrane fatty acids profile, and serum lipids with components of metabolic syndrome in type 2 diabetic patients. Res Pharma Sci. 2022;17(3):242–51.
Lampe L, Jentzsch M, Kierszniowska S, Levashina EA. Metabolic balancing by miR-276 shapes the mosquito reproductive cycle and Plasmodium falciparum development. Nat Commun. 2019;10(1):5634.
Shakeel M, Xu X, Xu J, Li S, Yu J, Zhou X, Xu X, Hu Q, Yu X, Jin F. Genome-wide identification of Destruxin A-responsive immunity-related MicroRNAs in diamondback moth Plutella xylostella. Front Immunol. 2018;9:185.
Zhang X, Zheng Y, Jagadeeswaran G, Ren R, Sunkar R, Jiang H. Identification and developmental profiling of conserved and novel microRNAs in Manduca sexta. Insect Biochem Mol Biol. 2012;42(6):381–95.
This work was supported by the National Natural Science Foundation of China (Grant No. 32272942 and 31770159) and the Natural Science Foundation of Chongqing, China (cstc2021jcyj-cxttX0005).
Ethics approval and consent to participate
Silkworms is an unprotected species and it was used as experimental material. Ethic approval is not applicable.
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Fig. S1. Validation of the microsporidia spores by H&E staining and TEM in N. bombycis congenitally infected silkworm embryos and larvae. 5-day embryos (A), 1-day larvae (B), 5-day larvae (C), and 10-day larvae (D) from the parent silkworms infected with N. bombycis were used to detect the microsporidia spores by H&E staining; 5-day embryos (E), 1-day larvae (F), 5-day larvae (G), and 10-day larvae (H) from the parent silkworms infected with N. bombycis were used to detect the microsporidia spores by TEM. Black arrow: yolks; Yellow arrow: muscle cells; Red arrow: adipocyte cells; Blue arrow: epidermic cell. Green arrow: spores. Fig. S2. Principal component analysis (PCA) was performed to evaluate a separation between different stages and study groups. A: PCA was performed based on mRNA data; B: PCA was performed based on lncRNA data; C: PCA was performed based on circRNA data; D: PCA was performed based on miRNA data. Fig. S3. KEGG enrichment analysis of significant differential microsporidia genes between 5-day embryos, 1-day larvae, 5-day larvae, and 10-day larvae challenged by N. bombycis congenital infection. A: KEGG enrichment analysis of cluster I microsporidia genes; B: KEGG enrichment analysis of cluster II microsporidia genes; C: KEGG enrichment analysis of cluster III microsporidia genes; D: KEGG enrichment analysis of cluster IV microsporidia genes; E: KEGG enrichment analysis of cluster V microsporidia genes. Fig. S4. KEGG enrichment analysis of significant differential microsporidia miRNA target genes between 5-day embryos, 1-day larvae, 5-day larvae, and 10-day larvae challenged by N. bombycis congenital infection. A: KEGG enrichment analysis of cluster I miRNA target genes; B: KEGG enrichment analysis of cluster II miRNA target genes; C: KEGG enrichment analysis of cluster III miRNA target genes; D: KEGG enrichment analysis of cluster IV miRNA target genes. Fig. S5 KEGG enrichment analysis of significant differential B. mori genes in silkworm embryos and larvae during the N. bombycis congenital infection. A: KEGG enrichment analysis of significant differential B. mori mRNAs in silkworm 5-day embryos exposed to N. bombycis congenital challenge; B: KEGG enrichment analysis of significant differential B. mori mRNAs in silkworm 1-day larvae exposed to N. bombycis congenital challenge; C: KEGG enrichment analysis of significant differential B. mori mRNAs in silkworm 5-day larvae exposed to N. bombycis congenital challenge; D: KEGG enrichment analysis of significant differential B. mori mRNAs in silkworm 10-day larvae exposed to N. bombycis congenital challenge. Fig. S6. KEGG enrichment analysis of DE lncRNA cis-target genes in silkworm embryos and larvae during the N. bombycis congenital infection. Fig. S7. KEGG enrichment analysis of DE circRNA source genes in silkworm embryos and larvae during the N. bombycis congenital infection. A: KEGG enrichment analysis of DE circRNA source genes in silkworm 5-day embryos during the N. bombycis congenital infection; B: KEGG enrichment analysis of DE circRNA source genes in silkworm 1-day larvae during the N. bombycis congenital infection; C: KEGG enrichment analysis of DE circRNA source genes in silkworm 5-day larvae during the N. bombycis congenital infection; D: KEGG enrichment analysis of DE circRNA source genes in silkworm 10-day larvae during the N. bombycis congenital infection. Fig. S8. KEGG enrichment analysis of DE miRNA target genes in silkworm embryos and larvae during the N. bombycis congenital infection. A: KEGG enrichment analysis of DE miRNA target genes in silkworm 5-day embryos during the N. bombycis congenital infection; B: KEGG enrichment analysis of DE miRNA target genes in silkworm 1-day larvae during the N. bombycis congenital infection; C: KEGG enrichment analysis of DE miRNA target genes in silkworm 5-day larvae during the N. bombycis congenital infection; D: KEGG enrichment analysis of DE miRNA target genes in silkworm 10-day larvae during the N. bombycis congenital infection.
Table S1. Novel genes of N. bombycis were identified in the microsporidia congenital infection in silkworm embryos and larvae. Table S2. Stage-specifically expressed genes of N. bombycis in the microsporidia congenital infection in silkworm embryos and larvae. Table S3. Differentially expressed mRNA, lncRNA, and miRNA of N. bombycis between four stages. Table S4. LncRNA-mRNA analysis (lncRNA antisense, cis and trans-regulation) of N. bombycis. Table S5. N.bombycis miRNA-target genes analysis. Table S6. LncRNA-miRNA-mRNA ceRNA networks of N. bombycis in the microsporidia congenitally infected silkworm embryos and larvae. Table S7. CircRNA -miRNA-mRNA ceRNA networks of N. bombycis in the microsporidia congenitally infected silkworm embryos and larvae. Table S8. The connection degree of each gene of N. bombycis in the ceRNA networks. Table S9. Novel genes of silkworm were identified in the microsporidia congenital infection in silkworm embryos and larvae. Table S10. The lncRNA-mRNA pairs including 99 silkworm lncRNAs and 105 mRNAs were predicted in the microsporidia congenitally infected silkworm embryos and larvae. Table S11. LncRNA-miRNA-mRNA networks of silkworm in the 5th day of N. bombycis congenitally infected silkworm embryos. Table S12. CircRNA-miRNA-mRNA networks of silkworm in the 5th day of N. bombycis congenitally infected silkworm embryos. Table S13. LncRNA-miRNA-mRNA networks of silkworm in the 1th day of N. bombycis congenitally infected silkworm larvae. Table S14. CircRNA-miRNA-mRNA networks of silkworm in the 1th day of N. bombycis congenitally infected silkworm larvae. Table S15. LncRNA-miRNA-mRNA networks of silkworm in the 5th day of N. bombycis congenitally infected silkworm larvae. Table S16. CircRNA-miRNA-mRNA networks of silkworm in the 5th day of N. bombycis congenitally infected silkworm larvae. Table S17. LncRNA-miRNA-mRNA networks of silkworm in the 10th day of N. bombycis congenitally infected silkworm larvae. Table S18. CircRNA-miRNA-mRNA networks of silkworm in the 10th day of N. bombycis congenitally infected silkworm larvae. Table S19. The connection degree of each gene of silkworm in the lncRNA ceRNA network during the N. bombycis congenital infection. Table S20. The connection degree of each gene of silkworm in the circRNA ceRNA network during the N. bombycis congenital infection.
About this article
Cite this article
Shen, Z., Yang, Q., Luo, L. et al. Non-coding RNAs identification and regulatory networks in pathogen-host interaction in the microsporidia congenital infection. BMC Genomics 24, 420 (2023). https://doi.org/10.1186/s12864-023-09490-3