SARS-CoV-2 causes a significant stress response mediated by small RNAs in the blood of COVID-19 patients.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has had a serious impact on the world. In this study, small RNAs from the blood of COVID-19 patients with moderate or severe symptoms were extracted for high-throughput sequencing and analysis. Interestingly, the levels of a special group of tRNA-derived small RNAs (tsRNAs) were found to be dramatically upregulated after SARS-CoV-2 infection, particularly in coronavirus disease 2019 (COVID-19) patients with severe symptoms. In particular, the 3'CCA tsRNAs from tRNA-Gly were highly consistent with the inflammation indicator C-reactive protein (CRP). In addition, we found that the majority of significantly changed microRNAs (miRNAs) were associated with endoplasmic reticulum (ER)/unfolded protein response (UPR) sensors, which may lead to the induction of proinflammatory cytokine and immune responses. This study found that SARS-CoV-2 infection caused significant changes in the levels of stress-associated small RNAs in patient blood and their potential functions. Our research revealed that the cells of COVID-19 patients undergo tremendous stress and respond, which can be reflected or regulated by small non-coding RNA (sncRNAs), thus providing potential thought for therapeutic intervention in COVID-19 by modulating small RNA levels or activities.

Introduction

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) outbreak has had a tremendous impact on health worldwide. SARS-CoV-2 infection causes coronavirus disease 2019 (COVID-19), which can be simply divided into four categories: mild, moderate, severe, and critical. To date, many researchers have been devoted to exploring the pathogenic mechanism of SARS-CoV-2, and it has been found that nearly 28% of deaths caused by COVID-19 are due to the cytokine storm and sepsis induced by SARS-CoV-2 infection. It has been reported that the molecular pathogenesis and pathophysiology of COVID-19 include direct viral invasion, renin-angiotensin-aldosterone system (RAAS) dysregulation, hypoxia, hyperinflammation, cytokine storm, endotheliopathy, and thrombosis. A subset of COVID-19 patients manifests a dysfunctional “hyperinflammatory” response with persistent fevers and elevated inflammatory marker levels. Severe COVID-19 patients show high levels of proinflammatory cytokines (such as interleukin-1 beta [IL-1β], IL-2, IL-6, IL-18, interferon gamma [IFN-γ], and tumor necrosis factor alpha [TNF-α]) and chemokines (CXCL8, CXCL10, macrophage inflammatory protein 1α [MIP1α], MIP1β, and monocyte chemoattractant protein-1 [MCP1])., However, a separate subset of severe patients manifests a dysfunctional “hyperinflammatory” response with an absence of fevers and, worse, multiple organ failure. A study showed significant changes in serum C-reactive protein (CRP) levels in critically ill patients with COVID-19. Despite considerable efforts to elucidate coronavirus mechanisms, our understanding of COVID-19 pathogenesis and disease processes is still greatly lacking.

A large number of studies have shown that microRNAs (miRNAs), a type of small non-coding RNA (sncRNA) with a length of 18–22 nt, play an important role in the process of virus infection. They induce the degradation or translational inhibition of target transcripts by binding to their 3′ UTR (untranslated region) to destroy receptor genes for virus entry, regulate virus replication or translation cofactors, affect initial or adaptive immune genes, silence genes triggering cell apoptosis, or promote/inhibit viral replication by targeting the viral genome sequence. To date, there have been many studies on the abnormal expression of miRNAs in infectious diseases and the pathogenic mechanism, including enterovirus 71 (EV71) infections, sepsis, chronic hepatitis C,11, 12, 13 and various infectious diseases, implying that they have potential clinical applications in infectious diseases in the future. Recently, a “nose cocktail” of miRNA was suggested to treat COVID-19 patients.

