Integrative Multiomics and Regulatory Network Analyses Uncovers the Role of OAS3, TRAFD1, miR-222-3p, and miR-125b-5p in Hepatitis E Virus Infection

The hepatitis E virus (HEV) is a long-ignored virus that has spread globally with time. It ranked 6th among the top risk-ranking viruses with high zoonotic spillover potential; thus, considering its viral threats is a pressing priority. The molecular pathophysiology of HEV infection or the underlying cause is limited. Therefore, we incorporated an unbiased, systematic methodology to get insights into the biological heterogeneity associated with the HEV. Our study fetched 93 and 2016 differentially expressed genes (DEGs) from chronic HEV (CHEV) infection in kidney-transplant patients, followed by hub module selection from a weighted gene co-expression network (WGCN). Most of the hub genes identified in this study were associated with interferon (IFN) signaling pathways. Amongst the genes induced by IFNs, the 2′-5′-oligoadenylate synthase 3 (OAS3) protein was upregulated. Protein-protein interaction (PPI) modular, functional enrichment, and feed-forward loop (FFL) analyses led to the identification of two key miRNAs, i.e., miR-222-3p and miR-125b-5p, which showed a strong association with the OAS3 gene and TRAF-type zinc finger domain containing 1 (TRAFD1) transcription factor (TF) based on essential centrality measures. Further experimental studies are required to substantiate the significance of these FFL-associated genes and miRNAs with their respective functions in CHEV. To our knowledge, it is the first time that miR-222-3p has been described as a reference miRNA for use in CHEV sample analyses. In conclusion, our study has enlightened a few budding targets of HEV, which might help us understand the cellular and molecular pathways dysregulated in HEV through various factors. Thus, providing a novel insight into its pathophysiology and progression dynamics.


Introduction
Animal-sourced viruses have posed a severe threat to public health globally, as seen by the recent appearance and spread of zoonotic viruses such as the Ebola virus and the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2). The Hepatitis E virus (HEV) ranks among the top 10 risk-ranking viruses with a high potential for zoonotic spillover. Therefore, considering its viral risks is a growing concern [1]. The current status of HEV corresponds to 1/8th of the global population. Among them are the recent and ongoing infections of 15-110 million individuals due to increased prevalence in developed countries [2]. HEV comprises around 7.2 kb quasi-enveloped, positive sense, single-stranded

HEV mRNA Expression Profile Extraction and Differential Expression Analysis (DEA)
We accessed the NCBI-GEO [22] (https://www.ncbi.nlm.nih.gov/geo/, accessed on 1 August 2022) utilizing "HEV" and "Hepatitis E Virus" as appropriate keywords for HEV-associated mRNA expression profile retrieval. The search results were further trimmed down in accordance with the following inclusion criteria: (1) the dataset(s) must be "expression profiling by array" type along with HEV infection being present in "Homo Sapiens" as host; (2) the dataset must include raw and processed microarray data; (3) the dataset must have control and infected sample types; (4) the submission date of the dataset must be within last ten years (i.e., 2012 to 2022); (5) the dataset must contain greater than 25 samples. We excluded studies devoid of case reports, review articles, abstracts, non-human samples, and cell-line-based experimental study designs. The series matrix expression file of the selected dataset was retrieved from GEO, followed by quality checks. ARSyNseq function (with unknown batch settings) was available within NOISeq to acquire batch-corrected expression values. Mapping of probe IDs to their corresponding HUGO Gene Nomenclature Committee (HGNC) symbols was handled via GEO2R. Averaging expression values for those gene(s) mapping to several probe IDs was conducted to avoid redundancy. p-values and log 2 (fold change) values of all genes across two sample groups were computed utilizing limma [23]. The genes corresponding to a p-value < 0.05 and |log 2 (fold change)| > 0.1 were considered as differentially expressed across two sample groups. HGNC multi-symbol checker (https://www.genenames.org/tools/multi-symbolchecker/, accessed on 5 August 2022) was queried thereafter to obtain DEGs with officially approved HGNC symbols.

HEV-Specific WGCN Construction and Hub Module Selection
WGCN formation and module(s) selection are discussed in detail in [24]. Module eigengene (ME) and MEdissimilarity (MEdiss) were calculated. ME dendrogram was examined based on Pearson correlation to merge module(s) with comparable high expression profiles. Standard intramodular connectivity (k.in) and module membership (MM) were computed for every module. The module(s) with a significantly higher correlation between k.in and MM were regarded as the hub module(s).

PPI Network Construction and Modular Analysis
The hub module DEGs of both Mo and WB groups were collectively given as input to Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/, accessed on 21 August 2022) v11.5 web-based tool [25] in order to establish a PPI network corresponding to highest confidence (i.e., interaction score > 0.9) and subsequently visualized utilizing Cytoscape v3.91 [26]. The parameters used for PPI hub module selection are discussed in detail in [27].

