Conserved MicroRNAs in Human Nasopharynx Tissue Samples from Swabs Are Differentially Expressed in Response to SARS-CoV-2

The use of high-throughput small RNA sequencing is well established as a technique to unveil the miRNAs in various tissues. The miRNA profiles are different between infected and non-infected tissues. We compare the SARS-CoV-2 positive and SARS-CoV-2 negative RNA samples extracted from human nasopharynx tissue samples to show different miRNA profiles. We explored differentially expressed miRNAs in response to SARS-CoV-2 in the RNA extracted from nasopharynx tissues of 10 SARS-CoV-2-positive and 10 SARS-CoV-2-negative patients. miRNAs were identified by small RNA sequencing, and the expression levels of selected miRNAs were validated by real-time RT-PCR. We identified 943 conserved miRNAs, likely generated through posttranscriptional modifications. The identified miRNAs were expressed in both RNA groups, NegS and PosS: miR-148a, miR-21, miR-34c, miR-34b, and miR-342. The most differentially expressed miRNA was miR-21, which is likely closely linked to the presence of SARS-CoV-2 in nasopharynx tissues. Our results contribute to further understanding the role of miRNAs in SARS-CoV-2 pathogenesis, which may be crucial for understanding disease symptom development in humans.


Introduction
Severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) is currently a global threat leading to considerable disease and mortality worldwide. SARS-CoV-2 is a close relative of SARS-CoV with 45-90% sequence similarity and has resulted in more than 8000 confirmed cases of severe acute respiratory syndrome in [1,2]. Coronaviruses are a diverse family of viruses associated with multiple respiratory diseases with different severities, such as the common cold, pneumonia, and COVID-19 [3].
SARS-CoV-2 belongs to the order Nidovirales, family Coronaviridae, subfamily Orthocoronavirinae, genus Betacoronavirus, and subgenus Sarbecovirus. It has a single-stranded positive-sense RNA genome, 26-32 kilobases (kb) in length [4]. Similar to an mRNA, the virus genome consists of a 5 cap structure together with a 3 poly(A) tail that translates its proteins. Most of the genome at the 5 end is occupied by the replicase gene, which encodes polyproteins ORF1a and ORF1b. ORF1a and ORF1b are further processed to generate nonstructural proteins (nsps). ORF1a contributes to the production of nsp1-nsp11, while the rest of the nsps (nsp12-nsp16) originate from ORF1b [5]. Additionally, viral structural proteins comprise surface (S), envelope (E), membrane (M), and nucleocapsid (N) proteins encoded by one-third of the genome at the 3 end [5,6].
Small noncoding endogenous RNAs called microRNAs (miRNAs) play a major role in posttranscriptional gene regulation related to diverse biological processes, including

RNA Samples, RNA Extraction, and Real-Time RT-PCR Detection of SARS-CoV-2
Total RNA was extracted from the nasopharynx using flocked swabs and a 1 mL viral transport medium (various vendors). RNA was extracted using an in-house method adapted from He et al. [17] using Sera-Mag beads (GE Healthcare, Chicago, IL, USA). Realtime RT-PCR was performed according to the diagnostic detection of Wuhan coronavirus 2019 by real-time RT-PCR [18]. RdRp, E, and human RNAseP assays were run in separate real-time RT-PCR reactions using a GB OneStep IPC Elite real-time RT-PCR Kit (Generi Biotech, Hradec Kralove, Czech Republic), which employs an internal positive control for the detection of PCR inhibition. For relative quantification of the viral titer in positive samples, the E gene assay results were normalized to hsa-mir-148a housekeeping miRNA. The reverse transcription and real-time PCR conditions for the hsa-mir-148a assay were, except for the fact that the RNA was not diluted, the same as described below in the validation of miRNA expression profiles by real-time RT-PCR section. The relative quantity of viral titer was calculated according to the 2 −∆∆C t method [19] by using the qBase+ software (Biogazelle, Zwijnaarde, Belgium). In total, the study was conducted on 10 SARS-CoV-2-positive (PosS) and 10 SARS-CoV-2-negative (NegS) RNA samples listed in Table 1.

