BIOCELL DOI:10.32604/biocell.2022.018166
Article

Identification of differential mRNA and lncRNA expression in AcMNPV-infected Sf9 cells

TIEJUN ZHAO, RIQIANG DENG*, MENGQIU CHEN and XUNZHANG WANG

School of Life Science, Sun Yat-sen University, Guangzhou, 510275, China
*Address correspondence to: Riqiang Deng, lssdrq@mail.sysu.edu.cn
Received: 04 July 2021; Accepted: 19 November 2021

Abstract: Sf9Sf9 are the ovarian cells of Spodoptera frugiperda that is the host of Autographa californica multiple nucleopolyhedrovirus (AcMNPV), and hence can serve as an effective test vehicle to understand the AcMNPV infection mechanism. In this study, through high-throughput sequencing technology using samples collected from Sf9 cells at different time points after AcMNPV infection, 3463 pieces of time-series differentially expressed RNA (1,200 mRNA and 2,263 lncRNA) are identified and justified by experimental verification of randomly selected samples from them, proving the validity of the bioinformatical analysis on this topic. Functional enrichment analysis and target prediction are performed on those differentially expressed RNA, from which the major functional enrichment distribution of those differentially expressed mRNA is derived. It has been found that the differential genes are mainly in the cellular anatomical entity and intracellular in terms of the cellular component, and in the binding and catalytic activity in terms of the molecular function. Also, the differential mRNA are mainly concentrated in global and overview maps, signal transduction, infectious diseases, and viral, etc. Moreover, those mRNA targeted by lncRNA are predicted. The correlation between those differentially expressed lncRNA and mRNA indicates that lncRNA is very likely playing an important role in the interaction between virus and host. Aided by an advanced co-expression analysis approach, the “hub” RNA is also identified. The study in this work pave the way for further analyzing and understanding how AcMNPV escapes from the host’s immunity, manipulates the host to realize the self-multiplication, and realizes the timely conversion between its two particle forms, laying the foundation for uncovering the host’s immune response process.

Keywords: AcMNPV; Sf9; lncRNA; mRNA; RNA-Seq

Abbreviations

Introduction

Baculoviruses are pathogenic microorganisms that infect insects and invertebrates, making them their host. The particle form of the virus is Baculoviridae (Vago et al., 1974). Its main hosts are Lepidoptera (corn borer, etc.), Hymenoptera (bees, etc.), and Diptera (mosquitoes, etc.) (Slack and Arif, 2007). They are divided into four genera, namely, Alphabaculovirus, Betabaculovirus, Gammabaculovirus, and Deltabaculovirus, according to the insect hosts from which they were isolated and their biological characteristics (Jehle et al., 2006). Autographa californica multiple nucleopolyhedrovirus (AcMNPV) is the archetypal species of the alphabaculovirus. During the infection cycle, it produces two enveloped virion phenotypes with different structures and functions, which play different roles in virus pathogenesis: budded virus (BV) and occlusion-derived virus (ODV) (Rohrmann, 2013). The nucleocapsid structures of these viruses are similar, but the origin and composition of their envelope and their roles in the virus life cycle are different.

Spodoptera frugiperda (Lepidoptera: Noctuidae) is a major migratory insect from tropical and subtropical areas of America and one of the major hosts for AcMNPV. As the ovarian cells of Spodoptera frugiperda, Sf9 is also the target to infect by AcMNPV. Many researches (Sparks, 1979; Johnson, 1987; Barr and Manning, 1997; Matindoost et al., 2015; Yin et al., 2008) on AcMNPV and its host has been started by infecting these Sf9 cells with AcMNPV.

The baculovirus establishes a strict parasitic and adaptive relationship with its insect host during the long evolutionary process, and the information of the interaction is recorded and solidified in the viral genome (Horton and Burand, 1993). For example, polyhedrins need to be alkali-interpreted by releasing virions to initiate infection in the pH-alkaline larval midgut, the alkaline environment of the insect host’s midgut and alkaline solubility of the viral polyhedrin protein crystal is the result of the long-term evolution of both.

