Correlation in Expression between LTR Retrotransposons and Potential Host Cis-Targets during Infection of Antherea pernyi with ApNPV Baculovirus

The published genome sequence of Antheraea yamamai (Saturnnidae) was used to construct a library of long terminal repeat (LTR)-retrotransposons that is representative of the wild silkmoth (Antherea) genus, and that includes 22,666 solo LTRs and 541 full-length LTRs. The LTR retrotransposons of Antheraea yamamai (AyLTRs) could be classified into the three canonical groups of Gypsy, Copia and Belpao. Eleven AyLTRs contained the env gene element, but the relationship with the env element of baculovirus, particularly A. yamamai and pernyi nucleopolyhedrovirus (AyNPV and ApNPV), was distant. A total of 251 “independent” full-length AyLTRs were identified that were located within 100 kb distance (downstream or upstream) of 406 neighboring genes in A. yamamai. Regulation of these genes might occur in cis by the AyLTRs, and the neighboring genes were found to be enriched in GO terms such as “response to stimulus”, and KEGG terms such as “mTOR signaling pathway” among others. Furthermore, the library of LTR-retrotransposons and the A. yamamai genome were used to identify and analyze the expression of LTR-retrotransposons and genes in ApNPV-infected and non-infected A. pernyi larval midguts, using raw data of a published transcriptome study. Our analysis demonstrates that 93 full-length LTR-retrotransposons are transcribed in the midgut of A. pernyi of which 12 significantly change their expression after ApNPV infection (differentially expressed LTR-retrotransposons or DELs). In addition, the expression of differentially expressed genes (DEGs) and neighboring DELs on the chromosome following ApNPV infection suggests the possibility of regulation of expression of DEGs by DELs through a cis mechanism, which will require experimental verification. When examined in more detail, it was found that genes involved in Notch signaling and stress granule (SG) formation were significantly up-regulated in ApNPV-infected A. pernyi larval midgut. Moreover, several DEGs in the Notch and SG pathways were found to be located in the neighborhood of particular DELs, indicating the possibility of DEG-DEL cross-regulation in cis for these two pathways.


Introduction
Sericulture is one of the great inventions of the ancient Chinese, and it has created enormous economic benefits for society by the art of silk production. Chinese oak silkworm, Antheraea pernyi

Identification and Characterization of LTR-Retrotransposons
To identify LTR-retrotransposons, the A. yamamai genome was used as a reference genome and downloaded from the NCBI database (project accession PRJNA383008 and PRJNA383025) [15]. LTR-retrotransposons were identified in the A. yamamai genome using the software LTRharvest (GenomeTools1.5.7) and LTRdigest [17,18]. Pairs of putative LTRs that were separated by 1-15 kb and flanked by target site duplications were screened by LTRharvest in A. yamamai genome. The threshold of LTR nucleotide similarity used in LTRharvest was set at higher than 80%; other parameters were set to their defaults. Internal features of LTR-retrotransposons, including protein domains, primer-binding sites, and polypurine tracts, were predicted using LTRdigest with default setting.
In the present study, LTR-retrotransposons which contain at least one relevant protein domain between the pairs of putative LTRs are called full-length LTR-retrotransposons. Correspondingly, LTR-retrotransposons lacking all protein domains are called solo LTR-retrotransposons.

Classification and Phylogenetic Analysis of LTR-Retrotransposons
The superfamily classification of LTR-retrotransposons of A. yamamai (AyLTR) is based on the homology between their reverse transcriptase (RT) domain sequences and the sequences of Peptidase_A17, RVT_1 and RVT_2 from the Pfam database [19]. Sequences of the RT domain from LTR-retrotransposons were used for multiple alignment using MUSCLE (v3.8.31) [20]. A phylogenetic tree built based on the neighbor-joining algorithm was generated from the RT domain alignment using MEGA5 with 1000 bootstrap replicates.
In addition, according to the positional relationship between full-length LTR-retrotransposons and host genes, LTR-retrotransposons are divided into three categories in our study. Full-length LTR-retrotransposons that do not overlap with gene exon sequences are defined as independent LTRs ("Stream"); full-length LTR-retrotransposons that partially overlap with a gene exon sequence are defined as partially overlapping LTRs ("Part"); and full-length LTR-retrotransposons that encompass gene exon sequences are defined as overlapping LTRs ("In").