In addition to miRNAs, tRNA-derived small RNAs (tsRNAs) are a novel kind of stable functional small molecule produced by reprocessing the precursors or mature transcripts of tRNA whose levels have been found to be induced by stress conditions (including temperature, starvation, virus infection, and cancer)., For instance, tsRNA has been reported to be generated in cells under the stress of virus infection and utilized by pathogens to improve their infection efficiency.19, 20, 21 Recent studies have found that tsRNAs can target primer-binding sites (PBSs) of endogenous retroviruses (ERVs) or human T cell leukemia virus type 1 (HTLV-1), thereby protecting the genome, or be utilized by viruses to promote their self-synthesis.

As the two most abundant small RNAs in somatic cells, miRNAs and tsRNAs are closely related to virus infection and the immune response. However, to date, the changes and functions of these small RNAs after SARS-CoV-2 infection are still unclear. To address this question, we performed high-throughput sequencing to identify the repertoire of small RNAs expressed in the blood of moderate/severe COVID-19 patients compared with that of healthy controls and predicted the potential post-transcriptional regulatory network of the significantly changed small RNAs.

Results

The vast majority of small RNAs in COVID-19 patients are 21 nt and 22 nt

We first collected blood samples from 10 COVID-19 patients and 4 healthy people. COVID-19 patients were divided into moderate and severe case groups according to the clinical classification, which is shown in Table S1. Next, we extracted small RNAs from each sample's cell lysates and performed high-throughput sequencing (Figure 1A). Raw data were subjected to adaptor trimming and quality control analysis (Table 1). The majority (nearly 99%) of reads originated from humans (Table 1), while small RNAs barely matched the SARS-CoV-2 genome.

Figure 1

Small RNA analysis and length distribution

(A) Flowchart of the small RNA extraction, preparation, sequencing, and analysis process. (B) Histogram of the length distribution of small RNAs in each group of samples, in which the types of small RNAs with lengths of 21 nt and 22 nt are displayed separately.

Table 1

Reads number and alignment rate in each sample

Sample ID	Clinical type	Total number of reads	Mapped to human genome	Overall alignment rate
H1	Healthy	26,130,020	25,932,012	99.24
H2	Healthy	33,220,438	33,081,293	99.58
H3	Healthy	25,507,591	25,417,197	99.65
H4	Healthy	27,511,432	27,383,911	99.54
M1	Moderate	42,539,283	42,341,314	99.53
M2	Moderate	32,905,763	32,716,049	99.42
M3	Moderate	23,063,364	22,936,188	99.45
M4	Moderate	42,938,961	42,383,435	98.71
M5	Moderate	28,662,930	28,534,787	99.55
M6	Moderate	24,586,678	24,432,128	99.37
S1	Severe	19,633,465	19,413,071	98.88
S2	Severe	22,074,028	21,882,645	99.13
S3	Severe	22,220,577	21,940,039	98.74
S4	Severe	27,059,405	26,923,097	99.5

We next analyzed the characteristics of these small RNAs and found that small RNAs with 21 nt and 22 nt lengths accounted for the largest proportion in each group (Figure 1B). Further results showed that 97.05% of all small RNAs were miRNAs. Interestingly, a small batch (accounting for 0.55%) of 22 nt tsRNAs gradually increased in abundance with the progression of the disease. We next investigated the changes in these small RNAs in different groups of samples.

The levels of stress-induced tsRNAs are significantly upregulated after SARS-CoV-2 infection

All the identified small RNAs were compared among the different groups. To our surprise, the most significantly changed small RNAs were tsRNAs, and their maximum change in expression exceeded 200-fold (Figures 2A and 2B). In addition, the expression of these tsRNAs was significantly upregulated in the patients after SARS-CoV-2 infection, which was even more dramatic in the severe group (Figures 2A and 2B). For example, the expression level of tRNA-Ala tsRNAs in healthy people was very low (mean RPM [reads per million mapped reads] = 33.82). However, in the COVID-19 patients, its average expression level (mean RPM value) rose to 321.61 (Figure 2C). It was more surprising that tsRNAs derived from tRNA-Val and tRNA-Asp were almost not expressed in either healthy controls or moderate cases (mean RPM = 0), but were significantly highly expressed in severe cases (mean RPM = 692.23 and 1735.43, p value = 1.26 × 10−3 and 1.38 × 10−3) (Figures 2D and 2E).