Small RNA Library Preparation and High-Throughput Sequencing
A small RNA library was constructed using the NEBNext ® Small RNA Library Prep Set for Illumina ® (NEB, Ipswich, MA, USA), and purification was performed with the TailorCut Gel Extraction Tool Set (SeqMatic, Fremont, CA, USA). The quality and quantity of the library were determined using the Agilent High Sensitivity DNA Kit (Agilent, Santa Clara, CA, USA). The quantity of the library was also determined by a Modulus™ Single Tube Multimode Reader (Turner Biosystems, Sunnyvale, CA, USA) using a Quant-iT™ dsDNA Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA) and finally with an MCNext™ SYBR ® Fast qPCR Library Quantification Kit (MCLAB, San Francisco, CA, USA) used with Rotor-Gene 3000 (Corbett Research, Sydney, Australia). All kits were used according to the manufacturer's instructions. The libraries were pooled at a concentration of 2 nM according to fluorimetry measurements, assuming that the final cloned small RNA products were~150 bp.
For each sequencing run, a final pooled library of small RNAs consisted of pooled samples. In total, 20 RNA samples were sequenced, 10 SARS-CoV-2-negative and 10 positive RNA samples, which were previously tested by real-time PCR. In total, six runs were performed on a MiniSeq instrument (Illumina, San Diego, CA, USA). The MiniSeq High Output Reagent Kit, 75 cycles (Illumina, San Diego, CA, USA) providing 36-nt long reads, was used. Table 1. List of the RNA samples. Ct values were generated with the E gene assay [18]. The relative quantity of the virus was generated with the E gene results normalized to hsa-mir-148a miRNA.

Bioinformatics and Data Evaluation
Sequence quality was controlled by using FastQC-0.10.1 [20]. Then, the reads were imported to CLC Genomics Workbench 6.5.1 (CLC Bio, Aarhus, Denmark) using the following parameters: sequence length 10-50 nucleotides, no ambiguous nucleotides, removal of the Illumina universal adapter sequence, and a Phred score assigned a Q score of 30 (Q30); reads were trimmed using the following parameters: removal of smallRNA_adapter (TGGAATTC), removal of sequences of length: minimum 19 nucleotides and maximum 25 nucleotides. The total number of known miRNAs was counted and annotated using miRbase-Release 22.1 (Homo sapiens) in CLC Genomics Workbench 6.5.1 (CLC Bio, Aarhus, Denmark). The statistical method to quantify differential expression in CLC Genomics Workbench 6.5.1 was used: Transcriptomic Analysis: Small RNA Analysis: Extract and Count and Annotate and Merge Counts. The numbers of miRNA sequences were normalized to 1 million reads (RPM) in order to enable a comparative analysis. Statistics as Kruskal-Wallis test, ANOVA, PCA, and scattered plot graphs were performed using PAST version 2.17c [21].

Validation of miRNA Expression Profiles by Real-Time RT-PCR
For validation of miRNA expression profiles, hsa-mir-21, hsa-mir-34b, hsa-mir-34c, and hsa-mir-342 and a housekeeping hsa-mir-148a were chosen. These miRNAs were selected according to the statistical significances of small RNA sequencing results. Prior to reverse transcription, the RNA samples were treated with DNase I (Sigma-Aldrich, Saint-Louis, MO, USA). From each sample, 1 ng of RNA was used for reverse transcription of each miRNA by using the TaqMan MicroRNA Reverse Transcription Kit Assay (Applied Biosystems, Foster City, CA, USA) according to the manufacturer's instructions. Oligonucleotides for reverse transcription and real-time PCR of each miRNA were supplied in TaqMan MicroRNA Assays (Applied Biosystems, Foster City, CA, USA), where the assay IDs for tested miRNAs were as follows: ID: 000397 for hsa-mir-21; ID: 000427 for hsa-mir-34b; ID: 000428 for hsa-mir-34c; ID: 002147 for hsa-mir-342; ID: 000470 for hsa-mir-148a. Real-time PCR amplification conditions and reactions using TaqMan Universal Master Mix II (Applied Biosystems, Foster City, CA, USA) were performed according to the TaqMan MicroRNA Assay (Applied Biosystems, Foster City, CA, USA) instructions. Each miRNA sample was amplified in triplicate and run in a qTower cycler (Analytik Jena, Jena, Germany). Prior to miRNA expression measurement, the primer pair efficiency (E) values were evaluated on the standard curves of serial dilutions of pooled cDNA for each miRNA. For miRNA expression normalization, the ∆∆Ct method [22] was used with hsa-mir-148a miRNA as a normalizer using qPCRsoft 3.4 software (Analytik Jena, Jena, Germany). The normalized expression values were then statistically evaluated by analysis of variance (ANOVA) (p-value ≤ 0.05) in Statistica 13 software (Tibco, Palo Alto, CA, USA). The same software was used in the analysis of correlation (Pearson's correlation coefficient) between relative viral quantity and normalized miRNA expressions.