Baculoviruses interact with their hosts at multiple levels, but the details remain unclear. The expression of baculovirus genes has a time sequence (Friesen and Miller, 1985, 1986, 1987; Nissen and Friesen, 1989). It has been reported that the host plays an important role in the expression of baculovirus genes. For example, the expression of early baculovirus genes depends on the host’s transcription system and factors (Ooi and Miller, 1988). After the baculovirus unshells and enters the nucleus, DNAs that have not yet been replicated are transcribed in the first place. The early transcription of the virus is inhibited by muscarine carnitine, so host RNA polymerase II mediates most early gene transcription, leading to the activation and expression of the late gene, thus making the expression of viral genes orderly (Tjia et al., 1979). Host RNA polymerase is the most important host component in the early gene expression of the virus and plays a dominant role in the life cycle of the virus. In addition, a host phosphoprotein with a relative molecular weight of 30 ku has been identified, and found to have high affinity with AcMNPV polyhedroid protein promoter (Burma et al., 1994). After extensive analysis, Williams and Faulkner (1997). found that the host nuclear membrane was related to the formation of viral envelope (Williams and Faulkner, 1997).

The baculovirus also affects its host. For example, one of the characteristics of baculovirus overexpression in early gene expression to late gene expression is significant inhibition of host transcription. Some studies (Knudson and Tinsley, 1974; Bilimoria et al., 1986) have shown that in Sf9 cells infected with AcMNPV, host protein synthesis begins to decline between 6 and 10 h after infection and seems to stop completely after 24 h. In addition, baculovirus can also inhibit the apoptosis of host cells and affect the cycle process of host cells. Bertin et al. (1996) demonstrated that the P35 protein of baculovirus inhibits apoptosis induced by baculovirus infection. One of the functions of AcMNPV’s very early gene, IE2, has been shown to encode proteins that block cell cycles in various cell lines.

The interaction between baculovirus and host is reflected at multiple levels, from simple physical adsorption (Granados, 1978; Summers, 1971), invasion to complex manipulation of host metabolic system (Funk and Consigli, 1993), and induction of host cell apoptosis (Zhang et al., 2002, 2019), etc. Many interesting and important host-virus interaction mechanisms have been identified from fundamental studies of the biochemistry and molecular biology of baculoviruses (Manji, 2000; Miller and Lu, 1997; Cory et al., 2005; Joshi et al., 2010; Miller, 1997). These studies, to some extent, provide new strategies for biological pest control and eukaryotic expression vector systems.

A large part of functional genes of baculoviruses have been identified thus far, but the research on the relationship between the two main focuses on viruses, and there is less research on the function of host factors, so the work in this field is worthy of further development (Rohrmann, 2011; Friesen and Miller, 1986). In addition, many previous studies (Mccarthy and Theilmann, 2008; Cochran, 1983; Fang et al., 2009; Wu et al., 2018) have carried out relevant findings from the perspective of host genes or viral gene functions, while from the perspective of bioinformatics, there is no relevant report on the changes in the expression level of host genes after baculovirus infection. In this study, after the baculovirus infecting Sf9 cells, infected cells were collected at different time points for sequencing, and the sequencing results were analyzed to identify those messenger RNA (mRNA) and long non-coding RNA (lncRNA) that show expression features of sequential difference. Moreover, functional enrichment, target and co-expression analysis were performed to shed light on the correlation of mRNA and lncRNA. The molecular interaction between viral genes (and their products) and host cytokines are so exquisitely coordinated that even subtle changes can inhibit viral replication. This study provides some candidate genes that are closely related to the molecular regulation mechanism of Sf9 cells infected with AcMNPV. Our results will provide a useful resource for further analysis of baculovirus-host interactions, and understand the immune escape mechanism of the virus.

Material and Methodology

Cells and insects

Sf9 (Spodoptera frugiperda IPLB-Sf21-AE clonal isolate 9) insect cells were cultured in Grace’s medium supplemented with 10% fetal bovine serum at 27°C. Groups of Sf9 cells were infected with AcMNPV at a multiplicity of infection (MOI) of 10, and subsequently, RNA samples were collected at 6, 12, 24, 36, 48, 60 and 72 h post-infection (hpi).

RNA isolation