Figure 2

The expression and change of tsRNAs in each group of samples

(A and B) Volcano plot of (A) COVID-19 moderate cases versus healthy controls and (B) COVID-19 severe cases versus healthy controls. (C–E) Expression value of significantly changed tsRNAs between different sample groups. The red dots are the mean expression value of tsRNAs in the samples of the corresponding group. (F–H) The type and length distribution of tsRNAs in each group of samples. The x axis represents the length of tsRNAs (unit: nt). The y axis represents the percentage of tsRNAs.

tsRNAs can be divided into five types according to their different sources, tRF-5 (5′ tRF), tRF-3 (3′ CCA tRF), tRF-1 (3′ trailer tRF), tRF-i (internal tRF), and tRH (tRNA half), each of which has a related function.24, 25, 26 To investigate whether the upregulated tsRNAs in COVID-19 patients are type specific, the source and length distribution of tsRNAs were further analyzed. To our surprise, the majority of 22 nt tsRNAs were derived from tRF-3 (3′CCA tRF), and their expression increased by approximately 20% in severe cases (Figures 2F–2H).

In particular, the expression level of tsRNAs derived from the 3′CCA region of tRNA-Gly (Figure 3A and 3B) gradually increased with the development of the disease (Figure 3C). Compared with those in healthy controls, their levels were upregulated more than 4-fold in severe patients according to the sequencing data (Figure 3C) and more than 200-fold according to qRT-PCR results (Figure 3D). More interestingly, we analyzed 30 clinical characteristics of patients (Table S1) and finally found that in patients with a high level of the inflammation marker C-reactive protein (CRP), this tsRNA was significantly upregulated (logFC = 1.28, p value = 7.9 × 10−3) (Figure 3E). The levels of tsRNA are often induced under stress conditions; thus we next wanted to know, in addition to these stress-induced tsRNAs, what role does the large number of miRNAs play in disease?

Figure 3

The characteristics of tRNA-Gly-derived tsRNAs

(A) The distribution and abundance of sequencing reads mapping to mature tRNA-Gly. The x axis represents the position on mature tRNA-Gly-GCC. The y axis represents the RPM value of the reads measured at this position. (B) The position of tRNA-Gly-tRF-3 in the secondary structure of tRNA-Gly. (C) Boxplot of the expression levels of tRNA-Gly-tRF-3 among samples determined by next-generation sequencing (NGS). (D) Bar plot of the expression levels of tRNA-Gly-tRF-3 among samples determined by qRT-PCR. (E) Boxplot of the expression levels of tRNA-Gly-tRF-3 in patients with increased CRP levels and normal patients.

Highly expressed miRNAs are significantly downregulated in severe COVID-19 patients

Principal component analysis results showed that, based on the expression level of miRNAs, 10 COVID-19 patients and 4 healthy people could be clearly divided into different groups (Figure 4A). For COVID-19 patients, moderate and severe cases could be basically grouped into two categories, indicating that the expression level of miRNAs in the blood was significantly different between COVID-19 patients and healthy controls, which was consistent with the clinical classification of patients. To illustrate that the expression level of these miRNAs is mainly driven by the disease rather than sex, treatment methods, or other factors, we conducted different regression modeling (Figure S1). Across the miRNAs tested, we observed no association between the miRNA expressions and sex and some slight γ-globulin treatment effects on the miRNA expressions (Figure S1A). Nevertheless, in a joint regression including these factors as covariates in the patients, the γ-globulin treatment effects on the miRNA expressions all can be captured by the disease status. Namely, the variation in the miRNA expressions is mainly driven by the disease status.

Figure 4

The expression profile of miRNAs in healthy controls, moderate patients, and severe patients