Pathway and Gene Ontology (GO) Term Enrichment Analyses for Hub Gene Selection
All the PPI hub module genes were given as input to the Enrichr web server [28] (https: //maayanlab.cloud/Enrichr/, accessed on 25 August 2022), and the top 10 significantly (p-value < 0.05) enriched pathways and GO terms were collected from BioPlanet, GO-Biological Process (BP), GO-Molecular Function (MF), and GO-Cellular Compartment (CC) libraries. The genes overlapping between pathway, GO-BP, GO-MF, and GO-CC genesets were regarded as our hub genes.

HEV-Specific Three-Node miRNA FFL Construction and Analysis
ChEA v3.0 database [29] was queried in order to fetch significant [corresponding to score (p-value) < 0.05] human TFs regulating our hub mRNAs. miRNAs repressing our hub mRNAs and TFs (compiled from ChEA) were retrieved from miRWalk v3.0 [30] (miRNAs with score > 0.95 and binding only on 3 UTR region only), StarBase v2.0 [31], and miRTargetLink v2.0 [32] (strongly validated miRNAs) databases to form miRNA-TF and miRNA-mRNA pairs. A union of miRNAs from all these databases was taken in order to encompass a wide spectrum of miRNAs. These miRNAs were validated via literature, and only HEV-associated miRNAs were retained for further analysis. To establish a closed three-node miRNA FFL, all three interaction pairs (i.e., TF-mRNA, miRNA-mRNA, and miRNA-TF) were altered in such a manner that the FFL is composed of only overlapping miRNAs, TFs, and mRNAs. The miRNA FFL was thereafter visualized using Cytoscape.

Identification of Secondary Structures and Protein-Protein Docking Studies
The α-fold structure of human 2 -5 oligoadenylate synthase (OAS3) protein was extracted from Research Collaborators for Structural Bioinformatics Protein Data Bank (RCSB PDB) database (ID: AF_AFQ9Y6K5F1; residues 1-1087) with global pLDDT (model confidence score) of 86.82 ( Figure S1) [33,34]. Only the third domain of OAS3 (referred to as OAS3_D3), having the three aspartic acid residues Asp816, Asp818, and Asp888, which corresponds to the catalytic triad in human OAS1, are conserved ( Figure S5) [35]. Therefore, the third domain of OAS3 (residues 743-1087) was selected for protein-protein docking studies. The NMR-determined structure of the TRAF-type zinc finger domain containing 1TRAFD1 or FLN29 genes was downloaded from the RCSB PDB database (PDB ID: 2D9K; residues 1-75) ( Figure S5).  [36] (https://wenmr.science.uu.nl/prodigy/, accessed on 5 September 2022). The final docked complex was analyzed utilizing Protein Interactions Calculator (PIC) to confirm the interactive residues within the hydrogen-bond distance (3.5 ) [37]. The structure and docked complex with their interacting residues were visualized using PyMol [38]. The inhibition constant (K i ) for the docked complex were estimated by utilizing the formula: where, R (universal gas constant) is 1.985 × 10 −3 kcal mol −1 K −1 , ∆G is docking energy, and T (temperature) is 298.15 K [39].

HEV mRNA Expression Profile Extraction and DEA
Based on the abovementioned inclusion and exclusion criteria, we selected the HEV mRNA expression profile with accession number GSE36539. The dataset comprised 16 WB patient samples (i.e., eight WB control and eight WB infected) and 16 Mo patient samples (eight Mo control and eight Mo infected). For post-batch correction, we bifurcated the two cell types' data, with each set comprising 16 samples. After averaging expression values corresponding to duplicate genes, a total of 15991 and 15990 unique gene symbols were retained for Mo and WB sample groups. Corresponding to a p-value < 0.05 and |log 2 (fold change)| > 0.1, we obtained 190 and 2749 DEGs corresponding to Mo and WB sample groups. After checking via the HGNC multi-symbol checker and filtering by Pigengene, we were finally left with 93 and 2016 DEGs corresponding to Mo and WB sample groups. Figure 1

HEV-Specific WGCN Construction and Hub Module Selection
Mo-associated (93) and WB-associated (2016) HEV DEGs were passed as input to WGCNA. The WGCN was generated at β = 15 (corresponding to R 2 = 0.89) and β = 9 (corresponding to R 2 = 0.8) for Mo and WB sample groups, respectively. The clustering tree (hierarchical) and dynamic tree cut algorithm resulted in sixteen (i.e., black, blue, brown, cyan, green, green-yellow, grey, magenta, pink, midnight-blue, purple, red, salmon, tan, turquoise, yellow) and two (blue, turquoise) color-coded modules corresponding to WB and Mo sample groups. The merging of modules with highly co-expressed gene patterns was possible only in the WB group (due to a large number of genes) by cutting the dendrogram at the height of 0.25. After merging, the initial sixteen modules of the WB group were amalgamated into nine color-coded modules (i.e., black, blue, cyan, green-yellow, grey, purple, red, salmon, and turquoise). The Mo group's two modules were retained due to no possible merging between modules. Since the grey module consists of unassigned genes, we discarded it for further analysis. Figures S2 and S3 show plots for β in consideration with scale-free topology (SFT) criteria for Mo and WB sample groups. The clustering tree (hierarchical) of WB-DEGs based on dissTOM and ME with original (16) and merged (9) color modules is shown in Figure 2A. Figure 3A shows the clustering tree (hierarchical) of Mo-DEGs based on dissTOM with the original two modules. Figures 2B and 3B show a multidimensional scaling (MDS) plot of all modules in three scaling dimensions in WB and Mo. WGCN depicting TOM among WB and Mo modules is shown in Figures 2C and 3C. Based on the most significant correlation between

HEV-Specific WGCN Construction and Hub Module Selection
Mo-associated (93) and WB-associated (2016) HEV DEGs were passed as input to WGCNA. The WGCN was generated at β = 15 (corresponding to R 2 = 0.89) and β = 9 (corresponding to R 2 = 0.8) for Mo and WB sample groups, respectively. The clustering tree (hierarchical) and dynamic tree cut algorithm resulted in sixteen (i.e., black, blue, brown, cyan, green, green-yellow, grey, magenta, pink, midnight-blue, purple, red, salmon, tan, turquoise, yellow) and two (blue, turquoise) color-coded modules corresponding to WB and Mo sample groups. The merging of modules with highly co-expressed gene patterns was possible only in the WB group (due to a large number of genes) by cutting the dendrogram at the height of 0.25. After merging, the initial sixteen modules of the WB group were amalgamated into nine color-coded modules (i.e., black, blue, cyan, greenyellow, grey, purple, red, salmon, and turquoise). The Mo group's two modules were retained due to no possible merging between modules. Since the grey module consists of unassigned genes, we discarded it for further analysis. Figures S2 and S3 show plots for β in consideration with scale-free topology (SFT) criteria for Mo and WB sample groups. The clustering tree (hierarchical) of WB-DEGs based on dissTOM and ME with original (16) and merged (9) color modules is shown in Figure 2A. Figure 3A shows the clustering tree (hierarchical) of Mo-DEGs based on dissTOM with the original two modules. Figures 2B and 3B show a multidimensional scaling (MDS) plot of all modules in three scaling dimensions in WB and Mo. WGCN depicting TOM among WB and Mo modules is shown in Figures 2C and 3C. Based on the most significant correlation between MM and k.in (Tables S1 and S2), the turquoise module was chosen as our hub module in the case of both Mo and WB. Figures 2D and 3D show the heatmap plot of turquoise module genes along with their corresponding ME levels in the case of WB and Mo. A total of 51 and 457 DEGs were present in the hub turquoise module within Mo and WB.
MM and k.in (Tables S1-S2), the turquoise module was chosen as our hub module in the case of both Mo and WB. Figures 2D and 3D show the heatmap plot of turquoise module genes along with their corresponding ME levels in the case of WB and Mo. A total of 51 and 457 DEGs were present in the hub turquoise module within Mo and WB.   Progressively darker and lighter red colors within the matrix indicate higher and lower topological overlap among genes. Dark-colored blocks along the diagonal signify the modules. (D) Expression heatmap of turquoise module genes where the rows and columns correspond to genes and samples. The red and green color bands within the heatmap imply higher and lower expression levels across turquoise module genes. In addition, the corresponding ME expression levels (along the y-axis) across all samples (along the x-axis) are displayed at the bottom panel of the module heatmap in the form of a barplot.

PPI Network Construction and Modular Analysis
A total of 506 hub genes within the turquoise module were inputted into the STRING database, of which 477 were mapped to their corresponding proteins. Figure  4A shows the PPI network comprising 114 nodes and 545 edges (corresponding to an interaction score > 0.9). Molecular Complex Detection (MCODE) revealed a total of four modules, out of which the first module had the highest score and was considered our PPI hub module. Figure 4B shows the PPI hub module comprising 28 nodes and 374 edges. Essential centrality measures such as node degree, betweenness, closeness, clustering coefficient, topological coefficient, and average shortest path length of the PPI network are shown in Figure S4.

PPI Network Construction and Modular Analysis
A total of 506 hub genes within the turquoise module were inputted into the STRING database, of which 477 were mapped to their corresponding proteins. Figure 4A shows the PPI network comprising 114 nodes and 545 edges (corresponding to an interaction score > 0.9). Molecular Complex Detection (MCODE) revealed a total of four modules, out of which the first module had the highest score and was considered our PPI hub module. Figure 4B shows the PPI hub module comprising 28 nodes and 374 edges. Essential centrality measures such as node degree, betweenness, closeness, clustering coefficient, topological coefficient, and average shortest path length of the PPI network are shown in Figure S4.

Pathway and GO Term Enrichment Analyses for Hub Gene Selection
A total of 25, 28, 20, and 21 genes within our PPI hub module participated in the top 10 significant pathways, GO-BP, GO-MF, and GO-CC terms, as shown in Tables S3-S6. The most significant pathway, GO-BP, GO-MF, and GO-CC terms were Interferon α/β signaling ( − value = 1.89 × 10 −58 ), cellular response to type I interferon ( − value = 2.24 × 10 −73 ), adenylyltransferase activity ( − value = 2.21 × 10 −08 ), mitochondrial envelope ( − value = 2.82 × 10 −05 ). A total of 13 hub genes overlapped between these top 10 significant pathways, GO-BP, GO-MF, and GO-CC genesets, shown by a Venn plot in Figure S4. The box-and-whisker plots showing the relative expression distribution of these HEV-hub genes are shown in Figure S5.

HEV-Specific Three-Node miRNA FFL Analysis
Our HEV-specific -node miRNA FFL, shown in Figure 5A, comprised 16 nodes and 51 edges. TF-mRNA, miRNA-mRNA, and miRNA-TF pairs embodied 31, 14, and 6 edges among the whole FFL. Amongst all the FFL nodes, 8, 4, and 4 belonged to mRNAs, miRNAs, and TFs. The degree range of miRNAs, mRNAs, and TFs in the FFL varied from 3 to 7, 5 to 7, and 8 to 10, respectively. The average degrees of TFs, mRNAs, and miR-NAs were 9.25, 5.625, and 4. The three pairs of regulatory relationships between miR-NAs, mRNAs, and TFs were summarized in Table S7. Out of all TFs, both TRAFD1 and ETV7 were equally regulated by the highest number of miRNAs (i.e., 2). Among all the hub DEGs, OAS3 was regulated by the highest number of miRNAs (i.e., 3). Tables S8-S10 show the top three miRNAs, mRNAs, and TFs within FFL ranked based on degree, betweenness, closeness, eigenvector, and radiality. Based on the observations from the table, the highest-order subnetwork motif comprised one TF (TRAFD1), two miRNAs (miR-222-3p and miR-125b-5p), and one mRNA (OAS3), as shown in Figure 5B.

Pathway and GO Term Enrichment Analyses for Hub Gene Selection
A total of 25, 28, 20, and 21 genes within our PPI hub module participated in the top 10 significant pathways, GO-BP, GO-MF, and GO-CC terms, as shown in Tables S3-S6. The most significant pathway, GO-BP, GO-MF, and GO-CC terms were Interferon α/β signaling (p-value = 1.89 × 10 −58 ), cellular response to type I interferon (p-value = 2.24 × 10 −73 ), adenylyltransferase activity (p-value = 2.21 × 10 −8 ), mitochondrial envelope (p-value = 2.82 × 10 −5 ). A total of 13 hub genes overlapped between these top 10 significant pathways, GO-BP, GO-MF, and GO-CC genesets, shown by a Venn plot in Figure S4. The box-andwhisker plots showing the relative expression distribution of these HEV-hub genes are shown in Figure S5.

HEV-Specific Three-Node miRNA FFL Analysis
Our HEV-specific-node miRNA FFL, shown in Figure 5A, comprised 16 nodes and 51 edges. TF-mRNA, miRNA-mRNA, and miRNA-TF pairs embodied 31, 14, and 6 edges among the whole FFL. Amongst all the FFL nodes, 8, 4, and 4 belonged to mRNAs, miRNAs, and TFs. The degree range of miRNAs, mRNAs, and TFs in the FFL varied from 3 to 7, 5 to 7, and 8 to 10, respectively. The average degrees of TFs, mRNAs, and miRNAs were 9.25, 5.625, and 4. The three pairs of regulatory relationships between miRNAs, mRNAs, and TFs were summarized in Table S7. Out of all TFs, both TRAFD1 and ETV7 were equally regulated by the highest number of miRNAs (i.e., 2). Among all the hub DEGs, OAS3 was regulated by the highest number of miRNAs (i.e., 3). Tables S8-S10 show the top three miRNAs, mRNAs, and TFs within FFL ranked based on degree, betweenness, closeness, eigenvector, and radiality. Based on the observations from the table, the highest-order subnetwork motif comprised one TF (TRAFD1), two miRNAs (miR-222-3p and miR-125b-5p), and one mRNA (OAS3), as shown in Figure 5B.

Protein-Protein Docking Analysis and Interaction Studies
The predicted binding residues obtained by CASTp v3.0 corresponding to TRAFD1 protein having solvent accessible surface area 10.735 Å 2 and for OAS3_D3 protein, the catalytic triad-Asp816, Asp818, and Asp888 were further analyzed for docking (Tables S11). The nine best clusters generated by docking for the OAS3_D3-TRAFD1 complex are shown in Table S12. The clusters with a low HADDOCK score, low Z-score, and low RMSD values were considered the best-docked complex. The model with a HADDOCK score of −91.7 +/− 6.1 Kcal/mol, RMSD of 0.7 +/− 0.5 Å, and Z-score of −2.4 was selected for OAS3_D3-TRAFD1 complex ( Figure 6A). The predicted binding affinity (ΔG) and dissociation constant ( K d ) of the OAS3_D3-TRAFD1 complex calculated by the PRODIGY server resulted in values −11.0 kcal/mol and 9.1E − 09 M at 25 ℃, respectively. Moreover, the inhibition constant of the OAS3_D3-TRAFD1 complex was computed to be 8.472E − 09. Analysis of docked complex (OAS3_D3-TRAFD1) using the PIC server showed the presence of interacting residues showing hydrogen bonding interactions between the side chain of OAS3_D3 and the main chain of TRAFD1 (Table S13a) and between side chains of OAS3_D3 and TRAFD1 (Table S13b). OAS3_D3 is shown to interact with TRAFD1 with buried surface area of 1706.5 +/− 71.8 Å 2 ( Figure 6B). Fourteen amino acid residues of OAS3_D3 are shown to form H-bonds with backbone atoms of TRAFD1 ( Figure 6C). Furthermore, seven amino acid residues OAS3_D3 form H-bonds with side chains of TRAFD1 ( Figure 6D). Interestingly, among these hydrogen bond interactions, the side chains of catalytic triad residues Asp814, Asp818, and Asp888 of OAS3_D3 jointly interact with side chains of Arg67 of TRAFD1, leading to the formation of salt-bridge-mediated interactions.

Protein-Protein Docking Analysis and Interaction Studies
The predicted binding residues obtained by CASTp v3.0 corresponding to TRAFD1 protein having solvent accessible surface area 10.735 2 and for OAS3_D3 protein, the catalytic triad-Asp816, Asp818, and Asp888 were further analyzed for docking (Table  S11). The nine best clusters generated by docking for the OAS3_D3-TRAFD1 complex are shown in Table S12. The clusters with a low HADDOCK score, low Z-score, and low RMSD values were considered the best-docked complex. The model with a HADDOCK score of −91.7 +/− 6.1 Kcal/mol, RMSD of 0.7 + /− 0.5 , and Z-score of −2.4 was selected for OAS3_D3-TRAFD1 complex ( Figure 6A). The predicted binding affinity (∆G) and dissociation constant (K d ) of the OAS3_D3-TRAFD1 complex calculated by the PRODIGY server resulted in values −11.0 kcal/mol and 9.1E-09 M at 25°C, respectively. Moreover, the inhibition constant of the OAS3_D3-TRAFD1 complex was computed to be 8.472E-09. Analysis of docked complex (OAS3_D3-TRAFD1) using the PIC server showed the presence of interacting residues showing hydrogen bonding interactions between the side chain of OAS3_D3 and the main chain of TRAFD1 (Table S13a) and between side chains of OAS3_D3 and TRAFD1 (Table S13b). OAS3_D3 is shown to interact with TRAFD1 with buried surface area of 1706.5 +/− 71.8 Å 2 ( Figure 6B). Fourteen amino acid residues of OAS3_D3 are shown to form H-bonds with backbone atoms of TRAFD1 ( Figure 6C). Furthermore, seven amino acid residues OAS3_D3 form H-bonds with side chains of TRAFD1 ( Figure 6D). Interestingly, among these hydrogen bond interactions, the side chains of catalytic triad residues Asp814, Asp818, and Asp888 of OAS3_D3 jointly interact with side chains of Arg67 of TRAFD1, leading to the formation of salt-bridge-mediated interactions. Genes 2023, 14, x FOR PEER REVIEW 10 of 15 Figure 6. Docked mode of OAS3_D3-TRAFD1 complex and its structural mapping of intramolecular interfaces. (A) Cartoon representation of OAS3_D3-TRAFD1 docked complex having low energy score determined by HADDOCK docking. OAS3_D3 is displayed in blue color, and TRAFD1 is represented in magenta color. (B) A surface representation of OAS3_D3-TRAFD1 complex. (C) Interacting residues between OAS3_D3 and TRAFD1 obtained through hydrogen bonding between the side chain of OAS3_D3 and the main chain of TRAFD1 are shown in licorice representation, OAS3_D3 residues (blue) and TRAFD1 residues (in magenta). (D) Interacting residues between OAS3_D3 and TRAFD1 obtained through hydrogen bonding between the side chain of OAS3_D3 and the side chain of TRAFD1 are shown in licorice representation, OAS3_D3 residues (blue) and TRAFD1 residues (in magenta). The atomic distances are shown in yellow dotted lines in angstrom (Å ).

Discussion
We primarily accessed GEO to extract a dataset containing mRNA from WB and Mo samples of kidney transplant patients, each from control and infected. Our analysis identified 93 and 2016 DEGs corresponding to Mo and WB sample groups, which were then rigorously screened using WGCNA to select hub modules containing 51 and 457 DEGs from Mo and WB sample groups, respectively. Furthermore, the PPI network modular analysis indicated four modules, out of which the first module was considered our PPI hub module. Pathway and GO term enrichment analyses of the PPI hub module revealed thirteen overlapping genes between the top 10 significant pathways, GO-BP, MF, and CC genesets. These thirteen upregulated genes-Guanylate BindingProtein 2 (GBP2),Interferon Regulatory Factor 1 (IRF1), Interferon Regulatory Factor 5 (IRF5), Interferon Regulatory Factor 7 (IRF7), Interferon Regulatory Factor 9 (IRF9), Interferon Stimulating Gene 20 (ISG20), MX Dynamin Like GTPase 1 (MX1), MX dynamin Like GTPase 2 (MX2), 2'-5'-Oligoadenylate Synthase 1 (OAS1), 2'-5'-Oligoadenylate Synthase 2 (OAS2), OAS3, SAM Domain And HD Domain-Containing Protein 1 (SAMHD1), and Signal Transducer And Activator Of Transcription 1 (STAT1) regarded as the hub genes were obtained. HEVassociated miRNAs and TFs related to these hub genes were also detected. Within our work, we utilized key transcriptional level interaction, i.e., miRNA→gene, miRNA→TF, and TF → gene regulation, for constructing a miRNA-centered FFL, which portrays a critical role in HEV pathogenesis. Interacting residues between OAS3_D3 and TRAFD1 obtained through hydrogen bonding between the side chain of OAS3_D3 and the main chain of TRAFD1 are shown in licorice representation, OAS3_D3 residues (blue) and TRAFD1 residues (in magenta). (D) Interacting residues between OAS3_D3 and TRAFD1 obtained through hydrogen bonding between the side chain of OAS3_D3 and the side chain of TRAFD1 are shown in licorice representation, OAS3_D3 residues (blue) and TRAFD1 residues (in magenta). The atomic distances are shown in yellow dotted lines in angstrom (Å).

Discussion
We primarily accessed GEO to extract a dataset containing mRNA from WB and Mo samples of kidney transplant patients, each from control and infected. Our analysis identified 93 and 2016 DEGs corresponding to Mo and WB sample groups, which were then rigorously screened using WGCNA to select hub modules containing 51 and 457 DEGs from Mo and WB sample groups, respectively. Furthermore, the PPI network modular analysis indicated four modules, out of which the first module was considered our PPI hub module. Pathway and GO term enrichment analyses of the PPI hub module revealed thirteen overlapping genes between the top 10 significant pathways, GO-BP, MF, and CC genesets. These thirteen upregulated genes-Guanylate BindingProtein 2 (GBP2), Interferon Regulatory Factor 1 (IRF1), Interferon Regulatory Factor 5 (IRF5), Interferon Regulatory Factor 7 (IRF7), Interferon Regulatory Factor 9 (IRF9), Interferon Stimulating Gene 20 (ISG20), MX Dynamin Like GTPase 1 (MX1), MX dynamin Like GTPase 2 (MX2), 2 -5 -Oligoadenylate Synthase 1 (OAS1), 2 -5 -Oligoadenylate Synthase 2 (OAS2), OAS3, SAM Domain And HD Domain-Containing Protein 1 (SAMHD1), and Signal Transducer And Activator Of Transcription 1 (STAT1) regarded as the hub genes were obtained. HEVassociated miRNAs and TFs related to these hub genes were also detected. Within our work, we utilized key transcriptional level interaction, i.e., miRNA → gene, miRNA → TF, and TF → gene regulation, for constructing a miRNA-centered FFL, which portrays a critical role in HEV pathogenesis.
A limited number of studies have been conducted concerning miRNA regulation during chronic HEV infection. Based on the centrality measures, we identified miR-222-3p (member of the miR-221/222 family) and miR-125b-5p as the top-ranked miRNAs engaged in the FFL regulation involving TRAFD1 and OAS3 in HEV. The miR-222-3p has already been implicated in regulating downstream acute HEV in vitro [40]. A study by Elghoroury et al. predicted the diagnostic value of miR-222 as a potential biomarker in diagnosing liver injuries or progression, cirrhosis, and Hepatocellular Carcinoma (HCC) caused by viral infections [41,42]. A case report conducted by Lin et al. discovered that re-infected chronic HEV patient promotes rapid development of HCC and cirrhosis without Hepatitis B Virus (HBV) [43]. This led to the fact that HEV re-infection might boost the progression of HCC in patients regardless of chronic HBV infection. Thus, showing the significance of miR-222 as a potential biomarker in the case of HEV. In addition, a recent study on the acute and chronic HEV-infected serum showed that the miR-125b-5p plays a pivotal role in the suppression of cell proliferation, replication, and apoptosis in HEV and can act as a biomarker in early detection and differentiation between acute and chronic infection [44] which in turn provide us with another potential miRNA as a candidate for biomarker.
Most of the hub genes detected in our work are associated with interferon (IFNs) signaling pathways. Among the genes induced by IFNs, oligoadenylate synthase (OAS) proteins have been identified as enzymes that mediate antiviral effector functions. In humans, the OAS family proteins OAS1, OAS2, and OAS3 are involved in RNase L activation pathways for the degradation of cellular and viral RNA. Recent reports suggested that the 110-kDa OAS3 protein gets activated at lower RNA concentrations than 33-kDa OAS1 and, on average, synthesizes longer 2 -5 -linked oligoadenylates (2-5As) to activate RNase L [35]. This is further supported by the previous studies that established the potential roles of OAS3 p100 in mediating the antiviral activity of the dengue virus [45] and HCV [46] in an RNase L-dependent manner. Interestingly, both Dengue virus and HCV are from the same family, Flaviviridae, and could induce hepatitis [47].
In addition to miR-125b-5p and miR-222-3p interaction with the OAS3, our network analysis also observed TRAFD1 (Tumor necrosis factor receptor-associated factor-type zinc finger domain containing 1) as a critical TF implicated in the regulatory aspects of HEV infection. TRAFD1 expression is inducible by interferon and suppresses Toll-like receptor 4-mediated NF-κB activation by binding to tumor necrosis receptor associated factor 6 (TRAF6) [48]. A study undertaken by Green et al. performed OAS1b-dependent immune transcriptional profiles of West Nile Virus (WNV) infection, which resulted in the hypothesis that TRAFD1 may contribute to innate immune protection mediated by the OAS1b network [49]. In lung adenocarcinoma (LUAD), IFN-γ response genes, TRAFD1 is reported to exhibit significant expression and a strong positive correlation with OAS3 [50]. As it has demonstrated relevance in non-viral diseases, it would be interesting to untangle its effects on viral infections, especially HEV. Since HEV infection in organ-transplanted patients can trigger the onset of immune responses, TRAFD1 can thus act as a master regulator to control the excessive innate immune response [51].
Our first hypothesis revolves around the association of miR-125b-3p, miR-222-3p, and OAS3. The miRNAs finalized based on the obtained results have a significant role in HEV, as described by the previous studies. It has been shown that miR-125b-3p upregulates in the case of chronic HEV infection, whereas miR-222 downregulates [40] by enhancing the OAS3 expression in HEV. In addition, the upregulation of OAS3 was observed and associated with miR-125b-5p, miR-222-3p, and TRAFD1 in our study. As miR-125b-5p and OAS3 are upregulated, we can hypothesize that miR-125b-5p positively regulates the OAS3 (having antiviral properties). This hypothesis was supported by an earlier study demonstrating that miRNAs can influence viral pathogenesis by either directly modifying viral gene expression or promoting cellular antiviral responses [52]. There is a strong possibility of miR-222-3p showing positive regulation with OAS3 as well, which is supported by the study showing miR-222 downregulation by enhancing the target gene in HEV [40], which needs further validation through in vitro and in vivo analysis.
Another hypothesis revolves around the connection of TRAFD1 with OAS3. TRAFD1 shows a negative feedback regulation related to excessive immune response [48]. It thus may be involved in the possible inhibition of upregulated OAS3 as OAS3 is an interferon stimulating gene (ISG) associated with innate immunity. To verify the above hypothesis, we conducted protein-protein docking where the inhibition constant (K i ) with 8.472E-09 value and dissociation constant (K d ) with −11.0 kcal/mol values showed a strong inverse association between TRAFD1 and OAS3 protein.
Our third hypothesis connects the miR-125b-3p and miR-222-3p with TRAFD1. Here as the miRNAs are showing a positive regulation with OAS3, there is a possibility that miRNAs are involved in inhibiting the TRAFD1 as TRAFD1 negatively regulates OAS3. The two miRNAs can work in combination or individually to affect OAS3 and TRAFD1.
Based on the above hypotheses, we can conclude that miR-125b-3p and miR-222-3p positively regulate the OAS3 and negatively regulates TRAFD1 resulting in higher expression of OAS3 against HEV. Although to prove this hypothesis, further validation is required through in vitro studies.

Conclusions
This study consigns the initial report of miRNA signatures determined in WB and Mo samples of kidney-transplanted chronic hepatitis E virus-infected patients. To our knowledge, it is the first time that miR-222-3p has been described as a reference miRNA for use in CHEV sample analyses. The three-node miR-222-3p/miR-125b-5p-OAS3-TRAFD1 regulatory network provides a novel insight into understanding the molecular mechanism of chronic HEV. Further experimental studies are needed to confirm the importance of their role in CHEV infection.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes14010042/s1, Figure S1: Three-dimensional model of human OAS3 and TRAFD1 proteins and their respective active site binding pocket residues; Figure  S2: Quality check plots to assess value of β with respect to scale-free topology (SFT) criteria for Monocytes (Mo) DEGs; Figure S3: Quality check plots to assess value of β with respect to scalefree topology (SFT) criteria for Whole Blood (WB) DEGs; Table S1: Module membership (MM) vs intramodular connectivity (k.in) correlation and p-values for Monocytes (Mo) modules; Table S2: Module membership (MM) vs intramodular connectivity (k.in) correlation and p-values for Whole Blood (WB) modules; Figure S4: Centrality measures showing node degree distribution, betweenness, closeness, clustering coefficient, average shortest path length, neighborhood connectivity, topological coefficient, and radiality of PPI network; Table S3: Top 10 significant pathways (based on p-values) and their gene counts associated with PPI hub module genes; Table S4: Top 10 significant GO-BP terms (based on p-values) and their gene counts associated with PPI hub module genes; Table S5: Top 10 significant GO-MF terms (based on p-values) and their gene counts associated with PPI hub module genes; Table S6: Top 10 significant GO-CC terms (based on p-values) and their gene counts associated with PPI hub module genes; Figure S5: Common 13 HEV-hub genes between top 10 significant GO term and pathway genesets and their relative expression distribution shown by box-and-whisker plots across control and infected patient Whole Blood (WB) samples; Table S7: Summary of regulatory relationships between HEV-associated mRNAs, miRNAs, and TFs; Table  S8: Top 3 miRNAs ranked based on centrality measures such as degree, betweenness, closeness, radiality, and eigenvector; Table S9: Top 3 mRNAs ranked based on centrality measures such as degree, betweenness, closeness, radiality, and eigenvector; Table S10: Top 3 TFs ranked based on centrality measures such as degree, betweenness, closeness, radiality, and eigenvector; Table S11: Active residues of domain 3 of OAS3 (AF_Q9Y6K5-F1) and FLN29 or TRAFD1 (PDB ID: 2D9K) proteins identified using literature and predicted using CASTp web server respectively; Table S12: Statistical analysis for HADDOCK generated OAS3-TRAFD1 docked complexes; Table S13: Interfacial residues in docked complex of OAS3-TRAFD1. Protein-protein main chain/side chain-side chain hydrogen bonds. Chain A represent OAS3_D3 protein and chain B represent TRAFD1 protein. Data Availability Statement: The dataset used in our work is available in NCBI-GEO at https: //www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36539, accessed on 1 August 2022, and can be accessed with GSE36539.