The Abundance of miRNAs Detected in Nasopharynx Tissues
The selected samples from the SARS-CoV-2 negative group (NegS) were assayed by using the SARBECO primer detection protocol [18] and were determined to be negative based on a Ct value threshold of 40. Samples from the SARS-CoV-2-positive group (PosS) were assayed by the same real-time RT-PCR method and the Ct values ranged from 24.82 (G9) to 33.75 (F5). The median Ct of the PosS group was 30.71. After normalization, the relative viral quantities ranged in the PosS group between 2.2 (F5) to 18529.1 (A7) showing almost a 10,000× difference in the viral quantity within the samples ( Table 1). The median relative viral quantity of the PosS group was 90.1.
In the present study, libraries representative of small RNA populations extracted from nasopharynx tissues and sequenced by Illumina SBS technology contained a quality of Q30 following the number of reads: in the NegS group ranging from 510,211 (3A) to 5,844,091 (13A) reads and in the PosS group ranging from 178,950 (D3) to 5,007,768 (H5) reads, the normalized numbers of reads are shown in Table 2. Sequencing data are deposited under BioProject acc. no. PRJNA747809, where the SRA experiments are available by acc. nos. SRX11490726-SRX11490735 (NegS group) and SRX11493797-SRX11493806 (PosS group).
The reads were normalized to reads per million values. The most abundant sRNAs belonged to the 21-and 22-nt classes in general (Figures 1 and 2). The highest median read value was for the 21-nt class sRNAs in both the NegS (9,598 reads) and PosS (6331 reads) groups. However, the average value in the NegS group (Figures 1 and 2) was the highest for the 22-nt class (13,261 reads), and that in the PosS group was the highest for the 21-nt class (15,088 reads) sRNAs (Table 2). The reads were normalized to reads per million values. The most abundant sRNAs belonged to the 21-and 22-nt classes in general (Figures 1 and 2). The highest median read value was for the 21-nt class sRNAs in both the NegS (9,598 reads) and PosS (6331 reads) groups. However, the average value in the NegS group (Figures 1 and 2) was the highest for the 22-nt class (13,261 reads), and that in the PosS group was the highest for the 21-nt class (15,088 reads) sRNAs ( Table 2).

Conserved miRNAs Identified in Nasopharynx Tissues
A total of 943 conserved miRNAs were identified in nasopharynx tissues (Supplementary Materials Tables S1 and S2). Among them, five conserved miRNAs were selected for detailed evaluation based on small RNA sequencing results (Figures 3 and 4). The distribution of conserved miRNAs within the NegS and PosS groups is depicted in Figure 5

Conserved miRNAs Identified in Nasopharynx Tissues
A total of 943 conserved miRNAs were identified in nasopharynx tissues (Supplementary  Materials Tables S1 and S2). Among them, five conserved miRNAs were selected for detailed evaluation based on small RNA sequencing results (Figures 3 and 4). The distribution of conserved miRNAs within the NegS and PosS groups is depicted in Figure 5 and Supplementary Materials Table S3. The most highly expressed miRNA in the NegS group was miR-148a, and the reads ranged from 3373 (sample 11A) to 38,424 (1A) with a median value of 22,341. The most highly expressed miRNA in the PosS group was miR-148, which ranged from 3209 (F7) to 53,305 reads (F5) with a median value of 7233 ( Figure 6, Supplementary Materials Table S4).