Analysis of the Env Genes from AyLTRs
AyLTRs possessing env-like elements were selected from the library of full-length LTR-retrotransposons. The fusogenic region of AyLTR envelope glycoproteins was analyzed using MegAlign (DNASTAR Lasergene.v7.1) and Weblogo (http://weblogo.berkeley.edu/logo.cgi). Phylogenetic analysis of LTR-retrotransposon env sequences from different insect species was performed using MegAlign and MEGA5 (5.0.1.102). In addition, sequence alignments were performed between AyLTR env genes and genes that code for envelope proteins (F and GP64) in Group I and Group II nuclear polyhedrosis viruses (NPVs) using MegAlign and MEGA5 with 1000 bootstrap replicates.

Identification of Genes Located in the Neighborhood of AyLTRs in the A. yamamai Genome
The UCSC Genome Bioinformatics tool was used to screen for genes located within 100 kb of upstream and downstream of full-length AyLTR elements. Target genes are defined as genes that have annotated exons (UTR and CDS) within the defined sequence space of 100 kb.
Blast2GO and WEGO were used to perform gene ontology (GO) classification. The GO terms included molecular function, cellular component and biological process. For the pathway enrichment analysis, the genes were mapped to Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The hypergeometric test was used to identify overrepresented KEGG pathway and GO terms with a significance level of p < 0.05.

Analysis of RNA-Seq Data
To analyze the LTR-retrotransposon transcriptome in A. pernyi infected with ApNPV, we downloaded the raw reads from NCBI with the accession numbers: SRR2919240, SRR2919241, SRR2919242 and SRR2919243 [3]. In the study that provided the transcriptome raw data for our analysis [3], A. pernyi larvae infected with 4.05 × 10 6 polyhedra/mL ApNPV for three days were the experimental group while the control group was treated with the same volume of 0.9% NaCl physiological saline. Each group included two independent biological replicates. The four cDNA libraries were designated as Ap_CK1, Ap_CK2, Ap_NPV1 and Ap_NPV2 [3].
These raw reads were mapped to the A. yamamai genome using TopHat2 (version 2.0.3.12) [22]. The mapped reads of each sample were assembled by software Cufflinks and Cuffmerge. Gene and LTR-retrotransposon abundances were quantified by software RSEM [23] and their expression level was normalized by FPKM (Fragments Per Kilobase of transcript per Million mapped reads).
To identify differentially expressed genes (DEGs) and differentially expressed LTR-retrotransposons (DELs), the edgeR package (http://www.rproject.org/) was used. Genes and LTR-retrotransposons with fold change of |log 2 FC | > 1 (FC: fold change) and a false discovery rate (FDR) <0.05 were considered as DEGs and DELs. The interaction networks of DELs and neighboring DEGs were imported to Cytoscape software for visualization.

De novo Detection of AyLTRs in the Antheraea yamamai Genome
AyLTRs were detected by LTRharvest and annotated with LTRdigest. A total of 23,207 AyLTRs were identified in the A. yamamai genome. Among these are 541 full-length AyLTRs and 22,666 solo AyLTRs. Furthermore, we identified 11 AyLTRs containing the env ORF. The sizes of the 541 full-length AyLTRs varied from 1127 to 15,130 bp and the lengths of the LTR ranged from 100 to 1000 bp. Details of the AyLTRs are shown in Supplementary Materials Table S1.

Phylogenetic Analysis and Classification of Full-Length AyLTRs
To classify the full-length AyLTRs, their RT domain sequences and those of other insect LTRs were used to build a multiple alignment and compute phylogenetic trees using the neighbor-joining method of MEGA 5. Among the 541 full-length AyLTRs, 189 AyLTRs had an RT domain that was sufficiently conserved for confident alignment during phylogenetic analysis. AyLTRs were classified into the three canonical groups Gypsy, Copia and Belpao ( Figure 1). Most of the AyLTRs included in the phylogenetic analysis belonged to the Gypsy group ( Figure 1).

Analysis of Env Genes from Full-Length AyLTRs
Eleven full-length AyLTRs containing the env gene element were identified in the A. yamamai genome (Figure 2A). Only four AyLTRs, namely AY-408, AY-35, AY-318 and AY-787, contained the complete structure (LTR-gag-pro-pol-env-LTR) of an insect retrovirus ( Figure 2A). We further found that the env amino acid sequences of AY-37, AY-58, AY-87 and AY-476 shared a region of similarity with the fusion proteins of known insect ERVs and Group II NPVs, such as the furin cleavage signal and the downstream fusion peptide ( Figure 2B). Logo representation of the furin-like consensus motif is RxxR and the peptide fusion consensus sequence is GxxxxxGxxxKxxxGxxDxxD ( Figure 2B). However, the furin cleavage site located in front of the fusion peptide in AY-58 and AY-87 is incomplete ( Figure 2B). On the other hand, the results suggest that the env genes of AY-37 and AY-476 have the potential to encode fusion proteins.

Analysis of Env Genes from Full-Length AyLTRs
Eleven full-length AyLTRs containing the env gene element were identified in the A. yamamai genome (Figure 2A). Only four AyLTRs, namely AY-408, AY-35, AY-318 and AY-787, contained the complete structure (LTR-gag-pro-pol-env-LTR) of an insect retrovirus ( Figure 2A). We further found that the env amino acid sequences of AY-37, AY-58, AY-87 and AY-476 shared a region of similarity with the fusion proteins of known insect ERVs and Group II NPVs, such as the furin cleavage signal and the downstream fusion peptide ( Figure 2B). Logo representation of the furin-like consensus motif is RxxR and the peptide fusion consensus sequence is GxxxxxGxxxKxxxGxxDxxD ( Figure 2B). However, the furin cleavage site located in front of the fusion peptide in AY-58 and AY-87 is incomplete ( Figure 2B). On the other hand, the results suggest that the env genes of AY-37 and AY-476 have the potential to encode fusion proteins. Phylogenetic analysis shows that the nucleotide sequences of env of AY-318, AY-35 and AY-408 are closely related with env of Trichoplusia ni TED virus. The AY-276 env was closely related to BmERV 94 env. In addition, close phylogenetic relationships of env genes were also observed between AY-37, AY-58, AY-476 and Drosophila melanogaster tirant and ZAM ERVs ( Figure 2C).
To explore the evolutionary relationship between AyLTR env genes and genes encoding the envelope protein (Env) in Group I and II NPVs, corresponding phylogenetic trees were generated. The results showed that the env sequences of AyLTRs were distantly related to Fa, Fb and gp64 genes from AyNPV, ApNPV, as well as other Group I and Group II NPVs ( Figure 2D, E). On the other hand, AY-483 env was clustered on the same branch of the phylogenetic tree with HaSNPVgOrf133 Fa ( Figure 2D). Phylogenetic analysis shows that the nucleotide sequences of env of AY-318, AY-35 and AY-408 are closely related with env of Trichoplusia ni TED virus. The AY-276 env was closely related to BmERV 94 env. In addition, close phylogenetic relationships of env genes were also observed between AY-37, AY-58, AY-476 and Drosophila melanogaster tirant and ZAM ERVs ( Figure 2C).
To explore the evolutionary relationship between AyLTR env genes and genes encoding the envelope protein (Env) in Group I and II NPVs, corresponding phylogenetic trees were generated. The results showed that the env sequences of AyLTRs were distantly related to Fa, Fb and gp64 genes from AyNPV, ApNPV, as well as other Group I and Group II NPVs ( Figure 2D, E). On the other hand, AY-483 env was clustered on the same branch of the phylogenetic tree with HaSNPVgOrf133 Fa ( Figure 2D).

Identification and Analysis of (Potentially Cis-Target) Genes that Occur in the Neighborhood of AyLTRs
To identify host genes that may be regulated by LTR-retrotransposons, we analyzed the positional relationship between AyLTRs and their neighboring genes. A total of 141 "In" AyLTRs were identified that encompassed exons of host genes and 141 "Part" AyLTRs were detected that overlapped with neighboring exons ( Figure 3A). Among these AyLTRs, six AyLTRs belong to both "In" and "Part" categories. In addition, 251 full-length "Stream" AyLTRs were identified as independent LTRs ( Figure 3A). 14 AyLTRs are located on scaffolds that do not contain genes.
"Stream" AyLTRs were further examined whether their expression was coordinated with neighboring cellular genes which could be indicative for regulation of expression of cellular genes in

Identification and Analysis of (Potentially Cis-Target) Genes that Occur in the Neighborhood of AyLTRs
To identify host genes that may be regulated by LTR-retrotransposons, we analyzed the positional relationship between AyLTRs and their neighboring genes. A total of 141 "In" AyLTRs were identified that encompassed exons of host genes and 141 "Part" AyLTRs were detected that overlapped with neighboring exons ( Figure 3A). Among these AyLTRs, six AyLTRs belong to both "In" and "Part" categories. In addition, 251 full-length "Stream" AyLTRs were identified as independent LTRs ( Figure 3A). 14 AyLTRs are located on scaffolds that do not contain genes.
"Stream" AyLTRs were further examined whether their expression was coordinated with neighboring cellular genes which could be indicative for regulation of expression of cellular genes in cis by the LTRs [24,25]. Applying a method for finding cis-target genes of lncRNAs [26], 406 genes located within 100 kb upstream and downstream of "Stream" AyLTRs were identified as potential cis-target genes of these AyLTRs. GO biological process analysis showed that most of the candidate cis-target genes were enriched in "cellular process", "metabolic process", "response to stimulus" and "single-organism process" ( Figure 3B). KEGG analysis illustrated mainly enrichment in "oxidative phosphorylation", "RNA transport", "Protein processing in endoplasmic reticulum", "mTOR signaling pathway" and several metabolic pathways ( Figure 3C); see Supplementary Materials Table S2 for more details). cis by the LTRs [24,25]. Applying a method for finding cis-target genes of lncRNAs [26], 406 genes located within 100 kb upstream and downstream of "Stream" AyLTRs were identified as potential cis-target genes of these AyLTRs. GO biological process analysis showed that most of the candidate cis-target genes were enriched in "cellular process", "metabolic process", "response to stimulus" and "single-organism process" ( Figure 3B). KEGG analysis illustrated mainly enrichment in "oxidative phosphorylation", "RNA transport", "Protein processing in endoplasmic reticulum", "mTOR signaling pathway" and several metabolic pathways ( Figure 3C); see Supplementary Data Sheet S2 for more details). According to the positional relationship between AyLTRs and their neighboring genes, AyLTRs were classified into three types, "In", "Part" and "Stream". "In" represent AyLTRs that encompass gene exon sequences. "Part" represents AyLTRs that partially overlap with exon sequences. "Stream" represents independent AyLTRs that do not overlap with host gene exons. (B) Histogram description of Gene Ontology enrichment of neighboring genes of "Stream" AyLTRs.
Genes were assigned to three categories: biological process (BP), cellular component (CC), and molecular function (MF). (C) The neighboring genes of "Stream" AyLTRs that were enriched in KEGG pathways. Genes located ∼100 kb upstream and downstream of "Stream" AyLTR were identified as "neighboring" genes.

