Next Article in Journal
Identification of Genes Involved in Resistance to High Exogenous 20-Hydroxyecdysone in Spodoptera litura
Previous Article in Journal
The Development and Evaluation of Insect Traps for the Asian Citrus Psyllid, Diaphorina citri (Hemiptera: Psyllidae), Vector of Citrus Huanglongbing
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Transcriptome Sequencing Highlights the Regulatory Role of DNA Methylation in Immune-Related Genes’ Expression of Chinese Oak Silkworm, Antheraea pernyi

1
State Key Laboratory of Silkworm Genome Biology, Key Laboratory of Sericultural Biology and Genetic Breeding, Ministry of Agriculture, Southwest University, Chongqing 400716, China
2
Cancer Center, Medical Research Institute, Southwest University, Chongqing 400716, China
*
Authors to whom correspondence should be addressed.
These authors contributed equally to this work.
Insects 2022, 13(3), 296; https://doi.org/10.3390/insects13030296
Submission received: 18 January 2022 / Revised: 10 March 2022 / Accepted: 10 March 2022 / Published: 17 March 2022

Abstract

:

Simple Summary

Antheraea pernyi is an important lepidopteran that has been used as a model insect for various physiological systems. DNA methylation has recently been discovered to control a variety of physiological processes in animals throughout their lives. We investigated the effect of DNA methylation on the innate immune system of A. pernyi and found that DNA methylation inhibition modifies gene expression, especially in immune-related genes. Finally, our data suggest that DNA methylation may regulate gene expression in insects, and microbial infection likely induces DNA methylation of the host genome.

Abstract

Antheraea pernyi is an important lepidopteran used as a model insect species to investigate immune responses, development, and metabolism modulation. DNA methylation has recently been found to control various physiological processes throughout the life of animals; however, DNA methylation and its effect on the physiology of insects have been poorly investigated so far. In the present study, to better understand DNA methylation and its biological role in the immune system, we analyzed transcriptome profiles of A. pernyi pupae following DNA methylation inhibitor injection and Gram-positive bacteria stimulation. We then compared the profiles with a control group. We identified a total of 55,131 unigenes from the RNA sequence data. A comparison of unigene expression profiles showed that a total of 680 were up-regulated and 631 unigenes were down-regulated in the DNA-methylation-inhibition-bacteria-infected group compared to the control group (only bacteria-injected pupae), respectively. Here, we focused on the immune-related differentially expressed genes (DEGs) and screened 10 genes that contribute to immune responses with an up-regulation trend, suggesting that microbial pathogens evade host immunity by increasing DNA methylation of the host genome. Furthermore, several other unigenes related to other pathways were also changed, as shown in the KEGG analysis. Taken together, our data revealed that DNA methylation seems to play a crucial biological role in the regulation of gene expression in insects, and that infection may enhance the host genome DNA methylation by a yet-unknown mechanism.

1. Introduction

DNA methylation is an epigenetic molecular mechanism that modulates various physiological processes in mammals and plants, such as gene regulation, DNA repair, development, tumorigenesis, and nutrigenomics [1,2,3]. Recent studies have also shown that DNA methylation has a crucial biological role in host immune responses. For example, viruses, especially DNA tumor viruses, can stimulate aberrant DNA methylation of the host genome to reduce the host immune responses and thereby evade antiviral immunity [4,5,6]. Viruses also affect host DNA methyltransferase expression patterns, the DNA methylation regulators for epigenetic dysregulation of immune-associated gene expression in host cells [7,8]. Interestingly, it has been shown that demethylation treatment of cancer cells can activate the viral RNA transcription from dormant endogenous retroviruses and trigger antiviral signaling [9,10]. In comparison to mammals, DNA methylation and its effect on different physiological activities of insects have been poorly investigated. Recent progress in DNA methylation research methods has prompted researchers to analyze the role of DNA methylation in insects. Whole-genome bisulfite sequencing studies on different insect species, e.g., Bombyx mori, Tribolium castaneum, Drosophila melanogaster, and Nasonia vitripennis [11,12,13,14], show that DNA methylation exists in the insect genome [15,16,17]. DNA methylation has been linked to immunity, aging, evolution, memory, and caste determination modulation in bees and other insect species [18,19,20,21,22,23,24].
The biological role of DNA methylation in innate immunity in insects has not been explored well. In the Aedes aegypti, genome-wide patterns of DNA methylation were disrupted after Wolbachia infection [25]. Similarly, Metarhizium anisopliae infection caused the differential expression of DNA methyltransferase (DNMT) genes in the larvae of Galleria mellonella [17]. Bombyx mori cytoplasmic polyhedrosis virus (BmCPV) infection in Bombyx mori induced the 27 gene differential expression and differential methylation in the midgut and fat body of infected larvae, respectively, indicating that BmCPV infection modifies the gene expression by mediating variations in DNA methylation [26]. These results suggest that DNA methylation may play a crucial biological role in the immune response of insects against pathogens. However, besides the importance of DNA methylation, this mechanism has poorly been studied, underscoring the indispensability of DNA methylation in insect physiology.
The Chinese oak silkworm, A. pernyi, has been used as a model insect species to study immune response, development, breeding technique, and metabolism regulation. In addition, it is cultured in various Asian countries, including China, Korea, and India, for silk production and as a highly nutritious food [27,28]. Thus, investigating different molecular mechanisms is important to understand the physiological responses of this species and that of other insects. In the present study, we construct a transcriptome-sequencing library from the fat bodies of A. pernyi after administrating DNA methylation inhibitors and a bacterial challenge and then compare it with the control group. The results are interesting, indicating the presence of DNA methylation in A. pernyi that regulates gene differential expression, especially of immune-related genes.

2. Materials and Methods

2.1. Experimental Insects

The Sericulture Research Institute in Henan, China, provided A. pernyi larvae, which were then fed fresh oak leaves and kept at the State Key Laboratory of Silkworm Genome Biology, Southwest University, Chongqing, China. They were kept at room temperature (21–26 °C), with a light/dark photoperiod of 10/14 h and 70% relative humidity. These larvae were allowed to transform into pupae, which were maintained at room temperature.

