Comparison of Transcriptomic Signatures between Monkeypox-Infected Monkey and Human Cell Lines

Monkeypox virus (MPV) is a smallpox-like virus belonging to the genus Orthopoxvirus of the family Poxviridae. Unlike smallpox with no animal reservoir identified and patients suffering from milder symptoms with less mortality, several animals were confirmed to serve as natural hosts of MPV. The reemergence of a recently reported monkeypox epidemic outbreak in nonendemic countries has raised concerns about a global outburst. Since the underlying mechanism of animal-to-human transmission remains largely unknown, comprehensive analyses to discover principal differences in gene signatures during disease progression have become ever more critical. In this study, two MPV-infected in vitro models, including human immortal epithelial cancer (HeLa) cells and rhesus monkey (Macaca mulatta) kidney epithelial (MK2) cells, were chosen as the two subjects to identify alterations in gene expression profiles, together with co-regulated genes and pathways that are affected during monkeypox disease progression. Using Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and MetaCore analyses, we discovered that elevated expression of genes associated with interleukins (ILs), G protein-coupled receptors (GPCRs), heat shock proteins (HSPs), Toll-like receptors (TLRs), and metabolic-related pathways play major roles in disease progression of both monkeypox-infected monkey MK2 and human HeLa cell lines. Interestingly, our analytical results also revealed that a cluster of differentiation 40 (CD40), plasmin, and histamine served as major regulators in the monkeypox-infected monkey MK2 cell line model, while interferons (IFNs), macrophages, and neutrophil-related signaling pathways dominated the monkeypox-infected human HeLa cell line model. Among immune pathways of interest, apart from traditional monkeypox-regulated signaling pathways such as nuclear factor- (NF-κB), mitogen-activated protein kinases (MAPKs), and tumor necrosis factors (TNFs), we also identified highly significantly expressed genes in both monkey and human models that played pivotal roles during the progression of monkeypox infection, including CXCL1, TNFAIP3, BIRC3, IL6, CCL2, ZC3H12A, IL11, CSF2, LIF, PTX3, IER3, EGR1, ADORA2A, and DUOX1, together with several epigenetic regulators, such as histone cluster family gene members, HIST1H3D, HIST1H2BJ, etc. These findings might contribute to specific underlying mechanisms related to the pathophysiology and provide suggestions regarding modes of transmission, post-infectious sequelae, and vaccine development for monkeypox in the future.