Transcriptome Analysis of LTR-Retrotransposons in A. pernyi Larval Midgut Samples
To obtain a global view of LTR-retrotransposon expression of the A.pernyi larval midgut in response to ApNPV infection, the A. yamamai genome was first used as a reference genome to extract LTR-retrotransposon sequences from the raw transcriptome data of the study described in [3]. After (A) According to the positional relationship between AyLTRs and their neighboring genes, AyLTRs were classified into three types, "In", "Part" and "Stream". "In" represent AyLTRs that encompass gene exon sequences. "Part" represents AyLTRs that partially overlap with exon sequences. "Stream" represents independent AyLTRs that do not overlap with host gene exons. (B) Histogram description of Gene Ontology enrichment of neighboring genes of "Stream" AyLTRs. Genes were assigned to three categories: biological process (BP), cellular component (CC), and molecular function (MF). (C) The neighboring genes of "Stream" AyLTRs that were enriched in KEGG pathways. Genes located ∼100 kb upstream and downstream of "Stream" AyLTR were identified as "neighboring" genes.

Transcriptome Analysis of LTR-Retrotransposons in A. pernyi Larval Midgut Samples
To obtain a global view of LTR-retrotransposon expression of the A. pernyi larval midgut in response to ApNPV infection, the A. yamamai genome was first used as a reference genome to extract LTR-retrotransposon sequences from the raw transcriptome data of the study described in [3]. After discarding adaptor and low-quality reads, 71,599,652, 57,723,648, 70,342,118 and 61,942,502 clean reads were obtained from four cDNA libraries of A. pernyi midgut (two controls and two ApNPV-infected samples; SRR2919240, SRR2919241, SRR2919242 and SRR2919243) ( Table 1) [3]. The clean reads were mapped onto the A. yamamai reference genome, and the mapping rate of each library ranged from 63.31-67.39% (Table 1). The analysis resulted in the identification of 93 full-length LTR-retrotransposons that were transcribed in ApNPV-infected A. pernyi larval midgut samples or their controls ( Figure 4A). The expression of these LTR-retrotransposons in the four A. pernyi midgut samples are presented in the heatmap of Figure 4B. Compared to uninfected controls, most of LTR-retrotransposons were up-regulated in ApNPV-infected samples ( Figure 4B). However, a considerable number of LTR-retrotransposons exhibit variable expression between the two samples in each category, especially in the samples after infection with ApNPV. LTRs that expressed inconsistencies in expression between samples in each category were removed, and finally only the LTR-retrotransposons that are consistently down-or up-regulated in each sample were used for subsequent analysis. Details of the expression of these LTR-retrotransposons are presented in Supplementary Materials Table S3.
Furthermore, on the basis of the differential expression analysis, 12 significant DELs were identified between the ApNPV-infected and uninfected midgut samples, of which six were up-regulated and six were down-regulated ( Figure 4C, Figure 5A; Table 2). In addition, compared to the uninfected control, 2963 up-regulated DEGs and 816 down-regulated DEGs were identified in ApNPV-infected A. pernyi larval midgut samples ( Figure 5B).       To further explore the connection between DELs and DEGs, we examined whether adjacent putative cis-targets of DELs could correspond to DEGs during ApNPV infection. On the basis of this analysis, a DELs-DEGs interaction network was constructed using Cytoscape ( Figure 6). In the network, 9 DELs-DEGs connections were positively correlated, whereas another 14 connections were negatively correlated ( Figure 6). The specific description of these 23 DEGs is shown in Table 3. To further explore the connection between DELs and DEGs, we examined whether adjacent putative cis-targets of DELs could correspond to DEGs during ApNPV infection. On the basis of this analysis, a DELs-DEGs interaction network was constructed using Cytoscape ( Figure 6). In the network, 9 DELs-DEGs connections were positively correlated, whereas another 14 connections were negatively correlated ( Figure 6). The specific description of these 23 DEGs is shown in Table 3.
Interestingly, DEGs involved in Notch signaling pathway including epsin-2 and lysine-specific demethylase lid (lsdl) were adjacent to AY_261_133645_148170 and AY_225_415711_427236, respectively, and might be cis-regulated by these DELs (Figures 6 and 7). Additionally, ataxin-2, an SG-related DEG, was also identified as a putative cis-target of AY_545_152779_164768 (Figures 6 and  7).  Figure 6, epsin-2 (evm.TU.AY_261.18) and lysine-specific demethylase lid (XLOC_008622) are considered as putative cis-targets of AY_261_133645_148170 and AY_225_415711_427236, respectively, in A. pernyi. Therefore, we speculate that AY-261 and AY-225 might participate in the regulation of the Notch signaling pathway by regulating the expression of epsin-2 and lysine-specific demethylase lid by a mechanism in cis. (B) Stress granule formation can occur as a result of eIF2α phosphorylation caused by eIF2a kinase. eIF2a kinase expression was activated during ApNPV infection in A. pernyi. Ataxin- Figure 7. DEGs that are located in the neighborhood (putative cis-targets) of DELs and that are associated with Notch signaling pathway and stress granule formation. (A) Some DEGs including Numb, Deltex, CBP, Epsin-2 and lysine-specific demethylase lid (lsdl) that are associated with Notch signaling are increased in expression during ApNPV infection. As also shown in Figure 6, epsin-2 (evm.TU.AY_261.18) and lysine-specific demethylase lid (XLOC_008622) are considered as putative cis-targets of AY_261_133645_148170 and AY_225_415711_427236, respectively, in A. pernyi. Therefore, we speculate that AY-261 and AY-225 might participate in the regulation of the Notch signaling pathway by regulating the expression of epsin-2 and lysine-specific demethylase lid by a mechanism in cis. Interestingly, DEGs involved in Notch signaling pathway including epsin-2 and lysine-specific demethylase lid (lsdl) were adjacent to AY_261_133645_148170 and AY_225_415711_427236, respectively, and might be cis-regulated by these DELs (Figures 6 and 7). Additionally, ataxin-2, an SG-related DEG, was also identified as a putative cis-target of AY_545_152779_164768 (Figures 6 and 7).