(A) Principal component analysis of miRNAs in different groups of samples. (B) (Left) Heatmap of the expression level of miRNAs in each sample. The unit of expression value is log2RPM. The black box in the figure indicates the group of highly expressed miRNAs. Clustering was performed based on the expression value of miRNAs and split into three gene sets. (Middle) Boxplot of the expression of each miRNA in all samples. The green background indicates the low-expression miRNAs (mean value less than the mid-range), and the yellow bottom indicates the high-expression miRNAs (mean value greater than the mid-range). (Right) The names of the miRNAs and the HI-DOWN miRNAs are highlighted.

Furthermore, these miRNAs were divided into high- and low-expression groups according to their average expression levels compared with the middle of the expression range of these miRNAs (Figure 4B). Interestingly, six highly expressed miRNAs all were downregulated in severe patients. We therefore named these miRNAs HI-DOWN miRNAs (highly expressed and significantly downregulated in severe patients).

Target molecules of differentially expressed miRNAs are enriched in the endoplasmic reticulum stress pathway

Based on miRNA function, we investigated 6,053 experimentally validated target genes of the differentially expressed miRNAs and built interaction networks among them. Furthermore, Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed that these genes were significantly enriched in the endoplasmic reticulum (ER) lumen (Figure 5A) and involved in the "protein digestion and absorption" and "protein processing in ER" pathways (Figure 5B).

Figure 5

Functional enrichment analysis of miRNA-target genes

(A) Cellular component analysis results of differentially expressed miRNA-target genes. (B) Results of KEGG pathway analysis of miRNA targets enriched in the ER lumen. (C) The interaction network between miRNAs and their target genes involved in the “protein digestion and absorption” and “protein processing in the ER” pathways. The circle in the figure denotes the miRNA, the triangle denotes the target gene, and the different colors represent the function of the target gene and the pathway they are involved in. Different edge types distinguish whether the MTI (miRNA-target interaction) is functional. Nodes of different sizes indicate the degree of interaction. (D) Boxplot of the expression levels of HI-DOWN miRNAs in each group of samples. The y axis indicates the log2RPM expression value.

Studies have found that coronavirus proteins can cause damage to the ER and affect the secretory pathway., In addition, the ER has been suggested as a potential therapeutic target for COVID-19 management., We therefore speculated whether the activation of ER-related pathways may be regulated by miRNAs. We thus established the regulatory network of target genes involved in key signaling pathways in the ER lumen. Surprisingly, the hub miRNAs (degree >4) in the network all were consistent with the HI-DOWN miRNAs (Figure 5C), whose expression was found to be downregulated in severe patients (Figure 5D). For instance, miR-335-5p may regulate 24 ER-related genes, including sensor genes in the three downstream pathways of the unfolded protein response (UPR), including the PERK pathway, ATF6 pathway, and IRE1 pathway (Figures 5 and S2). Furthermore, the enrichment analysis of the target genes of all differentially expressed miRNAs showed that they were involved in a series of downstream signaling pathways induced by the stress response, including inflammation and apoptosis pathways, such as the phosphatidylinositol 3-kinase (PI3K)-AkT pathway, p38 mitogen-activated protein kinase (MAPK) pathway, Janus kinase (JAK)-signal transducer and activator of transcription (STAT) signaling pathway, and transforming growth factor β (TGF-β) signaling pathway (Figure S3).

Discussion

SARS-CoV-2 has posed a tremendous threat to human health. This study revealed the key molecules of host cells in the response against viruses from the perspective of small RNAs. We found two types of small RNAs with significant changes, tsRNAs and miRNAs, which are both related to the stress response and downstream pathway changes, including inflammation, immunity, and even apoptosis. Upon virus infection, the internal and external environment of cells undergoes drastic changes, which might cause the cell to deviate from homeostasis or even induce programmed cell death. Our research indicates that small RNAs are involved in these severe stress and response processes. Although we are using whole-blood samples, which contains multiple cell types, we believe sncRNA-induced stress response is a common pathway in cells and the affection of these stress-related small RNAs should be shared by most types of cells. Future research can sort cells on this basis and then study the changes in these small RNAs separately to verify our conclusion. In the future, in-depth verification of the expression and function of these small RNAs will be needed to clarify the role and mechanisms of small RNAs in the occurrence and development of COVID-19.