2.2. 5-Azacytidine Administration, Bacterial Challenge, and Tissue Collection

We divided the pupae into two groups of nine. 5-azacytidine (5-AZA), a DNA methylation inhibitor, and Gram-positive bacteria (Bacillus cereus) were injected into first group, while only Gram-positive bacteria were injected into the second. The inhibition of DNMTs was also investigated by using qRT-PCR. Although a direct effect of 5-AZA on the bacteria used in this study cannot be ruled out, previous studies showed that 5-AZA has no impact on bacterial growth at the concentration used in the present study [23]. For a bacterial challenge, we injected the bacterial cells (4 × 103 cells/larva detected by the method of the colony-forming unit) and 1 μL 5-AZA (40 mM) into pupae after the fourth day since the larvae transformed to pupae. Injections were performed between the second and third tergite of the abdomen of the pupae. Bacterial infection was completed 24 h post second 5-AZA administration. The fat bodies of A. pernyi pupae were collected after 12 h of bacterial infection and then immediately stored in a refrigerator at −80 °C. Three pupal fat bodies were mixed together and taken as one sample, which was then repeated three times.

2.3. RNA Preparation, Library Construction, and Sequencing

Total RNA was extracted from pupal fat bodies using Trizol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The purity of the RNA was analyzed using a Nano spectrophotometer (Implen, Westlake Village, CA, USA), and the concentration was measured using the Qubit RNA Assay Kit in a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). A total of six RNA libraries were constructed using Illumina TruSeq RNA preparation kits, which were assembled in accordance with the instructions provided by the manufacturer. Further verification and quantification were carried out using the Qubit dsDNA BR Assay Kit (Q32850, Invitrogen™, Carlsbad, CA, USA) and the Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). To obtain raw sequencing data in this study, the HiSeq™ 2000 Sequencing System (Illumina, San Diego, CA, USA) was employed.

2.4. Transcriptome Assembly, Annotation, and Function Enrichment

The obtained raw sequencing data were processed to generate clean reads by removing reads containing adapters, poly-N, and low-quality reads using Cutadapt and Perl scripts in-house. The sequence quality, such as Q20, Q30, and GC-content of the clean reads, was confirmed by FastQC. All downstream analyses were based on clean data with high quality. Trinity software was used for de novo assembly of transcriptome data [28], and finally, unigene was obtained. The Salmon method [29] was utilized to execute the expression level for unigenes by calculating TPM [30]. The differentially expressed unigenes were selected with log2 (fold change) > 1 or log2 (fold change) < −1 and with statistical significance (p-value < 0.05) by the R package edger.

2.5. Unigene Annotation and Functional Classification

