Differential expression analysis and profiling of hepatic miRNA and isomiRNA in dengue hemorrhagic fever

Dengue virus causes dengue hemorrhagic fever (DHF) and has been associated to fatal cases worldwide. The liver is one of the most important target tissues in severe cases, due to its intense viral replication and metabolic role. microRNAs role during infection is crucial to understand the regulatory mechanisms of DENV infection and can help in diagnostic and anti-viral therapies development. We sequenced the miRNome of six fatal cases and compared to five controls, to characterize the human microRNAs expression profile in the liver tissue during DHF. Eight microRNAs were differentially expressed, including miR-126-5p, a regulatory molecule of endothelial cells, miR-122-5p, a liver specific homeostasis regulator, and miR-146a-5p, an interferon-regulator. Enrichment analysis with predicted target genes of microRNAs revealed regulatory pathways of apoptosis, involving MAPK, RAS, CDK and FAS. Immune response pathways were related to NF- kB, CC and CX families, IL and TLR. This is the first description of the human microRNA and isomicroRNA profile in liver tissues from DHF cases. The results demonstrated the association of miR-126-5p, miR-122-5p and miR-146a-5p with DHF liver pathogenesis, involving endothelial repair and vascular permeability regulation, control of homeostasis and expression of inflammatory cytokines.