The most recent studies have shown that miRNAs can be used as indicator molecules for detecting infections by a variety of pathogens, including Hendra virus infection, HIV infection, tuberculosis, malaria, and Ebola. Our study revealed that the expression profile of miRNAs can accurately distinguish different groups of COVID-19 patients. This finding was consistent with another study showing that miRNAs can predict the prognosis of mild/severe cases and can also identify asymptomatic infections in HBsAg carriers and chronic hepatitis B patients.

As the most important type of sncRNAs in cells, miRNAs are estimated to regulate 60% of all protein-coding genes at the post-transcriptional level and therefore play an important role in cell development, differentiation, and immune processes. In this study, we found that some HI-DOWN miRNAs in the blood of COVID-19 patients have been previously reported to have antiviral functions. For example, miR-181 can target the receptor gene CD163 to affect virus entry, miR-130 can regulate hepatitis C virus (HCV) cofactors and affect its replication, miR-24 can directly interact with the viral genome to inhibit its translation, and miR-29 is an antiviral factor induced by IL-21 in CD4 T cells. In addition, the expression level of serum miR-29 is negatively correlated with liver fibrotic stages and necroinflammation grades in patients with chronic hepatitis B virus (HBV) infection. Studies have confirmed that miR-335 directly regulates the stress response.,

The UPR is induced by ER stress caused by the accumulation of misfolded protein, which is highly activated after virus infection. Based on the results that miRNAs are enriched in the ER stress pathway, it can be speculated that after SARS-CoV-2 infection, cells may downregulate the expression of a series of miRNAs to activate the downstream pathway of the UPR. Previous studies have revealed that the p38 MAPK pathway is involved in viral replication and that the downregulation of miR-29 expression is involved in the activation of caspases (CASP) and affects cell apoptosis. In addition, the TGF-β signaling pathway, affected by miRNAs, can also induce apoptosis. The activation of these pathways further amplifies the stress response of cells. As all differentially expressed genes are involved in a series of downstream signaling pathways of stress, we believe that miRNA-target genes should control the specificity, timing, and severity of the cell response to SARS-CoV-2 infection stresses. Considering the extensive and cascade amplification effect of miRNAs on target regulation, the use of miRNAs for the treatment and monitoring of infectious diseases can be considered, which may also lead to a different way to treat COVID-19” would be appropriate.

Our research revealed that small RNAs from the blood of COVID-19 patients are mainly produced by the human genome, not by the SARS-CoV-2 genome. Growing evidence has shown that a variety of viruses can produce small viral RNAs (svRNAs), including influenza, EV71, hepatitis A virus (HAV), West Nile virus, and SARS-CoV. These svRNAs can use the host's RNAi mechanism to control the mRNAs in cells, which in turn affects the fate of both viruses and host cells. In addition, some studies have speculated on potential svRNAs derived from the SARS-CoV-2 genome by analysis of the virus genome sequence.56, 57, 58 In our transcriptome data, we could barely detect small RNAs from SARS-CoV-2, which may be caused by two hypothetical possibilities. First, SARS-CoV-2 cannot produce small RNAs after entering host cells. Second, SARS-CoV-2 can produce svRNAs only in the lungs or other invading tissues but cannot release them into the blood, or its production and release are affected by other factors. In the future, it will be necessary to detect the presence of SARS-CoV-2 small RNA in cell lines or lung tissues to verify the above hypothesis.