Total RNA was isolated from approximately 107 infected cells harvested using TRIzol reagent (Invitrogen, USA) according to the manufacturer’s protocol, and purified by using RNeasy mini kit (Qiagen, GmBH, Germany). Total RNA was dissolved in 50 μL RNase-free H2O. The RNA solution was digested for 30 min at 37°C with RNase-free DNase I (TaKaRa, Japan) to remove any contaminating genomic DNA. RNA integrity was analyzed by gel electrophoresis using a 1% agarose gel alongside an RNA marker. The samples were sent to Beijing Genomics Institute Shenzhen Co. (BGI-Shenzhen) (China) for sequencing using Solexa technology (Cuddapah et al., 2009; Nagalakshmi et al., 2008; Rougemont et al., 2008). In the subsequent sequencing test, RNA quality was checked at first. Only when the RNA integrity value (RIN) is larger than 8, indicating RNA quality and integrity meet the sequencing requirement, subsequent test could be continued. In this study, Illumina Novaseq 6000 was used for sequencing test, with a depth of 6G (20 M reads).

Raw data processing and transcriptome

Fastp (version 0.19.4) (Chen et al., 2018) was used to filter low-quality reads, cut adapters and quality control of raw FASTQ files to obtain clean reads. The clean reads of each sample were aligned against SLIVA (Quast et al., 2013) ribosomal RNA and AcMNPV genome (accession number NC_001623) with Bowtie2 (Langmead and Salzberg, 2012) to filter out ribosomal RNA and AcMNPV RNA.

Then, the unmapped reads of each sample were aligned against the Spodoptera frugiperda genome (NCBI, GCA_002213285.1) with HISAT2 (version 2.1.0) (Kim et al., 2015) and then subsequently assembled by StringTie (version 1.3.4d) (Pertea et al., 2015) separately. All the assemblies were merged into one transcriptome by TACO (Niknafs et al., 2017). Then, Salmon (version 0.11.2) (Patro et al., 2017) was applied to quantify the transcript expression of the transcriptome. To reduce noise, we include only transcripts with TPM > 0.1 and read count > 10.

BLASTP was used to compare the obtained sequences with the NR database, and the unmatched sequences were selected for subsequent analysis. PLEK (Li et al., 2014) was used to calculate the coding potential of the transcripts. Transcripts that had no coding potential, and harboured at least two exons, and had a length larger than 200 nt, were classified as lncRNAs. Transcripts that showed coding potential and/or had a similar sequence with proteins were classified as mRNAs. R package EBSeqHMM (Leng et al., 2015) was used to identify transcripts changing over time; FDR < 0.01 was considered as signature. The detailed flow is illustrated in Fig. 1.

Figure 1: LncRNA and mRNA sequences were separated and assembled by processing RNA-seq data through a series of analytical processes. First, the raw data were filtered to remove rRNA and viral RNA. The obtained clean data were then compared with the Sf9 genome. Sequences with TPM expression levels greater than 0.1 and read count greater than 10 were compared with the NR database. The unmatched sequences were further analyzed using PLEK to differentiate lncRNA and mRNA.

qRT-PCR analysis

Real-time quantitative Reverse Transcription-PCR (qRT-PCR) was performed to define the expression profile of mRNA and lncRNA. We compared 57,104 RNAs with eggNOG and other databases to find the sequence of selected genes and designed relevant primers based on the sequence information of selected genes. The mRNA and 5S rRNA were reversely transcribed using a random primer (Takara) and included in the respective sample. The qRT-PCR results were analyzed by utilizing the 2–ΔΔCT method. The experiment was repeated three times and averaged to show changes in the expression of host genes at different time points after infection.

GO and KEGG functional enrichment analysis methods

