Peripheral blood was collected in 10-mL heparin anticoagulation tubes from eight volunteers at Wuhan Sixth Hospital. All volunteers provided written informed consent. Sterile phosphate buffer (pH 7.4) was diluted in an equal volume and centrifuged at 3000 × g for 5 min. A pipette was used to aspirate the supernatant, which was rapidly frozen in liquid nitrogen for 1 h and then incubated at −80°C or stored in liquid nitrogen for long-term storage.

For RNA extraction, TRIzol was added to the blood samples at a 1:3 ratio and thoroughly mixed. Thereafter, the sample was thoroughly blown with a pipette and shaken until complete dissolution of the floc, followed by incubation at room temperature for 5 min to completely degrade the nucleoprotein body.

The RNA samples were quantified using various methods. Nanodrop was used for determining RNA purity (OD₂₆₀/₂₈₀ ratio) and Qubit for precisely quantifying RNA concentration. The Agilent 2100 bioanalyzer was used for precisely detecting RNA integrity and indicators such as RNA integrity number, 28S/18S, presence or absence of uplift at the baseline of the graph, and 5S peak. RNA electrophoresis was used to detect whether the sample was diffuse or contaminated with genomic DNA.

After passing the sample testing, the Small RNA Sample Pre Kit was used to construct libraries using the special structures of the 3′ and 5′ ends of the small RNA (sRNA), i.e., the complete phosphate group at the 5′ end and hydroxyl group at the 3′ end. Using total RNA as the starting sample, sRNAs were directly synthesized by adding connectors to both ends. Thereafter, reverse transcription was performed for cDNA synthesis. After polymerase chain reaction (PCR) amplification, the target DNA fragments were separated via polyacrylamide gel electrophoresis. The cDNA libraries were obtained by cutting and recovering the amplicons from the gel. The quality of the constructed libraries was tested, with all requirements met before the library was sequenced. The following instruments were used for quality testing: 1) Qubit 2.0 fluorometer for accurate quantification and 2) the Agilent 2100 bioanalyzer for library size detection. After the library size met the required expectations and the libraries passed the tests, they were sequenced using the Illumina platform based on the effective concentration and target downstream data volume.

Eight miRNA transcriptome sequencing libraries were constructed by the Frasergen Company for the healthy and PH groups. The quality of the raw sequencing data was evaluated, and the raw data (raw reads) were filtered using Seqtk software to obtain high-quality data (clean reads). These clean reads were compared with the reads in miRBase, a human genome database used to identify known miRNAs. The unknown sRNAs were compared with the RNAs in Rfam database, a database comprising noncoding RNAs (ncRNAs). Thereafter, the sRNA fragments were annotated and counted, and the unknown small rRNAs were predicted using miRDeep software based on the signature hairpin structures of the miRNA precursors.

The following analysis was performed by taking the average of the two sample groups, which were divided into the CK (healthy) and PH groups for comparison. The reads from each sample were compared using miRBase to determine miRNA expression (in counts per million [CPM]). The expression was screened using CPM to calculate the metrics (CPM). The DESeq software was used to perform differential expression analysis between the PH and CK groups. miRNAs with a p-value of ≤0.05 and log2 (foldchange [FC]) of ≥1 were selected as the differentially expressed miRNAs. The transcriptome data has been submitted to the National Genomics Data Center (https://ngdc.cncb.ac.cn/; GSA No.: HRA007288).

The differentially expressed miRNAs were analyzed. First, to visualize the distribution of the false discovery rate and differential FC of the differential miRNAs between the two sample groups, a volcano plot comparing the two sample groups was constructed. Second, to visualize the expression patterns of the differential miRNAs in all samples, expression pattern and clustering analyses were performed.

The online prediction software DAVID bioinformatics Resources 6.7 (https://david.ncifcrf.gov/) was used to functionally predict the target genes of the miRNAs, including GO terms and KEGG pathway enrichment analysis. The predicted target genes were enriched using GO and KEGG to analyze the functions and signaling pathways of the target genes and to elucidate the potential functions of the differentially expressed miRNAs (Ashburner et al., 2000; Kanehisa et al., 2004).

To verify the accuracy of transcriptome sequencing, four differentially expressed miRNAs in the PH group compared with the CK group were detected using qRT-PCR to verify the relative expression. 1) The upstream primers for qRT-PCR were designed based on the mature sequences of the four miRNAs, whereas the downstream primers were provided by the miRcute Plus miRNA qPCR Kit (Tiangen FP411). The qRT-PCR primer sequences were designed according to the instructions of the reagent kit and synthesized by Shanghai Bioengineering Biotechnology Service Co. 2) cDNA was prepared using 2× miRNA RT reaction, 10 μL buffer, 2 μL of miRNA RT Mix, 500 ng of RNA, and RNase-free water in sterilized Eppendorf tubes to a final volume of 20 μL. The product was reverse transcribed at 42°C for 60 min and 95°C for 3 min and then stored at −20°C. 3) For real-time quantitative analysis of the relative expression of the four miRNAs using the SYBR Green dye method, cDNA was used as the template and U6 was used as the internal reference gene. The total reaction volume for qRT-PCR was 10 μL and comprised the following components: 5 μL of 2× mix, 1 µL of the amplification primer 1 (2 μM), 1 µL of the reverse primer (2 μM), 1 μL of cDNA, and 2 µL of ddH₂O. The mixture was amplified under the following conditions: 50°C for 2 min; 95°C for 10 min; 40 cycles at 95°C for 15 s and 60°C for 1 min; and then 95°C for 15 s, 60°C for 15 s, and 95°C for 15 s. Three replicates were used for each sample. The relative expression of the target genes was calculated using the 2^−ΔΔCT method, where ΔCt = Ct target gene − Ct internal reference gene, CK-1 was used as the control sample, and ΔΔCt = ΔCt experimental sample −ΔCt control sample. The sequencing results were compared using Excel.