We conducted a “power analysis” to investigate whether the existing number of duplicate samples in our study can support our conclusion. The results of this type of analysis have been widely recognized., We first downloaded the data of Zheng et al. on miRNA analysis of COVID-19 patients during recovery, which contains peripheral blood mononuclear cell (PBMC) data from 6 mild, 7 moderate, and 5 severe patient blood samples. We performed the power analysis of the three group sets separately (moderate versus mild, severe versus mild, and severe versus moderate) using the ssize.fdr R package. The results showed that along with the increase in sample size, the power value also gradually increased (Figure S4A). When the sample size was greater than 5, the power value basically reached 1. When the sample size was 4, the power value of each group ranged from 0.81 to 0.92. In general, when the power value is greater than 0.8, the conclusion is considered plausible. It follows that for the analysis of blood small RNA data for COVID-19, to obtain reliable conclusions, the minimum number of samples in each group should be 4. The number of samples in each group in our study was 4 healthy, 6 moderate, and 4 severe. Theoretically, the power value obtained by each group should be greater than 0.8. To prove this, we next performed a power analysis between groups of our research data (Figure S4B). The results showed that the power values of our results were 0.95 (moderate versus healthy), 0.88 (severe versus healthy), and 0.97 (severe versus moderate), which are all greater than 0.8. Therefore, we believe that the sample size used in our study supports our conclusions. However, additional samples and experiments are needed before conducting clinical research. We will continue to follow up the progress in this area and hope that through the collaboration of researchers worldwide, we can finally reveal the mechanism of SARS-CoV-2 triggering the body's response from various perspectives.

Materials and methods

COVID-19 patient clinical information collection

The clinical data of patients with COVID-19 treated in the Fifth Affiliated Hospital of Sun Yat-sen University in Zhuhai from January 17, 2020, to March 2, 2020, were collected. In this retrospective study, the clinical, biochemical, and immunological characteristics of 10 patients with COVID-19 were analyzed. These patients were classified as severe (4 cases) or moderate (6 cases). We diagnosed COVID-19 in accordance with the criteria in the World Health Organization (WHO) interim guidelines and the Diagnosis and Treatment Plan for Novel Coronavirus Pneumonia by the National Health Commission (seventh trial version).

Total RNA isolation and small RNA sequencing

RNA was isolated after red blood cell lysis from peripheral blood by using TRIzol (Life Technologies, cat. 265,709, CA, USA). Total RNA was isolated from cell lysates and purified using TRIzol following the manufacturer's procedure. After quality inspection with an Agilent 2100 Bioanalyzer (Agilent, cat. G2939AA, CA, USA) and Qubit assay tubes (Life Technologies, cat. 1604220, CA, USA), 10 ng of total miRNA was used as the input for library preparation with a NEXTflex Small RNA-Seq Kit v3 (Bio Scientific Corporation, cat. NOVA-5132-05, TX, USA) following the manufacturer's instructions. We added a linker to the 3′ end and then added another linker to the 5′ end after removing and inactivating the excess 3′ linker. miRNAs with adaptors were reverse transcribed to generate cDNA by reverse transcription reagent. After purification by NEXTflex Cleanup Beads (Bio Scientific Corporation, cat. NOVA-5132-05, TX, USA), the ligated products were then amplified with PCR by the following conditions: initial denaturation at 95°C for 2 min; 22–25 cycles of denaturation at 95°C for 20 s, annealing at 60°C for 30 s, and extension at 72°C for 15 s; and final extension at 72°C for 2 min. The average insert size for the final cDNA library was 155 ± 5 bp. After purification by NEXTflex Cleanup Beads (Bio Scientific Corporation, cat. NOVA-5132-05, TX, USA), quality control of concentration and fragment size was performed with an Agilent 2100 Bioanalyzer (Agilent, cat. G2939AA, CA, USA) and Qubit assay tubes (Life Technologies, cat. 1604220, CA, USA). Finally, we performed 2 × 150 bp paired-end sequencing (PE150) on an Illumina NovaSeq 6000 (Guangzhou Huayin Health Technology, Guangzhou, China) following the vendor's recommended protocol.