Scientific Reports | (2021) 11:5554 | https://doi.org/10.1038/s41598-020-72892-w www.nature.com/scientificreports/ in liver leads to the development of an environment with multiple cells expressing IL-6 and TGF-beta, which may be related to apoptosis of hepatocytes 7 . Both IL-6 and TGF-beta are cytokines with key roles in inflammation and are probably associated to increased vascular permeability, plasma leaking and consequent worsening of dengue symptoms 8 . In dengue severe cases an increase of pro inflammatory and vasoactive cytokines as IFN-γ, TNF-α, IL-1, IL-6, IL-8 IL-10, CCL2 (MCP-1) and CCL5 is usually observed. Although not conclusive, hypothesis that hepatic inflammation during dengue infections is related to hepatocytes damage protection and tissue repairing to the homeostasis reestablishment is sustained 9 . During the immune system action in response to infections, there is an increase in the pathway communication among cells, as well as an intense cytoplasmic microenvironment modification in order to regulate gene expression in response to the agent. The microRNAs (miRNA) are RNA molecules that act in this scenario of RNA interference (RNAi), being key mediators of host response to infectious diseases. The miRNA sequences are single-stranded RNA, 22 nucleotides in length (average), associate to the RNA-induced silencing complex (RISC) and regulate gene expression by binding to the region between the base 2 and 8 (seed region). Its target is the miRNA recognition element (MRE) located at the 3 'UTR of mRNA 10 .
Alternative or imprecise biological processing of miRNA and post-transcriptional modifications may generate sequence variations, the isomiRNAs. The same miRNA locus can give rise to multiple distinct isomiRs, differing in length or sequence composition that might present different functionalities 11 .
Thus, the aim of this study was to characterize the differential miRNA expression profile of human liver tissue from fatal cases of DHF, identify the isomiR profile, and propose the regulatory mechanisms of these molecules in metabolic pathways related to the pathogenesis of the disease.
Most of the differentially expressed miRNAs were under-expressed in DHF libraries in comparison to the control condition. miR-122-5p showed the highest expression in all controls (log2FC = −6.59; FDR < 0.001), with average CPM = 225,797.2 in controls and 2,336.4 in DHF, that is, it is almost 100 times less expressed in DHF group. miR-10b-5p is the second with the highest expression level in the control group, decreased almost 100 times in DHF (log2FC = −6.167 (FDR < 0.001), followed by miR-146a-5p (log2FC = −2.6476 (FDR = 0.0387), both had high CPM values. The expression of miR-204-5p (log2FC = -5.632; FDR < 0.001), miR-148a-5p  Table S1). The expression of hepatic miRNA could correctly cluster the sample groups in DHF and controls, as showed by the MDS plot, with exception of the sample Dhf_299 which was excluded from the analysis (Supplementary Figure S1). Heatmap analysis based on log2(CPM) of also resulted in distinct clades of DHF and control groups (Fig. 1).
Target genes prediction and functional analysis. The target prediction step generated a list of 224 target genes of overexpressed miRNAs in liver were the input dataset to functional annotation. The enrichment analysis identified 26 pathways with significant p-value < 0.05 after Benjamini adjustment: nine pathways of apoptosis and cell cycle regulation; seven related to the immune response including pathways activated by NF-κβ, the other ten pathways are mainly related to cell signaling and signal transduction ( Fig. 2A, Supplementary Table S2). Functional analysis involving target genes of underexpressed miRNAs in DHF group was carried out with 570 target genes, from six downregulated miRNAs, obtained in the three databases and resulted in 208 pathways, with the predominance of cell death pathways (Fig. 2B, Supplementary Table S3).
isomiRNA differential expression analysis. At total, we detected 5500 differentially expressed isomiRs, however only 188 passed our Log2FC and p value criteria. Of these, 66 were canonical sequences and 122 were isomiRs. From the 8 miRNAs previously associated with DHF in hepatic profile: four presented isomiRs (hsa-miR-122-5p, hsa-miR-126-5p, hsa-miR-146a-5p and hsa-miR-423-5p), two (hsa-miR-204-5p and hsa-miR-148a-5p) presented only their canonical sequence and two (hsa-miR-133a-5p and hsa-miR-10b-5p) were not considered differentially expressed in our samples. The isomiRs fractions of the each four miRNAs associated with DHF, described here, is demonstrated in Fig. 3. Table 2. Top 30 microRNAs most expressed in hepatic tissue of dengue hemorrhagic fever (DHF) and controls, expression is presented by counts per million (CPM). *Statistically significant.  Table 3 as well as their respective categories and relevant statistical data. The counts were normalized in CPM and represent the group mean and the top was arranged from lowest to highest adjusted p values. Concerning the software isomiRs classification, sequences containing nt addition (uppercase) or deletion (lower case) in the UTR ends are classified as 5′ (t5) or 3′ isoforms (t3). Addition anywhere else is ad, substitution is mm and if there are a discrepancy in the 5′ end 2-8nt is seed represented by the nt position and the modification. The canonical sequence is ref.
For example, the mature sequence of miR-122-5p is if the sequence of the anotated miRNA is UGG AGU GUG ACA AUG GUG UUUG. The annotation miR-122-5p.iso.t3:T corresponds to an addition (T) in the 3′end, resulting in the sequence UGG AGU GUG ACA AUG GUG UUUG T.    miRNA 122-5p obtained 28 isoforms differentially expressed, and the miRNA 99a-5p eight found isoforms. The expression pattern and the divergences for the miRNA 122-5p are exposed in Fig. 4.
The majority of the miRNAs that were found differentially expressed in this study were downregulated in DHF. This was expected, since DENV can inhibit the expression levels of Dicer, Drosha, Ago1 and Ago2, key players in the miRNA biogenesis, resulting in the suppression several miRNAs 12 . Kakumani et al. detected the downregulation of 143 miRNAs and the overexpression of only 3 miRNAs. They could validate the suppression of some miRNAs by qRT-PCR during DENV infection, including miR-10b-5p, also found to be downregulated in the present study.
Kanokudom et al. also studied the regulation of selected miRNAs during DENV infection in HepG2 cells, and observed the downregulation of miR-146a 13 . Surprisingly, the miR-146a-5p is usually up-regulated in dengue fever, targeting the inflammatory signaling pathway to inhibit type I IFN-mediated antiviral defense by regulating NF-kb. It's up-regulation might serve as a mechanism of DENV to establish its initial infection and later dissemination 14 . Also, is well established that autophagy, a cellular degradation and recycling process 15 , is required for efficient DENV replication 16 . Pu et al. found a negative correlation between the levels of miR-146a and DENV-2 induced autophagy by targeting TRAF6, suggesting that this is necessary in order to suppress excessive inflammation and reduce pathological damage 16 . The downregulation of miR-146a in DHF, detected in this study, could be increasing DENV-2 autophagy and consequently promoting the inflammation response and pathological damage.
Only one study assessed the expression of miRNAs in patients with DHF 17 . In that study, the authors used mononuclear leukocytes derived from patients with and without DHF to evaluate eleven miRNAs that targets SOCS1, which in turn is a negative regulator of several cytokines involved in the immune response 18 . SOCS1 acts by inhibiting JAK, an activator of a family of transcription factors known as STATs that induces the transcription of several cytokines 18 . In the present study, miR-204-5p, previously described as a JAK2 inhibitor 19 , was found to be downregulated and miR-126, which activates the JAK1/STAT3 pathway by inhibiting ZEB1 20 was found www.nature.com/scientificreports/ to be upregulated. Taken together, these findings suggest that the tight regulation of cytokines by SOCS-1 and JAKs/STATs cascade might have a direct impact in the development of DHF. The miR-126-5p is located in the intronic Egfl7 gene. This region encodes the guide sequence miRNA-126-3p (named miR-126) and the passenger 5p. Both Egfl7 gene and miR-126-3p/5p are angiogenic factors involved in maintaining vascular integrity 21 . The miR-126-3p had one differentially expressed isomiR and the miR-126-5p had five, all up regulated in DHF. The growth and endothelial repair regulation by miR-126-5p (and miR-126-3p) are well established in literature, as well as the vascular permeability increase in dengue hemorrhagic pathogenesis, in this way is evident this possible relationship between miR-126-5p super-expression and DHF vascular symptoms 22 . It is possible that vascular permeability in DHF be regulated by miR-126-5p, however, due the lack of data related to gene expression and cytokine profile from samples in this study, further experimental analysis to elucidate the real mechanism that describe this association are necessary.
The miRNA miR-122-5p was the most expressed miRNA in the samples showing 28 differentially expressed isomiRs, 26 down-regulated in DHF group (as the canonical sequence) and two up-regulated. The miR-122-5p is a liver-specific miRNA that represents more than 50% of liver miRnome and plays important role in fat metabolism and homeostasis. Rats with this gene silenced developed fatty liver hepatitis, fibrosis and hepatocellular carcinoma 23 , that agrees with histopatologic results of steatosis, necrosis, hyperplasia and metabolic changes from DHF fatal cases 8 . The up-regulated isomiRs, were found in high levels in the DHF liver samples and in very low level in the control samples, this indicates that those isomiRs can play a different role in this condition, likely related to inflammation balancing acting towards reducing the production of vasoactive cytokines.
The enrichment analysis result of KEGG pathways from down regulated miRNA target genes in dengue (miR-122-5p, miR-10b-5p, miR-204-5p, miR-148a-5p, miR-146a-5p and miR-423-5p) described 29 pathways associated with cell death regulation, from which ten were found among the 30 best scoring pathways. This indicate a possible role of these miRNAs on DHF, once apoptosis is one of the main characteristics of dengue hemorrhagic pathogenesis 4 .
One sub-expressed, miR-148a-5p, was associated with cancer, and can also be related to cell proliferation control pathways 25 , but were never described in association to dengue pathogenesis as miR-126-5p, miR122-5p and miR146a-5p, an thus its role remains to be determined. Nevertheless, it showed statistically significant differences even with stringent analysis criteria and deserve further studies to clarify this relationship.
Of note, some miRNA present divergent expression pattern when compared to their respective isomiRs, suggesting that those divergent isomiRs might act differently from their canonical miRNA in these conditions. This analysis demonstrates that those expression results can be obscured by a DE analysis focusing only on total miRNA or canonical miRNA.

Conclusion
Our study provided insight into the hepatic miRNA and isomiR profile in DHF, made a deep differential expression analysis, characterized the isomiRs representation for previously miRNA associated with the disease, and among the eight miRNAs differentially expressed in DHF human liver. We highlight miRNAs miR-126-5p, miR-122-5p and miR-146a-5p, which seems to be involved in liver pathogenesis of DHF, because the revised function in growth and endothelial repair linked to dengue due the vascular permeability presented in bleeding cases. The central role of miR-122-5p in liver homeostasis and interferon regulation of miR-146a-5p associated to dengue by its immune response involvement.
MiRNAs and isomiRNAs are promising agents for anti-viral therapy and diagnostics development, experimental studies are necessary to identify target genes regulated by these miRNAs in liver during DENV infection, also isomiRs discovery may alter the current knowledge of immune modulation that leads to severe dengue and fatal cases.

Material and methods
This research study has been approved by the ethics committee of the Medical School of the University of São Paulo, Brazil (local number 253/12 CEP/FMUSP 08/2012). All methods were performed in accordance with the relevant guidelines and regulations of the University of São Paulo and Instituto Evandro Chagas.

Samples.
All procedures were carried out with the informed written consent of the next of kin. Liver tissues samples (n = 6) from patients who died from DHF during a DENV-2 outbreak in São Paulo State in 2010 were provided by the Department of Pathology from the Guilherme Álvaro Hospital (São Paulo, Brazil). Tissue samples were obtained by performing a viscerotomy within 12 h following death of patients and preserved using formalin-fixation and paraffin-embedding (FFPE). DENV infection was confirmed by clinical, serological (IgM and NS1 antigen detection by ELISA) and histological examination (hematoxylin-eosin staining detection of viral antigens by immunohistochemistry assay). Control FFPE liver tissue samples (n = 5) were obtained from patients that died of acute myocardial infarction. Histological analysis (performed within 12 h following death) showed no infection or hepatic injury (data not available). Data processing. Quality control of the raw reads was performed using FastQC (https ://www.bioin forma tics.babra ham.ac.uk/proje cts/fastq c) and FASTX-toolkit (https ://hanno nlab.cshl.edu/fastx _toolk it/index .html) was used to trim 3′ adapters, remove unresolved bases (N) and eliminate reads shorter than 17nt. Collapsing of redundant sequences and mapping against the human genome (version hg19-UCSC) was performed by the Mapper module available in miRDeep2 software 27 . miRDeep2 was also used to calculate read counts and to annotate the mature miRNA. By the end of the processing step, only miRNAs with 18 to 26nt length were selected. The SeqBuster-mirAligner v.1.2.2 28 allowed the identification of isomiRs by performing a survey on miRNA sequence variability in the samples, applying default parameters.

Differential expression analysis.
To determine the miRNAs expression levels, the miRNAs counts obtained by miRDeep2 were normalized using Counts per million (CPM). The differential expression analysis was performed with the edgeR package for R Bioconductor 29 . Exact Fisher and GLM (General Linear Model) statistical tests were applied with FDR (False Discovery Rate Benjamini-Hochberg) to adjust the p-value. True differential expressed miRNA was defined by a log2-transformed fold change (log2FC) ≥ 2 (or ≤ -2) and an FDR value < 0.05.
Prediction of target genes and functional analysis. TargetScan (www.targe tscan .org), miRDB (www. miRDB .org) and miRTarBase (www.miRTa rBase .mbc.nctu.edu.tw) databases were used to predict the targets of the differentially expressed miRNA. Only targets predicted by at least two of the mentioned databases considered valid and proceeded to functional annotation and enrichment analysis with DAVID 6.7 (https ://david .ncifc rf.gov/home.jsp)30. Only functional categories containing pathways with p-value (with Bejamini correction) < 0.05 were evaluated (Additional File 4). The interaction network of target genes and pathways was built with Cytoscape 3.3 platform (https ://www.cytos cape.org/).