The correlation of miRNA expression between samples is an important index to assess the reliability of an experiment and whether the sample selection is reasonable. The closer the correlation coefficient is to 1, the higher the similarity of expression patterns between samples. The correlation coefficients of intra- and inter-group samples are calculated based on the expression of all genes in each sample, and a heat map is drawn, which can directly demonstrate the sample differences between groups and the repetition of the samples within groups. The correlation between the miRNA expression in healthy controls and patients with PH is shown in Figure 1A. The correlation between the samples in the groups was good; however, the correlation between the groups was lower than that between groups, indicating that the data can used for differential miRNA analysis.

Determining TPM is a method for analyzing the expression of miRNA. It can determine the amount of gene expression by calculating the TPM value of each miRNA. TPM value refers to the percentage of each miRNA in the total RNA. It can reflect gene expression and is its commonly used unit. The advantage of the TPM method is that it can eliminate the differences between samples to compare the gene expression of different samples. Additionally, TPM can be used to compare the expression of different genes to determine which genes have the greatest expression in different samples. However, the miRNA TPM and TPM density distribution analysis of the eight samples revealed that the overall differences were not significant (Figures 1B, C).

We compared the clean reads of mature miRNAs, miRNA hairpins, mRNAs, mature and primary tRNAs, small nucleolar RNA (snoRNAs), rRNAs, other ncRNAs, and (optional) known RNAs using bowtie (Langmead and Salzberg, 2012). The ncRNAs were removed as much as possible to enable miRNA detection and prediction. The results showed that 66.27% of the sRNAs sequenced in the CK group were mapped to known miRNA sites in the genome, 14.45% to exons, and 4.27% to the snoRNA loci (Figure 2A). Of the sRNAs sequenced in the PH group, 73.3% were mapped to the known miRNA loci in the genome, 7.62% to exons, and 3.14% to the snoRNA loci (Figure 2D).

miRNA precursors are cleaved by Dicer to form mature miRNAs, and the restriction sites mark the first base of the mature miRNA sequence. Therefore, the distribution of the first base of the sRNAs and position of the mature miRNAs were compared and analyzed (Figures 2B, E). The results showed that the miRNAs in the CK and PH groups had highly similar first base preference.

For the clean reads of each sample, sRNAs of a certain length were selected for further analysis. The results showed that the sRNAs in the CK and PH groups had a highly similar length distribution (Figures 2C, F).

Generally, genes encoding miRNAs whose expression differs by > 2-fold in different groups are differentially expressed. Thus, setting the criteria for screening differential genes is critical. We used |log2 (FC)| >1 and Padj <0.05 as the criteria for identifying differentially expressed genes.

The volcanic map revealed 136 miRNAs with >2-fold differential expression in the PH group, of which 68 were upregulated and 68 were downregulated (Figure 3A). The numbers of upregulated and downregulated miRNAs showing differences between 2- and 4-fold, 4 and 8-fold, and >8-fold in the PH group were 42 and 46, 12 and 11, and 14 and 11, respectively (Figure 3B). A heat map (Figure 4A) and a line graph (Figure 4C) show the pattern of the differentially expressed miRNAs and their clustering in eight samples. Additionally, the difference in miRNA expression between the CK and PH groups was clustered using a heat map (Figure 4B).

We predicted the genes targeted by the differentially expressed miRNAs using MiRanda and qTar, and the genes at the intersection of the two software were shortlisted for GO classification and enrichment analysis (Figure 5). GO analysis (Figure 5A) revealed that the GO terms for the top five differentially expressed genes were signaling, cell communication, positive regulation of biological process, positive regulation of cellular process, and regulation of response to stimulus. Among the three categories of GO enrichment analysis, the top three GO terms for the biological process were cellular process, biological regulation, and regulation of biological process; the two GO terms enriched for the cellular component were cellular anatomical entity and protein-containing complex; and the top three GO terms for molecular function were binding, catalytic anatomical entity, and molecular function regulator (Figure 5B).

KEGG pathway analysis (Figure 6A) revealed that the top five pathways related to the differentially expressed genes were MAPK signaling, Ras signaling, mTOR signaling, cholinergic synapse, and NF-kappa B signaling pathways. Among the five categories of KEGG pathway enrichment analysis, the top three terms for the differentially expressed genes related to cellular processes were transport, catabolism, and cellular community-eukaryotes; the top three terms related to environmental information processing were signal transduction, signaling molecules and interaction, and membrane transduction; the top three terms related to genetic information processing were folding, sorting, and degradation; translation; and replication and repair; the top three terms related to metabolism were lipid, carbohydrate, and amino acid metabolism; and the top three terms related to organismal systems were immune, endocrine, and nervous systems (Figure 6B).

We selected two upregulated miRNAs (hsa-miR-1304-3p and hsa-miR-490-3p) and two downregulated miRNAs (hsa-miR-11400 and hsa-miR-31-5p) (Figure 3A) in the PH group based on their abundant and differential expression multiple (recognizable). We performed qRT-PCR for these four miRNAs to validate our transcriptome data and determine the feasibility of using them as clinical diagnostic markers for PH.

qRT-PCR confirmed that the expression of hsa-miR-1304-3p and hsa-miR-490-3p was significantly upregulated, whereas that of hsa-miR-11400 and hsa-miR-31-5p was significantly downregulated, in the peripheral blood of all patients with PH compared with the CK group (Figure 7). Collectively, these findings prove the reliability and validity of our data and provide evidence that these four differentially expressed miRNAs are potential diagnostic markers for PH.