Introduction
Monkeypox virus (MPV/MPXV) is commonly known as one of the rare systemic infections caused by one specie of the genus Orthopoxvirus in the family Poxviridae [1,2]. Other than the smallpox virus (variola virus (VARV)), the best-known species of this genus, monkeypox virus, together with cowpox virus (CPXV), horsepox virus (HSPV), raccoonpox virus (RCNV), camelpox virus (CMLV), mousepox virus (also referred to as Ectromelia virus (ECTV)), and Alaskapox (AKPV) virus, are species believed to infect healthy humans through different modes of zoonotic transmission [3][4][5]. After the eradication of smallpox, a result of universal vaccination, almost no natural outbreaks related to Orthopoxvirus have been recorded until recently, after a significant number of confirmed monkeypox cases were reported worldwide [6]. The first discovery of this virus in monkeys in the lab (in 1958), followed by the first human infection detected in Africa (in 1970), meant that monkeypox was no longer endemic in animals as previously assumed; at least 400 human cases scattered across the western and central parts of the African continent have been diagnosed with similar clinical characteristics to smallpox but causing milder symptoms and being less lethal [7][8][9]. As of May 2022, multiple cases of MPV have been documented in at least 11 regions outside of Africa, including the European Zone, Britain, Middle East, Australia, the North American continent (including the United States, Canada, and South America), resulting in unprecedented outbreaks in those developed countries [8][9][10]. Further, public health investigations revealed the approximately 20-fold increase in incidences that resulted from the discontinuation of smallpox vaccination since 1980, while the probable source of infection could be either an animal reservoir (mainly rodents) or cross-infection [11,12]. Transmission through person-to-person contact is claimed to primarily occur through close and prolonged contact with fluids from skin lesions, respiratory droplets, fomites of infected people, and even during sexual intercourse among populations of men who have sex with men (MSM) [13,14].
To date, definite diagnoses of monkeypox are largely based on either a polymerase chain reaction (PCR) of skin lesion specimens or immunohistochemical (IHC) staining of biopsies taken from those infected, as misdiagnoses of other blister-like diseases such as smallpox and chickenpox may occur in cases where only a visual skin examination and/or dermoscopy are solely employed [15,16]. Although the genomes of several variants of the MPV have been fully sequenced, there are limited reports on the underlying mechanisms of viral transmission and pathogenic pathways of this disease [17,18]. Moreover, proven treatment and pre-ventive vaccines specifically designed for monkeypox have not yet been formally approved by the agencies, thus clinical care only focuses on symptomatic treatment. Therefore, developing a research design that focuses on the underlying molecular mechanisms during disease progression would greatly contribute to improving the situation as it provides necessary background knowledge for subsequent studies, such as diagnostic test development, epidemiological tracking, designs of vaccines and antiviral drugs, and other related adjuvant treatments.
In this study, we attempted to provide a concise and informative overview of principal differences in gene expression profiles between two MPV-infected in vitro models, including human immortal epithelial cancer (HeLa) cells and rhesus monkey (Macaca mulatta) kidney epithelial (MK2) cells, together with co-regulated genes and pathways that are affected during disease progression. The Gene Expression Omnibus (GEO) database and bioinformatics approaches were employed to perform analyses of available samples related to these two models ( Figure 1). We also revealed several critical regulatory downstream networks and investigated the functions of potential genes to determine whether they can serve as clinical biomarkers for MPV. These findings may contribute to an understanding of specific underlying mechanisms related to the pathophysiology and provide suggestions regarding modes of transmission, post-infection sequelae, and vaccine development for monkeypox in the future.