Quality control and adaptor trimming

The raw FASTQ files were subjected to adaptor trimming and low-quality and short read (<18 bp) removal using Fastp 0.20.1. The human reference genome GRCH38 and SARS-CoV-2 reference genome NC_045512.2 were indexed with Bowtie2's default indexing parameters. For alignment, all clean reads were mapped to the human and SARS-CoV-2 reference genomes using Bowtie2 with default parameters. Bam files were transformed and sorted by Samtools 1.9.

miRNA expression quantification

We used miRDeep2 v.2.0.1.3 to process and quantify the clean reads. The human genome (assembly GRCh38) and human miRBase v.22 precursor and mature miRNA sequences were used as reference files. The clean reads were processed with the mapper.pl script. The reads shorter than 18 nt were discarded. Quantification and expression profiling of processed reads were performed by the quantifier.pl script. We allowed three mismatches when mapping reads to precursors. Counts of reads mapping to each miRBase mature miRNA were obtained from the miRDeep2 output and normalized by its default method.

tsRNA identification and quantification

tsRNA expression profiles were constructed using the tRFinder tool, which we established previously., Briefly, the clean data were mapped to the human genome (assembly GRCh38), and the unmatched data were mapped to tRNA precursors. According to the length and location of these tRFs, we classified them as tRF-5, tRF-3, tRF-1, or tRF-i. The read counts were normalized into RPM values. The RPM is evaluated using Equation 1:where R represents the number of reads mapped to a particular small RNA, and N represents the total number of all mapped reads from given library. All the log-transformed RPM values were performed using log2-transformed RPM plus 1.

Differential gene expression analysis

Differentially expressed RNAs were determined using the Limma package in R. The value of the expression matrix was RPM, and the matrix was transformed to log2-expression values. We selected 32 differentially expressed miRNAs (p value < 0.05). These miRNAs were split into three classes (class 1 [moderate high pattern], class 2 [healthy high pattern], and class 3 [severe high pattern]) using the ComplexHeatmap package in R.

miRNA target and enrichment analysis

For the host, a human miRNA-target interaction (MTI) file was obtained from miRTarBase v.8.0. All the MTIs were intersected with differentially expressed miRNAs, and the intersections were retrieved. Enrichment analysis of miRNA-target genes was performed using ClusterProfiler v.3.14.3 in R v.3.6.3. GO enrichment analysis was performed using the function enrichGO, and KEGG enrichment analysis was performed using the function enrichkegg. The interaction network between miRNAs and their target genes was constructed using Cytoscape software v.3.7.1.

Power analysis

Power analysis was performed with the ssize.fdr R package. The power calculation was performed with the ssize.twoSamp function. It was assumed that the actual difference between the average expression in the different group pairs as well as the standard deviations of the expression were identical for all miRNAs. The common value of the standard deviations was estimated from the data using the function sd_pooled (effectsize R package) and conservatively set to the 90th percentile of the residual standard deviations of gene expression. The proportion of non-differentially expressed genes (π0) was estimated using the function pi0est (qvalue R package) on the vector of p values obtained from the differential expression analysis in the different group pairs, and the mean of π0 was used for power calculation. Both datasets used 3 as the delta argument, and a false discovery rate (FDR) controlled at 10% was used. Other arguments of ssize.twoSamp were set to default.

Regression modeling

The function lm is used to fit linear models in R for regression. In univariable regressions of the miRNA expressions on γ-globulin treatment and multivariable regression, healthy individuals are excluded in order to clarify the effect of treatment on the disease. The visualization is done with function qqplot.

Ethical approval and consent to participate

The study was approved by the Research Ethics Committee at the Fifth Affiliated Hospital of Sun Yat-Sen University to allow retrospective access to patients' records and files (No. ZDWY2020 Lunzi No•K22-1).