We performed a Gene Ontology (GO) analysis for biological processes, cellular components and molecular function and a KEGG pathway analysis (Kyoto Encyclopedia of Genes and Genomes http://www.genome.ad.jp/kegg) via enrich R package.

RNAplex

RNAPlex software is used in this work to predict the interaction between lncRNA and mRNA. It mainly calculates the minimum free energy according to the thermodynamic structure to predict the optimal base-pairing relationship (Tafer and Hofacker, 2008).

WGCNA

The weighted gene co-expression network analysis (WGCNA) technique is designed to search gene modules that are co-expressed and explore the correlation between gene network and targeted expression types and the core genes in the network (Langfelder and Horvath, 2008). This analysis applies to complex data models. For example, the response at different time points after pathogen infection can be well analyzed. In this study, with the help of this analysis, the lncRNA with time-series difference and the “hub” RNA will be identified.

Results

Overview of RNA sequencing (RNA-seq)

After standardizing the raw data, a total of 24,230,852.75 clean reads that match the Spodoptera frugiperda genome were obtained through sequencing and filtering out ribosomal RNA and AcMNPV RNA of samples of the 8-time points.

Identification of lncRNA and mRNA

The obtained clean reads were assembled into 57,104 non-redundant RNA transcripts, including 19,982 lncRNA and 37,122 mRNA, as shown in Figs. 2A–2C. In general, GC content and length can be well differentiated.

Figure 2: Isograph of GC content, length, and gene expression of lncRNA and mRNA. RNA meeting the requirement that its transcripts present in at least 2 samples, counts > 10, TPM > 0.1 for each sample, was used. lncRNA and mRNA can be distinguished from each other in all three aspects. A: GC percentage distribution comparison between lncRNA and protein coding RNA transcripts. B: Length distribution comparison between lncRNA and protein coding RNA transcripts. C: Cumulative gene expression (TPM) distribution of lncRNA and protein coding RNA at 0 h. The horizontal axis is log10 (TPM) at 0 h, and the vertical axis is the accumulation of its probability density.

Sequence expression differences of infected Sf9 cells

Aided by further EBSeqHMM analysis, it was found that 1,200 mRNA and 2,263 lncRNA showed differences in temporal expression after AcMNPV infection. Furthermore, we randomly selected 4 RNA, including 1 mRNA and 3 lncRNA, for further experimental validation.

As shown in Fig. 3, the overall trend of the selected lncRNA and mRNA was the same whether the results were analyzed by EBSeqHMM or by qRT-PCR. After virus infection, XP_022825542.1, XP_022823371, XP_022827228.1, and XP_022814930.1 exhibit an overall downward trend in both the bioinformatics and qRT-PCR predictions. In contrast, XP_022827228.1 exhibit an overall upward trend in both the bioinformatics and qRT-PCR predictions.

Figure 3: Comparison between bioinformatics and qPCR results of lncRNA XP_022825542.1, XP_022823371, and XP_022827228.1, and mRNA XP_022814930.1, showing an overall downward trend for XP_022825542.1, XP_022823371, and XP_022814930.1, and an overall upward trend for XP_022827228.1, in both the bioinformatics and qRT-PCR predictions. The bioinformatics analysis was carried out by EBSeqHMM, and qPCR is the experimental methodology. After the cells were infected with the virus, cells at the time points involved in the bioinformatics analysis were collected to extract RNA for reverse transcription, and then qPCR was performed.

The result indicates that among the obtained 19,982 lncRNA, 2,263 showed differences in temporal expression, accounting for 11.33%, while only 1,200 of the obtained 37,122 mRNA showed differences in temporal expression, accounting for 3.23%, implying that many non-coding RNA played a role in the interaction between viruses and hosts.

GO and KEGG functional enrichment outcome

We analyzed GO and KEGG for the differential mRNAs, as shown in Fig. 4. In the GO analysis, the biological process difference of these differential mRNAs is mainly reflected on the biological regulation, and cellular and metabolic processes. There are 188, 177, and 64 genes exhibiting sequential differential expressions in the metabolic process, cellular process, and biological regulation, respectively. We found that the differential gene is mainly in the cellular anatomical entity and intracellular in terms of the cellular component, and in the binding and catalytic activity in terms of the molecular function. There are 159, 94, 209, and 186 genes exhibiting sequential differential expressions in the cellular anatomical entity, intracellular, binding, and catalytic activity, respectively.

Figure 4: A: GO enrichment analysis of differentially expressed mRNA. B: KEGG enrichment analysis of differentially expressed mRNA.

In the KEGG enrichment analysis, we found that the differential mRNAs were mainly concentrated in global and overview maps, signal transduction, and infectious diseases (viral), etc. There are 143, 141, and 89 genes exhibiting sequential differential expressions in the global and overview maps, signal transduction, and infectious diseases (viral), respectively.

RNAplex analysis outcome

In this study, differentially expressed lncRNAs and mRNAs were obtained through analysis, and these differentially expressed mRNAs may be directly regulated by lncRNAs. Thus, this study analyzed the targets of differential lncRNAs in these differential mRNAs with the help of RNAPlex. As shown in Fig. 5A, the target gene of the lncRNA TU187 (XP_022825542.1) is TU42049. However, due to incomplete annotation of SF9 gene, the function of this gene is not clear. Figs. 5B–5D shows that there are three target genes of the lncRNA TU23998 (XP_022823371.1), one is TU20268, which is a UDP-glycosyltransferase, the second is TU42142, which is an embryonic development factor, and the third is TU50365, which is a predicted homologue with ATPase. Therefore, its target is related to the energy metabolism and growth of the host. However, the target of XP_022827228.1 was not found in the mRNAs with time-series variation characteristics obtained in our analysis this time. This may be because we involved more time points and narrow conditions, and it was only found in the genes with differentially changed genes to find its target.

Figure 5: A: TU187 lncRNA and its target gene TU42049. B–D: TU23998 lncRNA and its 3 target genes TU20268, TU42142, and TU50365.

WGCNA analysis outcome

With the help of WGCNA, we analyzed the co-expression of these different lncRNAs and mRNAs. First, the results of WGCNA were analyzed, and the genes of each cluster were enriched by GO and Pathway analysis. Weight > 0.6 is used to filter the interaction relationship pairs, and then the interaction network graph is built up with Cytoscape software, as shown in Fig. 6. We can further see that some important RNAs are differentially expressed, such as TU32028 and TU33799 in mRNA. The former is a ribonucleoside diphosphate reductase subunit M2, which is an important gene associated with DNA synthesis, and the latter is a phosphodiesterase, which plays an important role in signal transduction and other processes. Among lncRNAs, TU29511 and TU42509 were co-expressed by GTP-binding proteins and ethanolamine kinase in the former, and casein kinase ATP-binding proteins in the latter.

Figure 6: WGCNA Analysis shows the Co expression network. Note that original vector file has been uploaded as a supplementary material.

Discussion

Spodoptera frugiperda is extremely harmful to crops. In 2016, for the first time, Spodoptera frugiperda was found in Africa and rapidly swept across 44 countries in the south of the Sahara within two years (Goergen et al., 2016; Chapman et al., 2017). According to the Centre for Agriculture and Bioscience International (CABI), a survey of 12 African maize-producing countries shows that Spodoptera frugiperda is responsible for losses of between 8.3 million and 20.6 million tons, or $2.5 billion to $6.2 billion USD, of maize production each year. In addition, more than $13 billion USD of crops are still at risk, causing major damage to agricultural production in Africa (Stokstad, 2017; Early et al., 2018). In 2017, the European Food Safety Authority (EFSA) listed Spodoptera frugiperda as a quarantine pest (Jeger et al., 2017). The Food and Agriculture Organization (FAO) of the United Nations has also issued a global warning, putting it on the world’s top 10 blacklisted plant pests (Food and Agriculture Organization of the United Nations, 2018). The work on prevention and control has not stopped, and biological control, as an important method, has also attracted researchers’ attention. Great progress in preventing and controlling Spodoptera frugiperda has been achieved by utilizing viruses, bacteria, nematode and other pathogenic microorganisms, such as Spodoptera frugiperda multiple nucleopolyhedrovirus (SFMNPV) and Beauveria bassiana et. (Roger et al., 2017; Casmuz et al., 2010). When the pathogenic microorganisms, such as B. bassiana, Metarhizium rileyi and Caenorhabditis elegans, were mixed with insecticides, the prevention and control effect of Spodoptera frugiperda could be significantly improved.

Researchers have successively found the encoded mRNA and function of the AcMNPV gene. AcMNPV-miR-1 is considered to target and down-regulate the expression of viral gene odv-e25, accelerate polyhedra formation, and promote viral infection efficiency in Trichoplusia ni larvae (Zhu et al., 2013, 2016). AcMNPV-miR-3 is located on the opposite strand of the viral gene ac101 coding sequence in the AcMNPV genome, and it can be detected at 6 h post-infection and accumulated to a peak around 12 h post-infection in AcMNPV-infected Sf9 cells. Five viral genes (ac101, ac23, ac25, ac86, and ac98) were verified to be regulated by AcMNPV-miR-3 (Jiao et al., 2019). As a parasite of Spodoptera frugiperda, AcMNPV is essential in the biological control of Spodoptera frugiperda in a safe and efficient way. Thus, in-depth study from gene level is an important step to understand their interaction.

In this study, by sequencing samples collected at different time points after infecting the host cells with AcMNPV, those time-series differentially expressed mRNA and lncRNA are identified. Four types of RNA are randomly selected and verified against experiment, proving that the experiment results is consistent with bioinformatics analysis, and hence that the bioinformatics analysis/methodology employed in this work is valid and reliable. Through GO enrichment analysis, we found that in the aspect of biological process, the differentially expressed genes are concentrated in the biological process, cellular process, metabolic process and regulation of biological process. This finding indicates that after the virus infecting hosts, the cell growth will be impacted. Knudson and Tinsley (1974) observed in their study that after the host is infected with baculovirus, DNA synthesis of the host is quickly inhibited. In the aspect of cellular component, they are mainly in the cellular anatomical entity, indicating that after infection with baculovirus, the host body composition is greatly affected. Volkman et al. (1986) and Monsma et al. (1996) pointed out that the transmission of BV, a particle form of baculovirus, from cell to cell in the host depends on endocytosis, during which the cell membrane, one of the physical components of the host cell, changes. Binding and catalytic activity are found to be the two major aspects of molecular function, indicating that after infection, the activities of many enzymes of the host are affected, thus affecting the normal life process of the host. For instance, it has been reported that the specific structure of AcMNPV’s homologous region 1 (hr1) enables it to bind to a 38 ku (kilo units) host nuclear protein, which may play an important role in enhancer function of hr1 (Kim et al., 2015).

In the KEGG enrichment analysis, it can be found that the time-series differentially expressed mRNA is enriched in global and overview maps, and signal transduction, etc., indicating that by regulating the signal transduction, the virus regulates the host so as to evade immunity and realize its own proliferation. The transformation of the two morphologic particles of the virus may also be closely related to signal transduction. As pointed out by many researchers, signal transduction is closely related to cell proliferation (Wiernas et al., 2010; Seuwen and Pouysségur, 1992). It has also been found in the study to some virus that signal transduction is essential in the growth and activity of the virus (Kieser, 2010; Vainionpää et al., 1991; Gupta and Vayuvegula, 1987; Hijikata et al., 1999). There has been limited research on the signal transduction of virus that infect insets. The results in this study indicates that signal transduction plays an important role in the self- proliferation after the host being infected by AcMNPV. The high enrichment result in this study indicates that the virus regulates the signal transduction of the host after it infects the host, thereby changing the cell function and directing it to the direction that is conducive to its own proliferation.

Based on GO and KEGG analysis, we obtained the differential expression of the host after virus infection at the macroscopic level. Instead of the conventional comparison of the difference before and after infection, we monitored and recorded the dynamically varying difference post infection. Those filtered differential genes could play important roles in the procedure of virus infecting hosts and realizing self-proliferation. For example, the 64 differential genes, enriched in the biological regulation, have high chance to impact the host life activity and hence the virus proliferation.

It has been reported that after AcMNPV infects the host, host DNA synthesis is quickly inhibited (Knudson and Tinsley, 1974), host transcription will be significantly suppressed (Ooi and Miller, 1988), and host protein synthesis begins to decline at 6 h–10 h after infection, and seems to stop completely after 24 h (Carstens et al., 1979; Wood, 1980). At later stages the cell ruptures, releasing polyhedra (Blissard, 1996). However, further studies and understanding are needed on what the status of the regulation network is, which mRNA and lncRNA are involved. For example, researchers already know that BV enters cells through endocytosis (Volkman and Goldsmith, 1985; Habib et al., 1996), but the details and molecular mechanisms of this process need to be further studied. In this study, added by RNAplex analysis, we obtained correlation between differentially expressed lncRNA and mRNA, laying the foundation for future exploration of the regulation network.

Aided by WGCNA, the co-expression of those differential lncRNA and mRNA was analyzed further, and hub RNA are identified, such as TU32028 in mRNA, which is essential in the DNA synthesis, and is part of the p53 signaling pathway, and TU29511 and TU42509 in lncRNA, which are related to GTP-binding protein and ethanolamine kinase etc., and casein kinase ATP binding protein, etc., respectively. Apparently, those differentially expressed RNA occupy an important position in the signal transduction. For the host, it also needs to rely on a variety of signal transduction pathways to meet its own needs.

How the virus uses these transduction pathways, so that the host evolves toward the direction that benefits the virus to achieve its proliferation is worth discussing. Bioinformatics analysis in this study provides a prerequisite for such a discussion.

The molecular interaction between virus and host is very delicate and coordinated. The smooth expression of virus genes is not only related to the virus itself, but also closely related to host factors. The understanding on how this regulation is conducted, and which gene plays the key role become the pre-requisite of uncovering the mechanism of this phenomena.

Although researchers have conducted many studies of AcMNPV and its host, at present the understanding of the regulating mechanism of a series of virus life process in the host cells is relatively limited, including the temporal transcriptional regulation mechanism of host genes during infection. The sequencing and annotation of host’s gene group were started relatively late. In 2014, Kakumani et al. (2014) sequenced the genome of Sf21, which is the separated ovarian cell line of Spodoptera frugiperda, by using the second-generation sequencing technique, which was the first genome data of Sf2 (Kakumani et al., 2014). In 2017, scientists resequenced maize and rice strains of the same ovarian cell line (Sf21), yielding a genome of about 396 Mb (Gouin et al., 2017). In the same year, Nandakumar et al. (2017) also sequenced the genome of another cell line (Sf9) of Spodoptera frugiperda, which obtained a total of 451 Mb of genome data. The assembly result and sequencing quality were better than the former two. In 2019, Liu et al. (2019) sequenced the entire genome of Spodoptera frugiperda that had invaded China and assembled the sequencing results to the chromosome level. The samples were collected from Yunnan and Guangdong provinces in China. This study obtained Spodoptera frugiperda genome size of about 536 Mb, with N50 reaching more than 14 Mb, and more than 80% of the genome data were assembled on the 31 chromosomes. In other words, the gene sequencing on hosts has not been paid sufficient attention until recent years, and annotation work lacks detailed report, which also prevents further in-depth analysis in this study. As the annotation work continues to improve, the data and analysis in this study serve as the first step toward a thorough understanding of the interaction between virus and hosts.

Conclusion

By high-throughput sequencing of RNA from Sf9 cells infected by AcMNPV, we found 1,200 mRNA and 2,263 lncRNA showing different temporal expression. Results from the experimental group are basically identical with those from the bioinformatic analysis. The overall trend is consistent, though there are differences at individual points in time, which may be associated with the variation in the repetition count of the experiment and sequencing analysis, and can be improved by more sequencing in further studies. Through GO and KEGG analysis, we obtained the enrichment information of differential genes, which enabled us to have a certain understanding of the classification of differential genes and the complex interrelationship between genes and metabolites from a macro perspective. This provides the opportunity to further study the genes that play a key role in host-virus interaction and contribute to a comprehensive understanding of the interaction between virus and host. We also analyzed the target genes of the obtained differential lncRNAs in the differential mRNAs by RNAplex. The last but not the least, we analyzed the co-expression of the differential lncRNAs and mRNAs by WGCNA, and identified those “hub” RNAs. This study paves the way for further in-depth analysis and understanding of the interaction between virus and host post AcMNPV infection, and the regulation mechanism of virus involved life processes.

Availability of Data and Materials: The datasets generated during and/or analysed during the current study are available in the NCBI database (PRJNA736178).

Authors’ Contribution: Tiejun Zhao and Mengqiu Chen performed the experiment and analysis. Tiejun Zhao, Riqiang Deng and Xunzhang Wang wrote the manuscript.

Funding Statement: The authors received no specific funding for this study.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.