Conserved miRNAs Identified in Nasopharynx Tissues
A total of 943 conserved miRNAs were identified in nasopharynx tissues (Supplementary Materials Tables S1 and S2). Among them, five conserved miRNAs were selected for detailed evaluation based on small RNA sequencing results (Figures 3 and 4). The distribution of conserved miRNAs within the NegS and PosS groups is depicted in Figure 5 and Supplementary Materials Table S3. The most highly expressed miRNA in the NegS group was miR-148a, and the reads ranged from 3373 (sample 11A) to 38,424 (1A) with a median value of 22,341. The most highly expressed miRNA in the PosS group was miR-148, which ranged from 3209 (F7) to 53,305 reads (F5) with a median value of 7233 ( Figure 6, Supplementary Materials Table S4).        A test for normal distribution showed that the obtained numbers of reads did not fit a normal distribution of the residuals. We used a nonparametric method, the Kruskal-Wallis test, and the p-values of the most abundant miRNA datasets ranged from 0.0113 (miR-100) to 0.04125 (miR-29a). Statistically significant differences were determined between NegS and PosS for miR-100, miR-34b, miR-200a, miR-34c, mir-342, let-7i, and miR-29a (Table 3, Figure 7), where the miRNAs were upregulated in samples of the PosS group. The abundance of conserved mir-148a and mir-21 was not significantly different between NegS and PosS. Table 3. Kruskal-Wallis test emphasizing statistically significant differences between the NegS and PosS groups of five selected miRNAs. All miRNA sequences with total abundance >2000 reads are included. The normalized read (RPM) abundances are listed. Significantly different expressions (p ≤ 0.05) are indicated with an asterisk. A test for normal distribution showed that the obtained numbers of reads did not fit a normal distribution of the residuals. We used a nonparametric method, the Kruskal-Wallis test, and the p-values of the most abundant miRNA datasets ranged from 0.0113 (miR-100) to 0.04125 (miR-29a). Statistically significant differences were determined between NegS and PosS for miR-100, miR-34b, miR-200a, miR-34c, mir-342, let-7i, and miR-29a (Table 3, Figure 7), where the miRNAs were upregulated in samples of the PosS group. The abundance of conserved mir-148a and mir-21 was not significantly different between NegS and PosS.

miRNA Expression Measured by Real-Time RT-PCR
Analysis of miRNA expressions normalized to mir-148a miRNA revealed that three miRNAs (miR-21, miR-34b, and miR-342) were upregulated in the PosS group (Table 4, Figure 8). Of these miRNAs, only miR-21 showed significantly higher expression in the PosS group than in the NegS group. One miRNA, miR-34c, was downregulated in the PosS group; however, there was no significant difference between the PosS and NegS groups in the expression of this miRNA (Table 4, Figure 8). Correlation analysis revealed that the expressions of mir-21 and mir-342 were in strong and very strong, respectively, positive correlation with viral relative quantity ( Table 4). The expressions of mir-34c and mir-34b showed only very weak and moderate, respectively, correlation with the viral relative quantity. Genes 2022, 13, x FOR PEER REVIEW 10 of 17

miRNA Expression Measured by Real-Time RT-PCR
Analysis of miRNA expressions normalized to mir-148a miRNA revealed that three miRNAs (miR-21, miR-34b, and miR-342) were upregulated in the PosS group (Table 4, Figure 8). Of these miRNAs, only miR-21 showed significantly higher expression in the PosS group than in the NegS group. One miRNA, miR-34c, was downregulated in the PosS group; however, there was no significant difference between the PosS and NegS groups in the expression of this miRNA (Table 4, Figure 8).   Correlation analysis revealed that the expressions of mir-21 and mir-342 were in strong and very strong, respectively, positive correlation with viral relative quantity ( Table 4). The expressions of mir-34c and mir-34b showed only very weak and moderate, respectively, correlation with the viral relative quantity.