Bioinformatics Application for Data Acquisition and
Processing. In this study, we investigated an in vitro system using gene expression data of two separate subjects, including human immortal epithelial (HeLa) cells and monkey Macaca mulatta kidney epithelial (MK2) cells to characterize the effects of physical monkeypox viral infection. Relevant microarray gene expression datasets from the GEO (http://www.ncbi.nlm.nih.gov/geo)-an international public functional genomics data repository affiliated with the National Center for Biotechnology Information (NCBI; https://www.ncbi.nlm.nih.gov/)-were retrieved for subsequent analyses. In particular, MPV-induced alterations in transcriptomic profiles of 7-hour monkeypox-infected Macaca mulatta kidney epithelial (MK2) cells were obtained from the GSE21001 dataset [19], whereas similar information of 6-hour monkeypox-infected human immortal epithelial (HeLa) cells was acquired from the GSE36854 dataset [20]. The following analyses of both studies are based on comparison of each subject with its corresponding mockinfected counterpart as control.

2
Journal of Immunology Research Before being mapped onto Ensembl features using packages attached to the biomaRt tool (vers. 2.26.1), the CLC Genomics Workbench vers. 10.1 (CLC bio, Aarhus, Denmark) was utilized to process and standardize the raw data [21][22][23]. Signals were subsequently presented as a clustered heatmap based on messenger (m)RNA expression patterns using pheatmap (vers. 1.0.12) in the R environment [24][25][26], as well as the online platform available at http:// www.bioinformatics.com.cn/srplot. For subsequent discoveries of enriched functionally related gene groups, the Database for Annotation, Visualization, and Integrated Discovery (DAVID), which is maintained by the Laboratory of Human Retrovirology and Immunoinformatic (LHRI) [27][28][29], was chosen as a powerful tool by the various functions it offers, especially in term of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). As previously described [30][31][32], the top 5% of upregulated genes derived from monkeypoxinfected groups were further processed as input data for the GO analyses [33][34][35], while a p value of <0.05 was set as the significance level. After biological networks and related pathways were constructed, a cut-off value of 0.05 was applied to select significantly enriched pathways or groups of annotated genes.

Pathway-Based Analyses and Network Enrichment
Analyses. MetaCore (Enrichment Analysis Workflow and analysis network; GeneGo, St. Joseph, MN, USA) was designed to identify biological processes related to gene microarray data. The top 5% of upregulated genes with substantial differences in transcriptome levels were processed as input to the MetaCore program to compare average levels of gene expressions in the two models. Signal transduction pathways were examined. A statistically significant difference was indicated by a p value of <0.05.

Construction of a Protein-Protein-Interacting (PPI)
Network using the STRING Analysis. The STRING database (vers. 11.0) (https://string-db.org/) serves as a search engine and a protein information resource that provides a considerable number of proteins and known interaction activities upon which differentially expressed genes (DEGs) were examined [36][37][38]. The Cytoscape stringApp was also employed in this work to construct a PPI-network, while

GO Analysis of Monkeypox-Infected Models.
For the monkey cell line model, transcriptomic data of the 7-hour monkeypox-infected Macaca mulatta kidney epithelial (MK2) cell retrieved from the GSE21001 dataset were matched with a mock treatment group for a comparative analysis of DEGs. GO analytical results associated with monkeypox infection are depicted in Figure 2(a). Pathways assigned to cluster 1 (c1) contained the most significant highly expressed genes and pathways in the monkeypoxinfected group, whereas cluster 2 (c2) exerted an opposite trend. The GO analysis suggested that c1 was related to inflammatory responses, chemokine activities, immune system processes, responses to peptidoglycans, cellular responses to lipopolysaccharide, nucleosome assembly, chromatin assembly or disassembly, responses to ionizing radiation, and responses to temperature stimuli. Pathways within c2 were associated with kinase binding, thymus development, responses to endoplasmic reticular stress, transfer (t)RNA aminoacylation for protein translation, and responses to amino acid stimulation. For the 6-hour monkeypox-infected human cell line model, transcriptomic data of monkeypox-infected human HeLa cells were retrieved from the GSE36854 dataset. Pathways assigned to cluster 1 (c1) were those most highly correlated with monkeypox infection, whereas cluster 2 (c2) exerted an opposite trend. The GO analysis suggested that pathways within c1 were related to growth factor activities, leukocyte activation, responses to organic cyclic compounds, inflammatory responses, regulation of the extracellular signal-regulated kinase 1 (ERK1) and the extracellular signal-regulated kinase 2 (ERK2) cascade, chemokine activities, and regulation of nitric oxide biosynthetic processes. Pathways in c2 were related to nucleosome assembly, protein heterodimerization activities, telomere maintenance, extracellular regions, the actin cytoskeleton, responses to protozoans, and coagulation.

Pathway Analysis of the Two Monkeypox-Infected
Models. Venn diagrams were employed to filter upregulated genes in common between the two models mentioned above. The top 5% of messenger (m)RNAs within the monkeypoxinfected Macaca mulatta kidney epithelial (MK2) cell model and the human immortal epithelial cancer (HeLa) counterpart were highly expressed compared to those within the mock control (Figure 3(a)). In total, 149 overlapping genes that were highly expressed in both monkey and human models were subsequently input to the MetaCore tool for a map analysis. Several standard maps were found to be related to monkeypox infection, including the "immune response_IL-1 signaling pathway," "NF-κB pathway in multiple myeloma," "immune response_CD40 signaling in B cells," "TNF-alpha-induced inflammation," "signaling in normal and asthmatic airway epithelium," "immune response_plasmin signaling," "glomerular injury in lupus nephritis," "immune response_IL-33 signaling pathway," "macrophage and dendritic cell phenotype shift in cancer," "renal tubulointerstitial injury in lupus nephritis," "signal transduction_nonapoptotic FasR (CD95) signaling," and "inflammatory mechanisms of pancreatic cancerogenesis" (Figure 3(b), Supplementary Table A1). Meanwhile, Endogenous Metabolic Networks corresponding to genes from the shortlist of upregulated genes shared by the two models revealed that pathways were associated with phospholipid biosynthesis, including the "Ceramide pathway," "N-acyl-sphingosine phosphate pathway," "2-arachidonoyl-glycerol 3-phosphocholine pathway," "phosphatidylethanolamine pathway," "(L)-valine pathways and transport," "carbohydrate metabolism_TCA and tricarboxylic acid transport," and "vitamin, mediator and cofactor metabolism_alpha-tocotrienol," which are metabolic signalling pathways dominated by these genes within both aforementioned models (Figures 3(c)-3(e)), Supplementary Table A2).
For the human immortal epithelial cancer (HeLa) cell model, the top 5% of upregulated genes compared to the mock treatment groups were input into the MetaCore platform for a comprehensive map analysis. The "immune response_IL-1 signaling pathway," "inflammatory mechanisms of pancreatic cancerogenesis," "G protein-coupled receptors signaling in lung cancer," "immune response_     Journal of Immunology Research plasmin signaling," "immune response_CD40 signaling in dendritic cells, monocytes, and macrophages," "interleukinsinduced inflammatory response in asthmatic airway fibroblasts," "TNF-alpha-induced inflammatory signaling in normal and asthmatic airway epithelium," "NF-κB-, AP-1-, and MAPKs-mediated proinflammatory cytokine production by eosinophils in asthma," "immune response_CD40 signaling in B cells," and "immune response_IL-17 signaling pathways" were among the most significant maps detected (Supplementary Figure A3, Table A5). Through the Metabolic Networks (Endogenous) corresponding to genes from the shortlist of upregulated genes shared by the monkeypox-infected human immortal epithelial cancer (HeLa) cells, data revealed that pathways associated with phospholipid biosynthesis included the "phosphatidylcholine pathway," "1-docosahexaenoyl-glycerol_3-phosphocholine pathway," "lysophosphatidic acid pathway," "1-icosatrienoylsn-glycero-3-phosphocholine pathway," "2-arachidonoyl-glycerol_3-phosphocholine pathway," and "GalNAcbeta1-3Gal pathway," which are metabolic signaling pathways dominated by these genes within this model (Supplementary Figure A4, Table A6). In addition to the pathways and network analysis for leading-edge genes generated by MetaCore as mentioned before, we also used the KEGG database to validate these data. The top enriched KEGG pathways were analyzed for common genes shared by the monkeypox-infected Macaca mulatta kidney epithelial (MK2) cell model and the human immortal epithelial cancer (HeLa) cell model, and pathways were ranked by p values (Supplementary Figure A5).  Journal of Immunology Research
Upregulated genes shared by these two models were also imported into the GO platform for cellular component annotations under monkeypox infection circumstances. In the case of the monkeypox-infected Macaca mulatta kidney epithelial (MK2) cell model, these genes mostly functioned in apical plasma membranes, neuronal cell bodies, the cytosol, nucleoplasm cytoplasm, extracellular exosomes, an integral component of presynaptic membranes, postsynaptic membranes, cytoplasmic exosomes, and the perikaryon ( Figure 5(a)). The same analytical flowchart was applied to the monkeypox-infected human immortal epithelial cancer (HeLa) cell model to identify cellular components under monkeypox infection circumstances. Genes with the greatest significance functioned in the extracellular space, receptor complexes, integral component of membranes, caveolae, external side of plasma membranes, extracellular matrices, voltage-gated calcium channel complexes, presynaptic active zones, glutamatergic synapses, and dendritic shafts ( Figure 5(b)). From such observations gained from the two in vitro monkeypox-infected models, we concluded that certain cellular mechanisms of monkeypox infection are almost identical but slightly different when comparing the two species ( Figure 5(c)).

Discussion
In most cases, viral infections only occur when a virus that causes disease successfully invades a host, beginning by penetrating the antiviral defenses, inserting its genetic materials, taking control of host cells' machinery to rapidly create a huge number of viral replications, and ultimately becoming transmissible between individuals. In the context of the world having gone through many epidemic zoonotic outbreaks that originated from animals [39][40][41], including the human immunodeficiency virus (HIV; which jumped from chimpanzees), influenza (from birds and pigs), bovine spongiform encephalopathy (BSE; from cows), Ebola (from bats), and at least three zoonotic coronavirus diseases in chronological order, including the severe acute respiratory syndrome coronavirus (SARS-CoV), the Middle Eastern respiratory syndrome coronavirus (MERS-CoV), and the severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2). Therefore, gaining a comprehensive understanding of the mechanisms of infection in both human and nonhuman hosts is urgently needed [42][43][44][45].
In attempts to prevent further outbreaks and deaths from these diseases, worldwide efforts to discover rapid early-detection methods, effective antiviral drugs, and preventive vaccines have been increasing over the years [46,47]. Other than conventional approaches, novel and advantageous technologies, such as high-throughput sequencing, searchable database repositories of gene expression data, and assistance from web-based or computer-based bioinformatics tools, have been employed to elucidate underlying molecular mechanisms during viral infection and explore how different host cells respond via performing analyses of alterations in gene expressions [48][49][50][51][52][53][54].
In the case of monkeypox infection, according to previous studies, a considerable number of immunomodulationrelated cytokines, including interleukin-1 (IL-1), tumor necrosis factors (TNFs), and interferons (IFNs), were reported to be involved in disease progression. The ILs superfamily (including IL-1, IL-18, and IL-33) and TNF signaling-mediated immune responses are important for antiviral activity [55,56]; however, the roles of these signaling pathways in monkeypox-infected somatic cells remain unclear. Interestingly, in case of other zoonotic viral diseases mentioned before, numerous cellular and molecular immunomodulatory pathways were also proven to play crucial roles. In particular, the G protein-coupled receptors (GPCRs) not only helps stimulate internalization during infection with the influenza A virus [57] but also takes part in Epstein-Barr virus-(EBV-) mediated immunosuppression and oncogenesis [58]. Studies on the molecular chaperonin HSP60 confirmed its interactions with Ebola virus infection [59] and human hepatitis B virus polymerase [60], while HSP70 plays a crucial role in mediating the progression of Zika virus in infected host cells [61]. Histamine, on the other hand, contributes to severe pneumonia in pigs infected with the H1N1 influenza virus [62], whereas levels of histamine and leukotrienes presented in acute dengue patients are closely related to disease severity [63].
Our research findings were consistent with previous studies on monkeypox infection as mentioned above; the GO and pathway analyses revealed similar enriched pathways for both monkeypox-infected monkey and human models, including IL-1, IL-33, IL-17, IL-18, NF-κB, MAPKs, and TNF-R2 signaling (Figure 3). However, interestingly, we also discovered relationships between monkeypox and GPCRs, HSP60/70, histamine, plasmin, and histone cluster-related signaling (Supplementary Figures A1, A2) that have rarely been reported before.
Regarding 24 common genes that were highly expressed in both in vitro models (Figure 4), when compared to the previous literature, similar expression patterns were also observed in case of CXCL1, TNFAIP3, BIRC3, IL6, CCL2, ZC3H12A, IL11, CSF2, LIF, PTX3, IER3, EGR1, ADORA2A, DUOX1, and several histone cluster family members such as HIST1H3D and HIST1H2BJ. According to a recent study, there were dramatic increases in serum concentrations of IL1-RA, IL-6, and IFN-γ in cynomolgus macaques after being infected with aerosolized MPV [64]. In myeloid cells, 10 Journal of Immunology Research   Journal of Immunology Research TNFAIP3 deficiency helps protect host cells against invasion by the influenza A virus [65]. Both ZC3H12A/MCPIP1 and CDKN1A/p21 may contribute to the synergistic effect on replication control and pharmacological manipulation of HIV-1 [66]. Infection with SARS-CoV-2 and several influenza viruses can lead to alterations in CSF2 and IL-11 expressions [67][68][69]. PTX3, the expression of which was reported to be upregulated by TNF-α during acute lung injury, was recently found to serve as a reliable prognostic indicator in predicting short-term mortality from SARS-CoV-2 [70]. Apoptosis and impacted viral replication were found to be induced by infection with the Venezuelan equine encephalitis virus; however, knockdown of early growth response 1 (EGR1) significantly inhibited these signaling pathways [71]. Two members of the NADPH oxidase (NOX) family, including the DUOX1 and DUOX2 proteins, are potential therapeutic targets against infectious diseases caused by the influenza A virus [72]. Previous research demonstrated that DNA viruses exploit host cellular epigenetic processes to their advantage during infection [73,74]. Most interestingly, several histone cluster family members were reported to be involved in differential responses of human fetal brain neural stem cells to Zika virus infection [75]. For example, three members of the histone cluster 1 H2B

13
Journal of Immunology Research (HIST1H2B) family, namely HIST1H2BB, HIST1H2BK, and HIST1H2BO, were among those hub genes and pathways that contributed to Zika virus infection [76]. However, only a limited number of studies have described such relationships in terms of monkeypox. Interestingly, results of our bioinformatics analysis revealed that expression level of certain histone cluster family members, such as HIST1H3D, HIST1H2BJ, HIST1H2AK, HIST1H2AD, HIST1H2AC, HIST1H1B, HIST2H2AB, HIST1H2BM, HIST1H2BH, and HIST1H2BB, were elevated in both monkeypox-infected cell models of monkey and human. These data suggested that some of these epigenetic regulators may also play important roles in monkeypox infection (Figure 4). Cross-validation performed in the 6-hour monkeypox-infected human immortal epithelial cancer (HeLa) cell model (GSE24125) [77] revealed similar patterns of elevated expressions of the aforementioned genes (Supplementary Figure A6). As the most challenging aspect of diagnosing monkeypox infection is distinguishing monkeypox vesicular rash eruption from those of other diseases [78], our study provides a prospect of the histone cluster family as potential regulators correlated with monkeypox infection. Since the current study only focused on signaling pathways regulated by DEGs caused by monkeypox infection, further investigation into the underlying pathogenesis is required for consolidated guidelines for monkeypox treatment in the future.

Conclusion
In the current study, we focused on investigating genetic signatures and their related pathways in two different monkeypox-infected models. Our research findings revealed the critical roles of multiple genes and their regulatory pathways in both monkeypox-infected models, which not only directly in line with previous literature but also provided novel functions as well as highly expressed genes related to MPV infection that had been less well revealed in the past. In the context of widespread epidemics of viral infectious diseases that are prevalent today, this work may contribute to bridging the traditional gap between bench research and clinical applications. However, with attempts to validate and examine the potentials of these genes as novel antiviral therapies, further expanding research on disease basics may provide more-concrete perspectives.

Data Availability
The datasets used and analyzed during the current study are available from the corresponding authors on reasonable request. Kaohsiung Medical University Hospital (KMUH110-0T06, KMUH109-9R82 and KMUH110-0M75). This work was financially supported by the Higher Education Sprout Project of the Ministry of Education (MOE) in Taiwan. The authors acknowledge the online platform for data analysis and visualization (http://www.bioinformatics.com.cn/), and we also thank the statistical/computational/technical support of the Clinical Data Center, Office of Data Science, Taipei Medical University, Taiwan. Figure S1: MetaCore-generated pathways regulated by top upregulated genes in the monkeypox-infected Macaca mulatta kidney epithelial (MK2) cell model. Figure S2: MetaCore-generated Endogenous Metabolic Network among top upregulated genes in the monkeypox-infected Macaca mulatta kidney epithelial (MK2) cells, performed on the GSE21001 dataset. Figure S3: MetaCore-generated pathways related to top upregulated genes in the human immortal epithelial cancer (HeLa) cell model. Figure S4: MetaCore-generated Endogenous Metabolic Network among top upregulated genes in the human immortal epithelial cancer (HeLa) cell model, performed on the GSE36854 dataset. Figure S5: Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathways in Formula A analysis in different models. Figure S6: heatmap displaying specific differentially regulated genes within the monkeypoxinfected human immortal epithelial cancer (HeLa) cell model (GSE24125) versus the mock-infected group. Table S1: top 50 upregulated pathways in the MetaCore analysis of potential maps. Comparison between monkeypox-infected Macaca mulatta kidney epithelial (MK2) cells and human epithelial (HeLa) cells from the GSE21001 and GSE36854 datasets. Table S2: top 50 upregulated pathways in the MetaCore analysis of potential Endogenous Metabolic Networks. Comparison between monkeypox-infected Macaca mulatta kidney epithelial (MK2) cells and human epithelial (HeLa) cells from the GSE21001 and GSE36854 datasets. Table S3: top 50 upregulated pathways in the MetaCore analysis of potential maps. Comparison between monkeypox-infected and mock-infectedMacaca mulatta kidney epithelial (MK2) cells from the GSE21001 dataset.