For annotation and functional classification, assembled unigenes were aligned against the non-redundant (Nr) protein database (http://www.ncbi.nlm.nih.gov/, accessed on 14 May 2021), Gene ontology (GO) (http://www.geneontology.org, accessed on 14 May 2021), Clusters of Orthologous Groups of proteins (KOG/COG), SwissProt (http://www.expasy.ch/sprot/, accessed on 14 May 2021), Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/, accessed on 14 May 2021), and eggNOG (http://eggnogdb.embl.de/, accessed on 14 May 2021) databases using DIAMOND with a threshold of E-value < 0.00001. GO and KEGG enrichment analysis was again performed on the differentially expressed unigenes by Perl scripts in-house.

2.6. Quantitative RT-PCR Analysis

Ten up-regulated DEGs were selected for quantitative RT-PCR to validate the transcriptome data. Total RNA was isolated from the 5-AZA+bacteria-injected samples and the control samples using Trizol solution according to the manufacturer’s instructions for cDNA synthesis, which was then used as a template for the qRT-PCR assay for the validation of RNA-Seq data. To eliminate any DNA contamination, all of the samples were treated with DNase I (Promega). The RNA concentrations were determined using a spectrophotometer, and a total of 2 µg of RNA was reverse-transcribed with Hifair® III 1st Strand cDNA Synthesis SuperMix for qPCR kit. The gene-specific primers were designed using the reference sequences, and 18s RNA was used for normalization (Supplementary Table S1). qRT-PCR was performed using a 20 μL reaction system with a procedure as follows: 95 °C for 30 s, followed by 39 cycles of 95 °C for 5 s, and 60 °C for 30 s. Melting curves were generated after each run to ensure a single PCR product. Each reaction was run in triplicate. The mRNA quantity of each gene was calculated with the 2−ΔΔCT method [31].

2.7. Statistical Analysis

In this study, all experiments were executed in triplicate, and the obtained data represented the means ± S.E. The one-way analysis of variance (ANOVA) and Tukey’s multiple range tests were used to evaluate the difference between groups.

3. Results

3.1. Summary of Transcriptome Sequence and Assembly

To investigate the effect of the DNA methylation inhibitor 5-AZA, two groups of pupae (5-AZA-treated and a control group) were sampled to collect fat body tissue after bacterial (B. cereus) treatment. We also examined the suppression of DNMTs prior to bacterial treatment (Figure 1A) and then sequenced the collected samples, with three replicates included in each group. After filtering the data, the transcriptome sequence generated 30,212,355, 24,056,953, and 25,478,892 clean reads from the 5-AZA-Gram-positive-bacteria (B. cereus)-infected group and 24,534,056, 26,678,456, and 23,536,596 from the control group (Table 1). The obtained sequenced data were assembled de novo using paired-end raw reads generated by the Illumina HiSeq2500 instrument. Quality evaluation was carried out in accordance with Illumina guidelines, and it was found that >95% and >92% of the sequencing data for each experimental group (treated with 5-AZA and bacteria) and control group (treated with only bacteria), respectively, had Q20 and Q30 quality scores. The contents for each experimental and control group were 43.95%, 43.00%, 42.21% and 40.95%, 42.52%, and 42.58%, respectively (Table 1). As shown in Table 2, 117,272 transcripts were generated using Trinity software with a mean length of 1425.323 bp and an N50 length of 2357 bp. A total of 55,131 unigenes were generated with an average length of 977.2239 bp and an N50 length of 1621 bp. Among these unigenes, 26,435 (47.95%) were in the range of 300−500 bp, 14,437 (26.19%) were 500−1000 bp, 7620 (13.82%) were 1−2 kbp, and 6638 (12.04%) were longer than 2 kbp. These observations propelled us to conclude that the obtained data are of high quality, and thus, the unigenes can be further used for annotation analysis. There were a total of 27,361 unigenes with 8187, 12,562, 13,351, 17,920, 19,325, 14,013, 25,402, and 23,492 annotated COG, GO, KEGG, Pfam, Swiss-Prot, eggNOG, and NCBI_nr, respectively (Supplementary Table S2).

3.2. Analysis of Differentially Expressed Genes

The expression patterns of genes in the 5AZA + bacteria (B. cereus)-injected group were compared with the control samples. In the comparison, it was found that by applying the standard threshold of p ≤ 0.001, false discovery rate ≤ 0.001, and |log fold change (FC)| ≥ 1, a total of 680 up-regulated and 631 down-regulated DEGs were screened out (Figure 1B and Figure 2).
To understand the functional classification of DEGs, the unigene sequences were aligned against the GO and COG databases. The DEGs were classified into cellular components, molecular functions, and biological processes. The biological process groups possessed 20 subcategories, while cellular component and molecular function groups represented 15 subcategories each. Among the biological processes, cellular processes, metabolic processes, and the single-organism metabolic process, the genes were greatly expressed, whereas the cell, cell part, membrane, and organelle genes were highly expressed in cellular components. Finally, catalytic activity and binding genes were the most abundant among the molecular functions (Figure 3, Supplementary Table S3).
The Eukaryotic Orthologous Groups (KOG) database categorizes homologous gene products, and each KOG protein is presumed to be derived from an ancestral protein. Based on the KOG functional classification, overall unigenes were mapped into 21 KOG categories, with general function prediction only. This was followed by signal-transduction mechanisms and posttranslational modifications, protein turnover, chaperones, and defense mechanisms. This suggests that unigenes may be involved in a variety of functions, especially immune defense in A. pernyi (Figure 1B). For the up-regulated unigenes, the most enriched pathways were phagosome, protein processing in the endoplasmic reticulum, oxidative phosphorylation, and the biosynthesis of amino acids (Figure 4A, Supplementary Table S4). The Fanconi anemia pathway, the lysosome, purine metabolism, and the apoptosis pathways were among the genes that were down-regulated (Figure 4B, Supplementary Table S4).

3.3. Validation and Reliability of the Transcriptome Data by qRT-PCR

First, we searched for the housekeeping gene actin in both the experimental (5-AZA+B. cereus) and control groups to ensure that the data were reliable before proceeding. We found that the expression of this housekeeping gene was the same in both groups. After confirming housekeeping gene expression, we selected 10 DEGs from A. pernyi related to the innate immune system and with different signaling pathways. In order to validate the expression profiles in both the experimental and control groups, we performed qRT-PCR analysis. The results revealed that the expression profiles of the candidate genes were consistent with those of the RNA-seq data (Figure 5), indicating that the RNA sequencing data were reliable. These findings will be useful in improving our understanding of the process of DNA methylation and its effect on gene expression.

4. Discussion

Unlike mammals, where DNA methylation has extensively been reported, many insect species have been shown to be devoid of the DNA-methylation process [22]. When it comes to most insect species, there is currently a dearth of literature relating to the existence of this epigenetic phenomenon. Therefore, our understanding of DNA methylation and its involvement in various physiological responses is extremely limited at this time. In mammals, the biological roles of DNA methylation have been extensively investigated; however, in insects, this information is still in its early stages. To investigate the biological roles of DNA methylation in A. pernyi innate immunity, the differential expression of genes related to immunity was evaluated after applying 5-AZA and bacteria. The suppression of DNMTs was confirmed by qRT-PCR. Interestingly, in addition to DNMT-1, DNMT-2 expression was also significantly reduced. Previous studies on Drosophila melanogaster, Dictyostelium discoideum, and Entamoeba histolytica showed that DNMT-2 could regulate DNA methylation [32,33,34] in addition to tRNA methylation [35]. There are multiple lines of evidence suggesting that 5-AZA, as a ribonucleoside, can also form covalent bonds with human and mouse DNMT-2 homologs and has the ability to inhibit its expression [35,36,37]. For this purpose, transcriptome sequencing was used to compare the genes in the DNA-methylation-inhibitor-bacteria infected with the control. In this study, we focused on the immune-related DEGs and reported that the inhibition of DNA methylation has an effect on the expression of genes in insects.
Insects have a complex immune system that is divided into humoral and cellular immunity components. Cellular immunity protects insects from pathogen infection through encapsulation, phagocytosis, and nodulation, while humoral immune responses include the production of antimicrobial peptides [38,39]. In the present study, from the RNA sequence data, we observed that several DEGs encoding immune-related proteins regulating cellular and humoral immunity were up-regulated (680 unigenes) with variable expression. Among these up-regulated, we selected 10 DEGs (c89459.graph_c0, c92066.graph_c0, c90624.graph_c0, c91406.graph_c2, c93759.graph_c1, c93774.graph_c0, c77163.graph_c1, c85165.graph_c0, c89772.graph_c1, c89772.graph_c1, c85340.graph_c0 encoding mitogen-activated protein kinase kinase kinase 4, integrin beta pat-3-like, GNBP, suppressor of cytokine signaling 2-like isoform X2, peptidoglycan recognition protein-like protein, cactus, spaetzle, cytochrome P450, serine protease inhibitor 13 precursor, hemolymph proteinase 9) representing different signaling pathways to validate the RNA sequence data and to highlight the importance of DNA methylation in the physiological processes of insects. A growing body of evidence suggests that DNA methylation, in particular cytosine (CpG) methylation, is a critical and reversible gene-regulation mechanism. DNA methylation controls transcriptional modulation by influencing the recruitment of regulatory factors to gene promoters [40]. Many studies have reported that microbial infection induces genome DNA methylation and that this is linked with host immune responses [41,42]. For example, Laurson [43] reported that human papillomavirus infection induces aberrant host DNA methylation through direct interaction of the viral protein E7 with the DNA methylation enzyme DNMT-1 [43]. Escherichia coli infection regulates the transcription patterns of DNMT-1, and inducing CpG hypermethylation reduces the expression of the cell-cycle inhibitor CDKN2A in human uroepithelial cells [44]. Similarly, Mycobacterium tuberculosis infection causes fast methylation of host DNA at distal enhancer elements and associated chromatin remodeling [45]. Kausar et al. [22] recently reviewed that, despite the fact that DNA methylation levels in insects are low, DNA methylation has the ability to modify a variety of physiological processes, including the innate immune system. From these epigenetic effects of microbial infections and subsequent modification of the innate immune response, it is clear that DNA methylation plays a role in the interactions between host and pathogen.
Based on our RNA sequencing data, we found that inhibiting DNA methylation could have an impact on the insect host–pathogen interaction. The Toll (GNBP, cactus, serine protease inhibitor protein, hemolymph protease 9, Spatzle), IMD (peptidoglycan recognition protein), and JAK-STAT (suppressor of cytokine signaling 2) signaling pathways are all important in the innate immunity of insects. In response to Gram-negative bacterial infections and viruses, IMD and JAK-STAT are activated, whereas the Toll pathway is activated in response to Gram-positive bacteria and fungi [46,47,48]. Following bacterial invasion, pattern-recognition receptors are sensed by GNBP, which then activates Spatzle, which plays a vital role in the activation of the Toll pathway. The activated Spatzle binds to Toll receptors for the activation of the downstream cascade [48]. Cactus is phosphorylated and cleaved off Dif late in this process, allowing Dif to translocate into the nucleus and modulate the transcription of effector genes [49]. Serine protease inhibitors and hemolymph protease control the melanization process of the hemolymph to attenuate pathogens [39]. The peptidoglycan-recognition protein has a strong affinity for Gram-positive bacteria and can also bind with Gram-negative bacteria to initiate signaling via IMD and the Toll pathway [49], whereas the suppressor of cytokine 2 is a well-known negative regulator of the JAK-STAT pathway, which checks the activity of this pathway during infection and development [50,51,52,53]. These pathways mainly attenuate microbial infection via inducing the production of antimicrobial peptides [39,49]. In the present study, when A. pernyi was exposed to a DNA methylation inhibitor and a Gram-positive bacterial challenge, we observed that genes associated with the Toll, IMD, and JAK-STAT pathways were increased in comparison to the control group. This suggests that DNA methylation inhibition induced robust immune responses in A. pernyi. 5-AZA, a commonly used DNA methylation inhibitor, seems to impair the production of DNA-methylation-associated genes, particularly DNMTs. This impairment inhibits the methylation status at the gene loci, which may affect the expression of genes [23].
Evidence suggests that these three pathways are interlinked and that they are activated in order to reduce microbial infection. For example, according to Li and Xiang [54], the Toll pathway in invertebrates is activated not only by Gram-positive bacteria but also by Gram-negative bacteria and the WSSV virus. In addition, Hedengren-Olcott et al. [55] suggested that the activation pathways seem to be specific to bacterial strains. They went on to say that some bacteria can induce both the Toll and IMD pathways, suggesting that there is some crosstalk between these immune-signaling pathways. Additionally, integrin, a cell surface receptor protein that has been reported to contribute to encapsulation, pathogen adherence, and phagocytosis, was increased in addition to these immune-pathway-related genes [56,57]. Here, we discussed our results very precisely because our data showed that the inhibition of DNA methylation has an impact on a variety of physiological processes, such as development. Altogether, these results suggest that pathogen infection may induce DNA methylation of the host genome, which may suppress immune responses by yet-unknown mechanisms. On the other hand, our data suggest that inhibiting DNA methylation could induce a robust immune response in insects, thereby shaping the host–pathogen interaction.

5. Conclusions

In conclusion, our findings imply that DNA methylation probably occurs in the Chinese oak silkworm, A. pernyi, similar to that of B. mori. We also identified a DNA methylation tool kit in this species (data not shown), which further ensured that the A. pernyi genome is methylated. In addition, in this species, inhibition of DNA methylation might result in the induction of the expression of genes, especially immune-related genes, which can influence the host–pathogen interaction. This study showed that DNA methylation inhibition caused modifications in a number of signaling pathways. Thus, this study provides a repository to investigate the level of DNA methylation in the genome of the insect, and subsequently, which genes promoted are methylated and how DNA methylation modifies different physiological processes, including development and immunity.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/insects13030296/s1, Table S1: Gene specific Primers used for qRT-PCR analysis.; Table S2: Summary Statistics of the A. pernyi Transcriptome Annotation.; Table S3: Table of GO Enrichment of All Genes and DEGs.; Table S4: Classification of the Unigenes Annotated in KOG.

Author Contributions

Conceptualization, investigation, writing—original draft preparation, S.K. and I.G.; software and formal analysis, R.L.; supervision, project administration and funding acquisition, M.N.A. and H.C. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the China Postdoctoral Science Foundation (2021M692676), the Chongqing Postdoctoral Science Foundation (cstc2021jcyj-bsh0128), the Doctoral Start-up Fund of Southwest University (SWU020023), the National Natural Science Foundation of China (31802142), the Fundamental Research Funds for the Central Universities (XDJK2019C089), and the Foundation of State Key Laboratory of Silkworm Genome Biology (SKLSGB-ORP202003).

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in article and Supplementary Material.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Park, J.; Peng, Z.; Zeng, J.; Elango, N.; Park, T.; Wheeler, D.; Werren, J.H.; Yi, S.V. Comparative analyses of DNA methylation and sequence evolution using Nasonia genomes. Mol. Biol. Evol. 2011, 28, 3345–3354. [Google Scholar] [CrossRef] [PubMed]
  2. Baylin, S.B.; Jones, P.A. A decade of exploring the cancer epigenome-biological and translational implications. Nat. Rev. Cancer 2011, 11, 726–734. [Google Scholar] [CrossRef] [PubMed]
  3. Riccio, A.; Alvania, R.S.; Lonze, B.E.; Ramanan, N.; Kim, T.; Huang, Y.; Dawson, T.M.; Snyder, S.H.; Ginty, D.D. A nitric oxide signaling pathway controls CREB-mediated gene expression in neurons. Mol. Cell 2009, 21, 283–294. [Google Scholar] [CrossRef] [PubMed]
  4. Cicchini, L.; Blumhagen, R.Z.; Westrich, J.A.; Myers, M.E.; Warren, C.J.; Siska, C.; Raben, D.; Kechris, K.J.; Pyeon, D. High-risk human papillomavirus E7 alters host DNA Methylome and represses HLA-E expression in human keratinocytes. Sci. Rep. 2017, 7, 3633. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Dongt, S.M.; Lee, H.G.; Cho, S.G.; Kwon, S.H.; Yoon, H.; Kwon, H.J.; Lee, J.H.; Kim, H.; Park, P.G.; Kim, H.; et al. Hypermethylation of the interferon regulatory factor 5 promoter in Epstein-Barr virus-associated gastric carcinoma. J. Microbiol. 2015, 53, 70–76. [Google Scholar] [CrossRef]
  6. Di Bartolo, D.L.; Cannon, M.; Liu, Y.F.; Renne, R.; Chadburn, A.; Boshoff, C.; Cesarman, E. KSHV LANA inhibits TGF-beta signaling through epigenetic silencing of the TGF-beta type II receptor. Blood 2008, 111, 4731–4740. [Google Scholar] [CrossRef]
  7. Kang, M.S.; Kieff, E. Epstein-Barr virus latent genes. Exp. Mol. Med. 2015, 47, e131. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  8. Wu, J.; Xu, Y.; Mo, D.; Huang, P.; Sun, R.; Huang, L.; Pan, S.; Xu, J. Kaposi’s sarcoma-associated herpesvirus (KSHV) vIL-6 promotes cell proliferation and migration by upregulating DNMT1 via STAT3 activation. PLoS ONE 2014, 9, e93478. [Google Scholar] [CrossRef]
  9. Roulois, D.; Loo Yau, H.; Singhania, R.; Wang, Y.; Danesh, A.; Shen, S.Y.; Han, H.; Liang, G.; Jones, P.A.; Pugh, T.J.; et al. DNA-Demethylating agents target colorectal Cancer cells by inducing viral mimicry by endogenous transcripts. Cell 2015, 162, 961–973. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Chiappinelli, K.B.; Strissel, P.L.; Desrichard, A.; Li, H.; Henke, C.; Akman, B.; Hein, A.; Rote, N.S.; Cope, L.M.; Snyder, A.; et al. Inhibiting DNA methylation causes an interferon response in Cancer via dsRNA including endogenous retroviruses. Cell 2015, 162, 974–986. [Google Scholar] [CrossRef] [Green Version]
  11. Capuano, F.; Mulleder, M.; Kok, R.; Blom, H.J.; Ralser, M. Cytosine DNA methylation is found in Drosophila melanogaster but absent in Saccharomyces cerevisiae, Schizosaccharomyces pombe, and other yeast species. Anal. Chem. 2014, 86, 3697–3702. [Google Scholar] [CrossRef]
  12. Xiang, H.; Zhu, J.; Chen, Q.; Dai, F.; Li, X.; Li, M.; Zhang, H.; Zhang, G.; Li, D.; Dong, Y.; et al. Single base-resolution methylome of the silkworm reveals a sparse epigenomic map. Nat. Biotechnol. 2010, 28, 516–520. [Google Scholar] [CrossRef] [PubMed]
  13. Feliciello, I.; Parazajder, J.; Akrap, I.; Ugarkovic, D. First evidence of DNA methylation in insect Tribolium castaneum: Environmental regulation of DNA methylation within heterochromatin. Epigenetics 2013, 8, 534–541. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Beeler, S.M.; Wong, G.T.; Zheng, J.M.; Bush, E.C.; Remnant, E.J.; Oldroyd, B.P.; Drewell, R.A. Whole-genome DNA methylation profile of the jewel wasp (Nasonia vitripennis). G3 (Bethesda) 2014, 4, 383–388. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Glastad, K.M.; Hunt, B.G.; Yi, S.V.; Goodisman, M.A. DNA methylation in insects: On the brink of the epigenomic era. Insect Mol. Biol. 2011, 20, 553–565. [Google Scholar] [CrossRef]
  16. Standage, D.S.; Berens, A.J.; Glastad, K.M.; Severin, A.J.; Brendel, V.P.; Toth, A.L. Genome, transcriptome and methylome sequencing of a primitively eusocial wasp reveal a greatly reduced DNA methylation system in a social insect. Mol. Ecol. 2016, 25, 1769–1784. [Google Scholar] [CrossRef] [PubMed]
  17. Wu, P.; Jie, W.; Shang, Q.; Annan, E.; Jiang, X.; Hou, C.; Chen, T.; Guo, X. DNA methylation in silkworm genome may provide insights into epigenetic regulation of response to Bombyx mori cypovirus infection. Sci. Rep. 2017, 7, 16013. [Google Scholar] [CrossRef] [Green Version]
  18. Bonasio, R.; Li, Q.; Lian, J.; Mutti, N.S.; Jin, L.; Zhao, H.; Zhang, P.; Wen, P.; Xiang, H.; Ding, Y.; et al. Genome-wide and caste-specific DNA methylomes of the ants Camponotus floridanus and Harpegnathos saltator. Curr. Biol. 2012, 22, 1755–1764. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  19. Yan, H.; Bonasio, R.; Simola, D.F.; Liebig, J.; Berger, S.L.; Reinberg, D. DNA methylation in social insects: How epigenetics can control behavior and longevity. Annu. Rev. Entomol. 2015, 60, 435–452. [Google Scholar] [CrossRef] [PubMed]
  20. Biergans, S.D.; Giovanni, G.C.; Judith, R.; Charles, C. Dnmts and Tettarget memory-associated genes after appetitive olfactory training in honey bees. Sci. Rep. 2015, 5, 21656. [Google Scholar] [CrossRef] [Green Version]
  21. Herb, B.R.; Wolschin, F.; Hansen, K.D.; Aryee, M.J.; Langmead, B.; Irizarry, R.; Amdam, G.V.; Feinberg, A.P. Reversible switching between epigenetic states in honeybee behavioral subcastes. Nat. Neurosci. 2012, 15, 1371–1373. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Kausar, S.; Abbas, M.N.; Cui, H. A review on the DNA methyltransferase family of insects: Aspect and prospects. Int. J. Biol. Macromol. 2021, 186, 289–302. [Google Scholar] [CrossRef] [PubMed]
  23. Baradaran, E.; Moharramipour, S.; Asgari, S.; Mehrabadi, M. Induction of DNA methyltransferase genes in Helicoverpa armigera following injection of pathogenic bacteria modulates expression of antimicrobial peptides and affects bacterial proliferation. J. Insect Physiol. 2019, 118, 103939. [Google Scholar] [CrossRef]
  24. Heitmueller, M.; Billion, A.; Dobrindt, U.; Vilcinskas, A.; Mukherjee, K. Epigenetic mechanisms regulate innate immunity against uropathogenic and commensal-like Escherichia coli in the surrogate insect model Galleria mellonella. Infect. Immun. 2017, 85, e00336-17. [Google Scholar] [CrossRef] [Green Version]
  25. Ye, Y.H.; Woolfit, M.; Huttley, G.A.; Rances, E.; Caragata, E.P.; Popovici, J.; O’Neill, S.L.; McGraw, E.A. Infection with a virulent strain of Wolbachia disrupts genome wide-patterns of cytosine methylation in the mosquito Aedes aegypti. PLoS ONE 2013, 8, e66482. [Google Scholar] [CrossRef] [PubMed]
  26. Vilcinkas, A. The role of epigenetics in host–parasite coevolution: Lessons from the model host insects galleria mellonella and Tribolium castaneum. Zoology 2016, 119, 273–280. [Google Scholar] [CrossRef] [PubMed]
  27. Xin, Z.Z.; Liu, Q.N.; Liu, Y.; Zhang, D.Z.; Wang, Z.F.; Zhang, H.B.; Ge, B.M.; Zhou, C.L.; Chai, X.Y.; Tang, B.P. Transcriptome-Wide Identification of Differentially Expressed Genes in Chinese Oak Silkworm Antheraea pernyi in Response to Lead Challenge. J. Agric. Food Chem. 2017, 65, 9305–9314. [Google Scholar] [CrossRef]
  28. Grabherr, M.G.; Haas, B.J.; Yassour, M.; Levin, J.Z.; Thompson, D.A.; Amit, I.; Adiconis, X.; Fan, L.; Raychowdhury, R.; Zeng, Q.; et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 2011, 29, 644–652. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Patro, R.; Duggal, G.; Love, M.I.; Irizarry, R.A.; Kingsford, C. Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods 2017, 14, 417–419. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  30. Mortazavi, A.; Williams, B.A.; McCue, K.; Schaeffer, L.; Wold, B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 2008, 5, 621–628. [Google Scholar] [CrossRef]
  31. Livak, K.J.; Schmittgen, T.D. Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−ΔΔCT Method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef] [PubMed]
  32. Kunert, N.; Marhold, J.; Stanke, J.; Stach, D.; Lyko, F. A Dnmt2-like protein mediates DNA methylation in Drosophila. Development 2003, 130, 5083–5090. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Fisher, O.; Siman-Tov, R.; Ankri, S. Characterization of cytosine methylated regions and 5-cytosine DNA methyltransferase (Ehmeth) in the protozoan parasite Entamoeba histolytica. Nucleic Acids Res. 2004, 32, 287–297. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Katoh, M.; Curk, T.; Xu, Q.; Zupan, B.; Kuspa, A.; Shaulsky, G. Developmentally regulated DNA methylation in Dictyostelium discoideum. Eukaryot. Cell 2006, 5, 18–25. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Liu, K.; Wang, Y.F.; Cantemir, C.; Muller, M.T. Endogenous assays of DNA methyltransferases: Evidence for differential activities of DNMT1, DNMT2, and DNMT3 in mammalian cells In vivo. Mol. Cell. Biol. 2003, 23, 2709–2719. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Schaefer, M.; Hagemann, S.; Hanna, K.; Lyko, F. Azacytidine inhibits RNA methylation at DNMT2 target sites in human cancer cell lines. Cancer Res. 2009, 69, 8127–8132. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Geyer, K.K.; Rodríguez López, C.M.; Chalmers, I.W.; Munshi, S.E.; Truscott, M.; Heald, J.; Wilkinson, M.J.; Hoffmann, K.F. Cytosine methylation regulates oviposition in the pathogenic blood fluke Schistosoma mansoni Kathrin. Nat. Commun. 2011, 2, 424. [Google Scholar] [CrossRef] [Green Version]
  38. Walderdorff, L.; Laval-Gilly, P.; Bonnefoy, A.; Falla-Angel, J. Imidacloprid intensifies its impact on honeybee and bumblebee cellular immune response when challenged with LPS (lippopolysacharide) of Escherichia coli. J. Insect Physiol. 2018, 108, 17–24. [Google Scholar] [CrossRef]
  39. Kausar, S.; Abbas, M.N.; Zhao, Y.; Cui, H. Immune strategies of silkworm, Bombyx mori against microbial infections. Invertebr. Surviv. J. 2019, 16, 130–140. [Google Scholar]
  40. Jones, P.A. Functions of DNA methylation: Islands, start sites, gene bodies and beyond. Nat. Rev. Genet. 2012, 13, 484–492. [Google Scholar] [CrossRef]
  41. De Monerri, N.C.S.; Kim, K. Pathogens hijack the epigenome: A new twist on host-pathogen interactions. Am. J. Pathol. 2014, 184, 897–911. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. McGrath-Morrow, S.A.; Ndeh, R.; Helmin, K.A.; Chen, S.Y.; Anekalla, K.R.; Abdala-Valencia, H.; D’Alessio, F.R.; Collaco, J.M.; Singer, B.D. DNA methylation regulates the neonatal CD4+ T-cell response to pneumonia in mice. J. Biol. Chem. 2018, 293, 11772–11783. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Laurson, J.; Khan, S.; Chung, R.; Cross, K.; Raj, K. Epigenetic repression of E-cadherin by human papillomavirus 16 E7 protein. Carcinogenesis 2010, 31, 918–926. [Google Scholar] [CrossRef] [Green Version]
  44. Tolg, C.; Sabha, N.; Cortese, R.; Panchal, T.; Ahsan, A.; Soliman, A.; Aitken, K.J.; Petronis, A.; Bägli, D.J. Uropathogenic, E. coli infection provokes epigenetic downregulation of CDKN2A (p16INK4A) in uroepithelial cells. Lab Investig. 2011, 91, 825–836. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Pacis, A.; Tailleux, L.; Morin, A.M.; Lambourne, J.; MacIsaac, J.L.; Yotova, V.; Dumaine, A.; Danckaert, A.; Luca, F.; Grenier, J.C.; et al. Bacterial infection remodels the DNA methylation landscape of human dendritic cells. Genome Res. 2015, 25, 1801–1811. [Google Scholar] [CrossRef] [Green Version]
  46. Lemaitre, B.; Nicolas, E.; Michaut, L.; Reichhart, J.M.; Hoffmann, J.A. The dorsoventral regulatory gene cassette spätzle/Toll/cactus controls the potent antifungal response in Drosophila adults. Cell 1996, 86, 973–983. [Google Scholar] [CrossRef] [Green Version]
  47. Kajla, M.; Choudhury, T.P.; Kakani, P.; Gupta, K.; Dhawan, R.; Gupta, L.; Kumar, S. Silencing of Anopheles stephensi heme peroxidase HPX15 activates diverse immune pathways to regulate the growth of midgut bacteria. Front. Microbiol. 2016, 7, 1351. [Google Scholar] [CrossRef] [PubMed]
  48. Gupta, K.; Dhawan, R.; Kajla, M.; Misra, T.; Kumar, S.; Gupta, L. The evolutionary divergence of STAT transcription factor in different Anopheles species. Gene 2017, 596, 89–97. [Google Scholar] [CrossRef]
  49. Yu, B.; Sang, Q.; Pan, G.; Li, C.; Zhou, Z. A Toll-Spätzle Pathway in the immune response of Bombyx mori. Insects 2020, 11, 586. [Google Scholar] [CrossRef]
  50. Kingsolver, M.B.; Huang, Z.; Hardy, R.W. Insect antiviral innate immunity: Pathways, effectors, and connections. J. Mol. Biol. 2013, 425, 4921–4936. [Google Scholar] [CrossRef] [Green Version]
  51. Hoffmann, J.A. The immune response of Drosophila. Nature 2003, 426, 33–38. [Google Scholar] [CrossRef]
  52. Abbas, M.N.; Kausar, S.; Zhao, E.; Cui, H. Suppressors of cytokine signaling proteins as modulators of development and innate immunity of insects. Dev. Comp. Immunol. 2020, 104, 103561. [Google Scholar] [CrossRef]
  53. Abbas, M.N.; Kausar, S.; Gul, I.; Ke, X.X.; Dong, Z.; Lu, X.; Cui, H. Suppressor of cytokine signalling 6 is a potential regulator of antimicrobial peptides in the Chinese oak silkworm, Antheraea pernyi. Mol. Immunol. 2021, 140, 12–21. [Google Scholar] [CrossRef] [PubMed]
  54. Li, F.; Xiang, J. Signaling pathways regulating innate immune responses in shrimp. Fish Shellfish. Immunol. 2013, 34, 973–980. [Google Scholar] [CrossRef] [PubMed]
  55. Hedengren-Olcott, M.; Olcott, M.C.; Mooney, D.T.; Ekengren, S.; Geller, B.L.; Taylor, B.J. Differential activation of the NF-kappaB-like factors Relish and Dif in Drosophila melanogaster by fungi and Gram-positive bacteria. J. Biol. Chem. 2004, 279, 21121–21127. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Lavine, M.D.; Strand, M.R. Haemocytes from Pseudoplusia includes express multiple alpha and beta integrin subunits Insect. Mol. Biol. 2003, 12, 441–452. [Google Scholar]
  57. Li, C.; Zhang, K.; Pan, G.; Zhang, L.; Hu, X.; Zhao, G.; Deng, C.; Tan, M.; Li, C.; Xu, M.; et al. Bmintegrin β1: A broadly expressed molecule modulates the innate immune response of Bombyx mori. Dev. Comp. Immunol. 2021, 114, 103869. [Google Scholar] [CrossRef]