Discussion
In the present study, we attempted to explore the role of LTR-retrotransposons in the process of ApNPV infection in the wild silkmoth A. pernyi. However, genome-wide detection of LTR-retrotransposons in A. pernyi has not been performed due to the deficiency of genomic data.
Fortunately, the annotated genome sequence of A. yamamai, the first published genome within the Saturniidae family, was released recently [15]. A. pernyi and A. yamamai are sibling species within the Antheraea genus and the Saturniidae family [27]. Given that no A. pernyi genome is currently available, the genome of A. yamamai is used as reference genome in this study. Indeed, the clean reads of A. pernyi transcriptome could be mapped confidently onto the A. yamamai reference genome. The mapping rate of each library ranged from 63.31-67.39%.
For our analysis, a library of LTR-retrotransposons of wild silkworm (Antherea) was first constructed after screening of the A. yamamai genome with LTRharvest and LTRdigest. The obtained library contained 541 full-length AyLTRs and 22,666 solo AyLTRs. Based on a comparative analysis of the conserved RT domain, AyLTRs could be grouped into the canonical phylogenetic clades of Gypsy, Copia, and Belpao [28].
However, among 541 full-length AyLTRs, only eleven AyLTRs contained the env element. The infectious ability of insect LTR-retrotransposons is thought to be associated with the expression of the envelope protein encoded by the env gene [29]. Thus far, only the well-known gypsy and ZAM elements of D. melanogaster have been shown to possess infectious properties [30,31]. Furthermore, it was reported that the env gene of LTR-retrotransposons was evolutionary related with the envelope fusion protein lineage of baculovirus [32,33]. The region of the envelope fusion protein with the highest sequence similarity with insect LTR-retrotransposons includes the furin hydrolysis signal and a fusion peptide located downstream [33]. In our study, we observed that the envelope fusion protein of AY-476 and AY-37 possessed a complete furin cleavage site with a downstream fusion peptide, whereas in AY-58 and AY87 only a partial furin cleavage site was detected. These results suggest that particular AyLTRs with env gene elements, such as AY-476 and AY-37, might possess infectious properties. It is believed that some LTR-retrotransposons which lacked the env gene, became integrated into the dsDNA genome of a baculovirus from which they could "capture" the env gene [32]. However, we found that the env sequences of AyLTRs only had a distant relationship with genes encoding F (Fa and Fb) and GP64 protein from Group I and Group II NPVs in general, and more specifically from ApNPV and AyNPV. These results indicate that the env elements were not derived from the envelope gene of ApNPV and AyNPV.
It is well documented that LTR-retrotransposons can affect the physiopathology of host cells at multiple levels. For instance, LTR-retrotransposons can modulate the expression of adjacent host genes in the human genome [34]. While the expression of LTR-retrotransposon proteins with conventional retroviral functions can influence the host's physiological or pathological states [35], it was also reported that non-coding LTR-retrotransposons can be biologically active [36]. To investigate interactions between LTR-transposons and cellular genes in cis, we first identified the host genes that are adjacent to AyLTRs within the A. yamamai genome. It was assumed that some of these AyLTRs may be involved in host biological responses through regulation of the expression of host genes.
Based on the information on AyLTRs and their neighboring genes, it was attempted to decipher the biological function of AyLTR during ApNPV infection of A. pernyi by expression analysis. Due to the absence of genomic data of A. pernyi, the published transcriptome analysis of the A. pernyi infected with ApNPV was de novo assembled using the Trinity platform in the published study [3]. In this study, on the other hand, the transcriptome of LTR-retrotransposons and genes was analyzed after mapping of the reads to the A. yamamai genome. Using this approach, 93 full-length LTR-retrotransposons were found to be transcribed in ApNPV-infected A. pernyi larval midgut samples or their uninfected controls. Interestingly, six LTR-retrotransposons were significantly up-regulated, and six LTR-retrotransposons were significantly down-regulated in ApNPV-infected A. pernyi midgut. Whether these 12 DELs can play a functional role during ApNPV infection of A. pernyi will require further experimentation. Moreover, interactions between 7 DELs and their 23 adjacent host DEGs were identified during ApNPV infection. We speculate that the expression of the 23 identified DEGs is under the regulatory control of the seven DELs by a mechanism of interaction in cis along the chromosome during ApNPV infection in A. pernyi. Indeed, an involvement in host-virus interaction as well as immune response processes is very well documented in the case of endogenous retroviruses, a special category of LTR-retrotransposons [37][38][39].
For a few of these 23 DEGs, a role in the response of the virus-infected host was observed in the silkworm Bombyx mori. The expression of reverse transcriptase (RT) was up-regulated during BmNPV infection in BmNPV resistant silkworm strains [40]. The up-regulation of RT in A. pernyi midgut during ApNPV infection confirms its possible role in host antiviral response. Conversely, the expression of Tret1 was down-regulated in A. pernyi midgut during ApNPV infection, which fits with previous research documenting that down-regulation of trehalose transporter Tret1-like might enhance infection of BmNPV in B. mori [40].
For a more in-depth analysis of a subset of our data, attention was focused on a possible role of the Notch signaling pathway and stress granule regulation during ApNPV infection of the midgut of A. pernyi and the possible involvement of their regulation by LTR-retrotransposons.
The Notch signaling pathway is important in development, tissue homeostasis, as well as disease [41], including viral disease [42,43]. In this study, several up-regulated DEGs were found to be associated with Notch signaling during ApNPV infection ( Figure 7A). Among these DEGs, the homolog of deltex 1 is involved in T cell immunity in mammals [44]. CREB-binding protein (CBP) encodes a coactivator protein (histone acetyltransferase or HAT) that plays a role in the innate antiviral immunity pathway in vertebrates [45] and is also targeted by virus-encoded suppressor mechanisms [46]. In the LNX/Numb/Notch pathway, Numb protein (DEG during ApNPV infection) was also identified as a target of the Np9 protein encoded by the human endogenous retrovirus K [47]. Thus, we speculate that Notch signaling-related DEGs such as deltex, CBP and Numb might be associated with host responses in A. pernyi against ApNPV.
In addition, epsin-2 and lysine-specific demethylase lid were up-regulated after ApNPV infection in our analysis (Table 3). Coincidentally, these two genes were found to play an important role in Notch signaling [48,49]. Moreover, epsin-2 and lysine-specific demethylase lid are putative cis-targets of AY_261_133645_148170 and AY_225_415711_427236, respectively in A. pernyi. Based on these interesting connections, we forward the hypothesis that expression of important genes in Notch signaling in the host can be regulated by LTR-retrotransposons, which in turn are responsive to ApNPV infection.
In mammalian systems, ataxin-2-like protein encodes a component of stress granules (SG) that also regulates p-body (PB) formation [50]. SGs and PBs provide cell homeostasis and mRNA stability during the stress response induced by viral infection [51]. More importantly, SGs play a critical role in the host antiviral immune response. SG formation can interfere with viral replication because all viruses require the host translation machinery to synthesize viral proteins. Accordingly, many viruses have evolved diverse mechanisms to inhibit SG formation. Antiviral functions of SGs include the establishment of an antiviral state by limiting viral protein accumulation and the regulation of signaling cascades that affect virus replication and immune responses [52]. In mammals, SG formation can occur as a result of eIF2α phosphorylation caused by diverse eIF2a kinases activated by different stress conditions [51]. Interestingly, the expression of eIF2a kinase was also up-regulated in A. pernyi after ApNPV infection ( Figure 7C). Accordingly, we speculate that differential expression of ataxin-2-like protein and eIF2a kinase reflects the regulation of SG formation during ApNPV infection in A. pernyi ( Figure 7B). Moreover, given that ataxin-2 has been identified as a putative cis-target of AY_545_152779_164768 (Table 3), the possibility is considered that AY_545 is involved in the regulation of stress granule formation during ApNPV infection.
As summarized in Figure 8, a wild silkworm LTR-retrotransposons library was established based on the first draft genome in the family of Saturniidae, which included 22,666 solo LTRs and 541 full-length LTRs. Using the A. yamamai genome as a reference genome [15] and published raw data of the transcriptome of ApNPV-infected A. pernyi midgut [3], 93 full-length LTR-retrotransposons were found to be transcribed in ApNPV-infected A. pernyi larval midgut samples or their uninfected controls. Candidate DEGs were identified, including epsin-2, lysine-specific demethylase lid and ataxin-2-like protein that are involved in Notch signaling and stress granule formation, that could be subject to regulation in cis by DELs in A. pernyi during ApNPV infection. Further experimentation is required to verify whether LTR-retrotransposons are involved in the response of A. pernyi to ApNPV by regulating particular genes associated with Notch signaling and stress granule formation. The functional relevance of other DEL-DEG combinations identified by bioinformatics likewise needs experimental validation in future studies.
As summarized in Figure 8, a wild silkworm LTR-retrotransposons library was established based on the first draft genome in the family of Saturniidae, which included 22,666 solo LTRs and 541 full-length LTRs. Using the A. yamamai genome as a reference genome [15] and published raw data of the transcriptome of ApNPV-infected A. pernyi midgut [3], 93 full-length LTRretrotransposons were found to be transcribed in ApNPV-infected A. pernyi larval midgut samples or their uninfected controls. Candidate DEGs were identified, including epsin-2, lysine-specific demethylase lid and ataxin-2-like protein that are involved in Notch signaling and stress granule formation, that could be subject to regulation in cis by DELs in A. pernyi during ApNPV infection. Further experimentation is required to verify whether LTR-retrotransposons are involved in the response of A. pernyi to ApNPV by regulating particular genes associated with Notch signaling and stress granule formation. The functional relevance of other DEL-DEG combinations identified by bioinformatics likewise needs experimental validation in future studies.  [15], which is the first draft genome in the family of saturniid moths, a wild silkworm LTRretrotransposons library was constructed of 22,666 solo LTRs and 541 full-length LTRs. Using this information, published RNA sequence raw data of ApNPV-infected A. pernyi larval midgut and uninfected control samples [3] were analyzed to identify 12 DELs and 3779 DEGs with respect to ApNPV infection. In more detailed analysis, several DEGs, that are associated with Notch signaling and stress granule formation, are considered to be regulated in cis by DELs during ApNPV infection (Dotted line). Our study indicates a potential role for LTR-retrotransposons to regulate the host gene response during ApNPV infection of A. pernyi.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Table S1: Properties of identified AyLTRs. Table S2: GO and KEGG pathway analysis of neighboring host genes (potential cis-targets) of "Stream" AyLTRs. Table S3: Expression data of LTR-retrotransposon in ApNPV-infected A. pernyi larval midgut samples and their uninfected controls.
Author contributions: M.F. participated in the design of the study, collected and analyzed data and drafted the manuscript. F.R., Y.Z., N.Z. and Q.L. helped with the data analysis. L.S. and J.S. participated in the design and coordination of the study and drafted the manuscript. All authors read and approved the final manuscript. Figure 8. General overview of the study and its conclusions. Based on the genome of Antheraea yamamai [15], which is the first draft genome in the family of saturniid moths, a wild silkworm LTR-retrotransposons library was constructed of 22,666 solo LTRs and 541 full-length LTRs. Using this information, published RNA sequence raw data of ApNPV-infected A. pernyi larval midgut and uninfected control samples [3] were analyzed to identify 12 DELs and 3779 DEGs with respect to ApNPV infection. In more detailed analysis, several DEGs, that are associated with Notch signaling and stress granule formation, are considered to be regulated in cis by DELs during ApNPV infection (Dotted line). Our study indicates a potential role for LTR-retrotransposons to regulate the host gene response during ApNPV infection of A. pernyi.
Author Contributions: M.F. participated in the design of the study, collected and analyzed data and drafted the manuscript. F.R., Y.Z., N.Z. and Q.L. helped with the data analysis. L.S. and J.S. participated in the design and coordination of the study and drafted the manuscript. All authors read and approved the final manuscript.