Discussion
This is the first attempt to use a small RNA high-throughput sequencing technique to identify conserved miRNAs differentially expressed in human nasopharynx tissues in response to SARS-CoV-2 infection. The experimental strategy of this study was designed to investigate the profile of human miRNAs in two groups with 10 samples that tested negative (NegS) and positive (PosS) for SARS-CoV-2 by real-time RT-PCR. Human-pathogen interactions were completely evaluated at the level of total extracted RNA with no additional information concerning the examined individuals. miRNA profiling and prediction associated with SARS-CoV-2 in positive human RNA was performed in previous studies by Arisan et al. [23], Chen et al. [24], Chow and Salmena [25], Demongeot and Seligman [26], Hosseini Rad and McLellan [27], and Sardar et al. [28]; however, these studies used theoretical prediction based on the available datasets without performing

Discussion
This is the first attempt to use a small RNA high-throughput sequencing technique to identify conserved miRNAs differentially expressed in human nasopharynx tissues in response to SARS-CoV-2 infection. The experimental strategy of this study was designed to investigate the profile of human miRNAs in two groups with 10 samples that tested negative (NegS) and positive (PosS) for SARS-CoV-2 by real-time RT-PCR. Human-pathogen interactions were completely evaluated at the level of total extracted RNA with no additional information concerning the examined individuals. miRNA profiling and prediction associated with SARS-CoV-2 in positive human RNA was performed in previous studies by Arisan et al. [23], Chen et al. [24], Chow and Salmena [25], Demongeot and Seligman [26], Hosseini Rad and McLellan [27], and Sardar et al. [28]; however, these studies used theoretical prediction based on the available datasets without performing practical sequencing experiments. The study of Li et al. [29] first used a high-throughput sequencing approach for profiling the miRNAs in the peripheral blood from patients with SARS-CoV-2 infection. Lu et al. [30] used the real-time RT-PCR method to quantify SARS-CoV-2-associated miRNAs in mouse cardiomyocytes targeting the SARS-CoV-2 entry receptor ACE2.
Based on the hypothesis that miRNAs associated with SARS-CoV-2 infection are present in nasopharynx tissues, we identified 943 conserved miRNAs. We revealed that the most abundant small RNAs were the 21-(PosS) and 22-nt classes (NegS). Fang et al. [31] mentioned that the majority of mature human miRNA sequences consist of 22 nucleotides. We revealed that the reads of the 22-nt miRNA class were not normally distributed among NegS and PosS, and the values were not significantly different (p-value 0.7624). The reads of the 21-nt class did not show any significant difference, similar to the 22-nt class (p-value 0.4057). The overexpression of 21-nt miRNA class reads in PosS could be associated with the disease course. According to published results, duplexes of 21 nucleotide short interference RNA (siRNA) with 2-nucleotide in the 3 -overhangs are the most efficient triggers of nucleotide sequence-specific mRNA degradation [32,33].
According to the differences in abundances between the NegS and PosS groups, we focused on the five following conserved miRNAs: miR-148a, miR-21, miR-34c, miR-34b, and miR-342. The most abundant miRNA was miR-148a in the NegS group and in the PosS group, with no significant difference, but a higher expression was observed in the NegS group. MiR-148a in human tissues plays a key role in many biological processes. MiR-148a has ordinary functions shared by many miRNA classes, including a role in cellular differentiation and development. Porstner et al. [34] proved that miR-148a expression increased after pre-B cell activation. Moreover, upregulated miR-148a expression helps to the differentiation of activated B cells to plasma cells and so it helps the survival of plasma cells by constraining various transcription and proapoptotic factors. This miRNA is a muscle-derived miRNA that may facilitate the differentiation of myoblasts and skeletal muscle cells by targeting regulatory Rho-associated coiled-coil containing protein kinase 1. MiR-148a may also promote the growth of myoblasts in the G1 phase of the cell cycle and shorten the S phase. This role prompts myoblast differentiation into myotubes [35]. Van Wijnen et al. [36] showed that miR-148a is involved in osteoclast formation. MiR-148a induces the transformation of monocytes to osteoclasts by inhibition of the transcription factor V-maf musculoaponeurotic fibrosarcoma oncogene homolog B. Gailhouste et al. [37] showed that enhanced expression of miR-148a might induce the differentiation and maturation of liver cells by inhibiting DNA. MiR-148a also can promote primary adipocytes to differentiate into mature adipocytes [38]. MiR-148a may affect the development of the nervous system by targeted regulation too [39]. It has been shown that miR-148a regulates several phenotypes, including those present in embryonic stem cells [40]. Downregulated expression of miR-148a can control the phenotype of mesenchymal stem cells. This is done by helping the expression of the endothelial PAS domain that contains protein 1 transforming factor [41]. The downregulated expression of miR-148a can be detected in many types of cancers, including gastric, colorectal, pancreatic, liver, oesophageal, breast, non-small cell lung, and urogenital system cancers. Nonetheless, upregulated expression of miR-148a may be seen also in glioma and osteosarcoma. Moreover, the expression levels of miR-148a have been clearly and strongly linked to the clinical classification, efficacy, and prognosis of cancer [42].
Second, the most abundant miRNA in both investigated groups was miR-21. miR-21 was significantly differentially expressed between the NegS and PosS groups according to real-time RT-PCR (p-value 0.03156). In contrast to miR-148, based on the results miR-21 probably has a direct link to SARS-CoV-2 infection. No significant difference (p-value 0.1306) was observed between the NegS and PosS groups by small RNA sequencing, but higher expression of miR-21 was observed in the PosS group according to both detection methods (Tables 3 and 4, Figures 4-6). The scatter plot ( Figure 6) also showed that higher abundances of miR-21 were typical for the PosS group. This means that the trend is the same as was found out by real-time RT-PCR. Moreover, a strong positive correlation between miR-21 expression and COVID-19 titer (Pearson's correlation coefficient 0.73, Table 4) indicates a tight connection between miR-21 expression and infection with COVID-19. Farr et al. [43] have recently posted the same results, showing that miR-21 was upregulated in COVID-19. miR-21, a SARS-CoV-2-binding microRNA, has confirmed four binding sites on the SARS-CoV-2 genome. mir-21 is one of the better-known miRNAs whose expression increases in many pathological conditions, including asthma, pulmonary fibrosis, and viral infection [44,45]. There are no reports on the direct binding of miR-21 to other human viral genomes in recent years, and current reports about the involvement of miR-21 in viral infections are very limited to host transcripts modulating. For example, the positive role of miR-21 in influenza A replication was attributed to the miR-21-host HDAC8 interaction [46]. In addition, it has been shown that miR-21 reduces the antiviral NF-KB pathway by binding to IRAK1 and TRAF6 transcripts in HIV and HCV infections [47,48]. According to the results of Jafarinejad-Farsangi et al. [49], miR-21 has two binding sites in the spike protein-coding regions. In addition, miR-21 was one of the top miRNAs that targeted the upregulated host DEGs in response to SARS-CoV-2, which agrees with our results, and it showed normalized median reads of NegS (3531) versus PosS (5772). One of the miR-21 targets is CXCL-10, which is a biomarker for viral, bacterial, fungal, and parasitic contamination [50].
In the study of Jafarinejad-Farsangi et al. [49], high levels of CXCL-10 were observed in the lungs of COVID-19-infected patients compared to healthy patients. miR-21 expression was increased in the COVID-19 group compared to healthy controls [51].
The rest of the comprehensively quantified miRNAs in this study showed significant differences between the NegS and PosS groups according to small RNA sequencing. The third most expressed miRNA was mir-34c. The number of normalized miR-34c reads was significantly different according to a p-value of 0.02334, and miR-34c was upregulated in PosS (3350 reads) compared to NegS (1444 reads); however, this phenomenon was not confirmed by real-time RT-PCR. This is only one miRNA where the trend is different between small RNA sequencing and real-time RT-PCR results. The results of small RNA sequencing are in contrast with a case-control study [52] in which respiratory syncytial virus (family Paramyxoviridae)-infected patients showed low levels of miR-34c expression compared with controls. qPCR is often considered a gold standard in the detection and quantitation of various gene expressions. Nevertheless, the rapid increase in the number of miRNAs renders qPCR inefficient on a genomic scale, and it is probably better used as a validation rather than a discovery tool [53]. Small RNA sequencing comes into wider use and is unmatched for the discovery and experimental validation of novel or predicted miRNAs [53]. miR-34c is a 77 bp long noncoding RNA. This miRNA is located on human chromosome 11 belonging to the miR-34 family. The miR-34 family may regulate cell processes by binding to target gene sequence fragments. miR-34c inhibits Bcl2 by binding to the 3 untranslated region (UTR) of the Bcl2 gene, thus downregulating the viability of laryngeal cancer cells and inducing apoptosis [54]. miR-34c is also associated with emphysema severity and thus modulates SERPINE1 expression [55]. SERPINE1 encodes a member of the serine proteinase inhibitor superfamily. SERPINE1 is the principal inhibitor of tissue plasminogen activator and urokinase. The final protein also functions as a component of innate antiviral immunity and high concentrations of the gene product are associated with thrombophilia [56].
miR-34b was detected by small RNA sequencing by 49 reads in the NegS group and by 124 reads in the PosS group, and the p-value was 0.01258. The results were confirmed by realtime RT-PCR but not with significant values (p-value 0.22542, Pearson's correlation coefficient 0.48). In total, we determined seven miRNAs with statistically significant expression level differences (miR-100, miR-34b, miR-200a, miR-34c, miR-342, let-7i, and miR-29a). The results are not in line with the current study by Demiray et al. [57]. Statistically significant expression level differences (p < 0.05) were detected by real-time RT-PCR in nine miRNAs in COVID-19 patients and healthy controls. Seven miRNA (let-7d, miR-17, miR-34b, miR-93, miR-200b, miR-200c, and miR-223) expression levels were found to be significantly decreased, and the expression levels of two miRNAs (miR-190a and miR-203) were significantly increased compared to healthy controls [57]. In total, the findings that have been captured by Demiray et al. [57] can be validated by our experiment only in the case of miR-203 and miRNAs let-7d, miR-17, miR-34b, miR-93, miR-200b, miR-200c, miR-223, and miR-190a showed the opposite expression. Moreover, it is not clear why Demiray et al. [57] focused on 21 miRNAs expressed in serum and how selected these 21 miRNAs.
The lowest level of expression of the selected miRNAs was observed for miR-342; in the NegS group, it was detected by 18 reads, and in the PosS group, it was detected by 55 reads (p-value 0.0411). Its upregulation was also observed by real-time RT-PCR, and the relative expression values were 0.042 (NegS) and 0.065 (PosS group) with a significance p-value of 0.24653. Interestingly, there was a very strong positive correlation between miR-342 expression and the COVID-19 titer (Pearson's correlation coefficient 0.87, Table 4). However, the correlation result may not be relevant as, although significantly upregulated in the PosS group unveiled by small RNA sequencing, the miR-342 expression did not differ in PosS and NegS in realtime RT-PCR. The computational approach revealed an exciting hypothesis that miR-342 is involved in MERS-CoV pathogenesis [58]. mir-342-5p also suppresses coxsackievirus B3 biosynthesis by targeting the 2C-coding region [59].
The identified miRNAs were expressed in both RNA groups, NegS and PosS: miR-148a, miR-21, miR-34c, miR-34b, and miR-342. The most differentially expressed miRNA was miR-21, which is likely closely linked to the presence of SARS-CoV-2 in nasopharynx tissues. Upregulation of miR-21 in the PosS group was almost two times higher and could also be linked with cardiac fibroblast and endothelial cell dysfunction in COVID-19 patients [51]. Due to the nature of our results, it is probable that a rapid increase in the number of miRNAs renders qPCR inefficient on a genomic scale, and it is maybe better to use small RNA sequencing as a discovery tool. To our knowledge, this is the first study reporting important miRNAs detected in SARS-CoV-2-positive RNA extracted from nasopharyngeal tissues. Although this study encompasses the role of RNAs, a description of the symptoms caused by SARS-CoV-2 is missing and future studies should include this information. These findings should serve as a foundation for larger studies and should contribute to evaluating the long-term course of COVID-19 patients.