Figure 1. The inhibition of DNMTs and the comparison of the number of differentially expressed genes. (A): The inhibition of DNMTs analyzed by the qRT-PCR. (B): The levels of gene expression in 5-AZA- and bacterial-challenged groups compared to those in the control group. The differentially expressed genes were screened out by |log fold change (FC)| ≥ 1. Bars show the mean ± S.E. (n = 3), and asterisks indicate significant differences (* p < 0.05).
Figure 1. The inhibition of DNMTs and the comparison of the number of differentially expressed genes. (A): The inhibition of DNMTs analyzed by the qRT-PCR. (B): The levels of gene expression in 5-AZA- and bacterial-challenged groups compared to those in the control group. The differentially expressed genes were screened out by |log fold change (FC)| ≥ 1. Bars show the mean ± S.E. (n = 3), and asterisks indicate significant differences (* p < 0.05).
Insects 13 00296 g001
Figure 2. Volcano plot of the differentially expressed genes. The fold change of gene expression between the 5-AZA + bacterial treated group and the control group is represented by the horizontal ordinate. The vertical ordinate exhibits the statistical significance of the change in gene expression. Each gene is represented by a point on the plot, with the blue and red points representing the genes that were significantly down-regulated and up-regulated, respectively.
Figure 2. Volcano plot of the differentially expressed genes. The fold change of gene expression between the 5-AZA + bacterial treated group and the control group is represented by the horizontal ordinate. The vertical ordinate exhibits the statistical significance of the change in gene expression. Each gene is represented by a point on the plot, with the blue and red points representing the genes that were significantly down-regulated and up-regulated, respectively.
Insects 13 00296 g002
Figure 3. The top 25 GO terms in biological process cellular components, and molecular function categories with enriched in the differentially expressed gene. The Y-axis corresponds to the number of DEGs, whereas the X-axis shows different gene functions.
Figure 3. The top 25 GO terms in biological process cellular components, and molecular function categories with enriched in the differentially expressed gene. The Y-axis corresponds to the number of DEGs, whereas the X-axis shows different gene functions.
Insects 13 00296 g003
Figure 4. (A): The top 21 enriched pathways for the up-regulated genes. (B): The top 21 enriched pathways for the down-regulated genes.
Figure 4. (A): The top 21 enriched pathways for the up-regulated genes. (B): The top 21 enriched pathways for the down-regulated genes.
Insects 13 00296 g004
Figure 5. The expression of selected differentially expressed genes for the validation of data. (A): Expression of 10 selected unigenes expression (FPKM value) in the experimental and control group. (B): Quantitative RT-PCR verification of RNA sequence data. X-axis and Y-axis show the 10 selected differentially expressed genes and the relative expression (B. cereus-injected vs. 5-AZA+B. cereus-injected), respectively. The sequences include c89459.graph_c0, c92066.graph_c0, c90624.graph_c0, c91406.graph_c2, c93759.graph_c1, c93774.graph_c0, c77163.graph_c1, c85165.graph_c0, c89772.graph_c1, c85340.graph_c0 encoding mitogen-activated protein kinase kinase 4, integrin beta pat-3-like, GNBP, suppressor of cytokine signaling 2-like isoform X2, peptidoglycan recognition protein-like protein, cactus, spaetzle, cytochrome P450, serine protease inhibitor 13 precursor, and hemolymph proteinase 9. (* p < 0.05, ** p < 0.01).
Figure 5. The expression of selected differentially expressed genes for the validation of data. (A): Expression of 10 selected unigenes expression (FPKM value) in the experimental and control group. (B): Quantitative RT-PCR verification of RNA sequence data. X-axis and Y-axis show the 10 selected differentially expressed genes and the relative expression (B. cereus-injected vs. 5-AZA+B. cereus-injected), respectively. The sequences include c89459.graph_c0, c92066.graph_c0, c90624.graph_c0, c91406.graph_c2, c93759.graph_c1, c93774.graph_c0, c77163.graph_c1, c85165.graph_c0, c89772.graph_c1, c85340.graph_c0 encoding mitogen-activated protein kinase kinase 4, integrin beta pat-3-like, GNBP, suppressor of cytokine signaling 2-like isoform X2, peptidoglycan recognition protein-like protein, cactus, spaetzle, cytochrome P450, serine protease inhibitor 13 precursor, and hemolymph proteinase 9. (* p < 0.05, ** p < 0.01).
Insects 13 00296 g005
Table 1. Statistical analysis of the transcriptome sequence data.
Table 1. Statistical analysis of the transcriptome sequence data.
Clean BasesClean Reads Mapped Reads Mapped RatioQ20%Q30%GC (%)
5-AZA and bacteria injected8,970,065,91430,212,35523,677,53578.3797.8093.8143.95
5-AZA and bacteria injected7,182,585,60824,056,95319,384,24280.5897.5393.4143.00
5-AZA and bacteria injected7,613,303,87625,478,89220,665,63981.1197.5093.3442.21
Only bacteria injected7,333,560,22224,534,05619,994,00481.4997.4893.2040.95
Only bacteria injected7,653,093,91426,678,45620,316,29079.1297.9894.2142.52
Only bacteria injected7,018,186,15423,536,59618,881,66280.2297.7393.6942.58
Table 2. Distribution of splicing length.
Table 2. Distribution of splicing length.
Length Range Transcript Unigene
300–50035,064 (29.90%)26,435 (47.95%)
500–100028,596 (24.38%)14,437 (26.19%)
1000–200025,629 (21.85%)7620 (13.82%)
2000+27,981 (23.86%)6638 (12.04%)
Total Number117,27255,131
Total Length167,150,43153,875,332
N50 Length23571621
Mean Length1425.323977.2239
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Kausar, S.; Liu, R.; Gul, I.; Abbas, M.N.; Cui, H. Transcriptome Sequencing Highlights the Regulatory Role of DNA Methylation in Immune-Related Genes’ Expression of Chinese Oak Silkworm, Antheraea pernyi. Insects 2022, 13, 296. https://doi.org/10.3390/insects13030296

AMA Style

Kausar S, Liu R, Gul I, Abbas MN, Cui H. Transcriptome Sequencing Highlights the Regulatory Role of DNA Methylation in Immune-Related Genes’ Expression of Chinese Oak Silkworm, Antheraea pernyi. Insects. 2022; 13(3):296. https://doi.org/10.3390/insects13030296

Chicago/Turabian Style

Kausar, Saima, Ruochen Liu, Isma Gul, Muhammad Nadeem Abbas, and Hongjuan Cui. 2022. "Transcriptome Sequencing Highlights the Regulatory Role of DNA Methylation in Immune-Related Genes’ Expression of Chinese Oak Silkworm, Antheraea pernyi" Insects 13, no. 3: 296. https://doi.org/10.3390/insects13030296

APA Style

Kausar, S., Liu, R., Gul, I., Abbas, M. N., & Cui, H. (2022). Transcriptome Sequencing Highlights the Regulatory Role of DNA Methylation in Immune-Related Genes’ Expression of Chinese Oak Silkworm, Antheraea pernyi. Insects, 13(3), 296. https://doi.org/10.3390/insects13030296

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop