Isolation and characterization of plasma-derived extracellular vesicles
Plasma EVs were characterized for their size and morphology using TEM. Particles are mostly in exosomal range nonetheless large sized vesicle populations were also observed. Further, nanoparticle tracking analysis did not show any significant difference either in the size range or concentrations across all the groups observed (Fig. 1 B-C). Immunoblot analysis showed the presence of exosomal markers CD63 and CD81 and negative for calnexin and any endoplasmic reticulum contamination.
Input read alignment and small RNA biotype mapping
The reading obtained from sequencing were used for alignment and mapping to the human genome after clipping and quality filtering. The input reads in plasma exosomes were between 5.58 million reads to 21.75 million reads between NS (14031359 ± 2019014), CS (10688617 ± 917325), WP (8765589 ± 768710), E-cig users (9570944 ± 455585) and DS (8988900 ± 798445). Reads were mapped to human rRNA to exclude rRNA sequences before mapping to human genome. The percentage of input reads alignment from each subject in individual group was given in Fig. 2. The input reads were significantly lower in WP and DS in comparison to NS subjects (P<0.05).
The reads mapped to human genome were classified into various small RNA biotypes. After mapping and excluding rRNA, the microRNAs mapped along with all other RNA transcripts in GENCODE. The microRNA was between 78-81% of all biotype counts in all NS, CS, WP, E-cig and DS groups. There was significant lower counts of Mt-tRNA observed in CS, WP, E-cig and DS (P<0.05) in comparison to NS. Further, snoRNA counts were also significantly (P<0.05) lower in all groups in comparison to NS (Fig. 3 A, B).
Comparison of microRNAs expression profiles between non-smokers, E-cig users, cigarette smokers, waterpipe smokers and dual smokers
Multidimensional Scaling (MDS) plot
The MDS plot was generated using the microRNAs that exhibited largest distances across all the samples from subjects of the NS, CS, WP, E-cig and DS groups. Majority of the NS samples cluster separately whereas all the smokers (CS, WP and DS), and E-cig users cluster separately suggesting these effects may be due to common tobacco use (Fig. 4).
Volcano plot
Majority of the differentially expressed microRNA from NS vs. CS, NS vs. WP, NS vs. E-cig users, NS vs. DS (Fig. 5 A-D) and CS vs. E-cig, CS vs. WP, CS vs. DS, WP vs. DS (Suppl. Fig. 1 A-D) are presented as volcano plot. The volcano plots were made by plotting the -log10 of adjusted p-values on the y-axis, and the log2 fold change between two groups on the x-axis showing the up- and down regulation appearance similar distance from the center. Some microRNAs with extreme log2 fold changes were not included in the volcano plot in order to make the scales of the volcano plots consistent among different comparisons. The microRNA values plotted show two regions with highest magnitude of fold change and high statistical significance. The highlighted spots (red color) with microRNA names are with greatest difference (at least two fold changes) in expression and statistically significant (adjusted P<0.05) after correction for multiple testing.
Hierarchical Clustering of microRNAs
The hierarchical cluster analyses of differentially expressed microRNAs were performed in all the groups. The heat maps were generated using normalized values from individual samples from NS vs. CS, NS vs. WP, NS vs. E-cig users, NS vs. DS (Fig. 6 A-D) and CS vs. E-cig, CS vs. WP, CS vs. DS, WP vs. DS (Suppl. Fig. 2 A-D). In the heatmap, each row represents individual microRNA and each column an individual sample. The microRNA clustering on the left indicates hierarchical clustering of significant microRNA. The color scale at the right side in panel A indicates the relative expression level of microRNA in all samples. The red color indicates lower level than the mean and green a level higher than the mean.
Differentially expressed microRNAs in E-cig versus non-smokers group
The microRNAs that were most differentially expressed in plasma exosomes of E-cig users in comparison to NS are presented in Suppl. Table 1. There were a total of 17 microRNAs that were differentially altered in E-cig users in comparison to NS. Among them 4 were downregulated and 13 were upregulated as shown in Suppl. Table 1. The upregulated microRNAs include hsa-miR-365a-3, hsa-miR-365b-3p, hsa-let-7f-5p, hsa-miR-1299, hsa-miR-21-5p, hsa-let-7i-5p, has-let-7a-5p, hsa-miR-30a-5p, hsa-miR-193b-3p, hsa-miR-100-5p, hsa-miR-423-3p, hsa-miR-30c-5p, hsa-miR-143-3p and hsa-miR-224-5p,. The downregulated expression of 4 microRNAs includes hsa-miR-362-5p, hsa-miR-29b-3p, hsa-miR-451a and hsa-miR-30e-5p.
Differentially expressed microRNAs in cigarette smokers versus non-smokers group
The differentially expressed microRNAs that were significantly altered between CS versus NS were given in Suppl. Table 2. A total of 24 microRNAs were significantly altered of which 16 were upregulated and 8 were downregulated as shown in Suppl. Table 2. The maximum fold change in upregulated microRNAs are hsa-mir-149-5P (20.29 fold), hsa-miR-532-5p (19.79) andhsa-miR-2355-5p (19.65), while downregulated include hsa-miR-29b-3p (23.57), hsa-miR-29a-3p (2.58), and hsa-miR-320b (-2.13).
Differentially expressed microRNAs in waterpipe smokers versus non-smoker group
Sixteen miRNAs were differentially expressed in WP versus NS comparisons as shown in Suppl. Table 3. The maximum fold change observed in upregulated microRNAs are hsa-miR-2355-5p (39.81 fold), hsa-miR-149-5p (29.22) and hsa-miR-582-5p (23.15and downregulated are hsa-miR-362-5p (45.35), hsa-miR-29b-3p (24.62) and hsa-miR-320d (-4.57).
Differentially expressed microRNAs in dual smokers versus non-smokers
The DS vs NS groups have shown total 20 differential expressed microRNAs that are shown in Suppl. Table 4. The top fold change upregulated microRNAs are hsa-miR-149-5p (29.29 fold), hsa-miR-139-5p (16.41), hsa-miR-424-3p (16.01), and downregulated microRNAs are hsa-miR-362-5p (44.11 fold), hsa-miR-29b-3p (21.54), and hsa-miR-144-3p (2.36).
Differentially expressed microRNAs in E-cig users versus cigarette smokers
Total 9 microRNAs were altered significantly in E-cig users versus CS of which 5 were upregulated and 4 were downregulated (Supplementary Table 5). Among the top fold changed microRNAs, upregulated in E-cig are hsa-miR-362-5p (19.67), hsa-miR-2355-5p (19.62), hsa-miR-532-5p (19.41), and downregulated microRNAs are hsa-miR-365a-3p (-24.12), hsa-miR-1299 (-24.01), and hsa-miR-193b-3p (-7.85).
Differentially expressed microRNAs in waterpipe smokers versus cigarette smokers
Differentially expressed microRNAs in WP versus CS are given in Supplementary Table 6. Out of six microRNAs differentially expressed, the upregulated are hsa-miR-532-5p (21.30 fold), hsa-miR-362-5p (20.47), hsa-miR-144-5p (19.42), and downregulated are hsa-miR-1299 (-23.53), hsa-miR-582-5p (-22.63), hsa-miR-1-3p (-6.72) in WP compared to CS.
Differentially expressed microRNAs in dual smokers versus cigarette smokers
The 5 microRNAs that were differentially expressed in DS versus CS are shown in Supplementary Table 7. hsa-miR-424-3p (20.26) is upregulated, while hsa-miR-144-5p (21.40), hsa-miR-532-5p (20.74), hsa-miR-2355-5p (19.26), and hsa-miR-362-5p (19.22) were downregulated in CS compared to DS .
Differentially expressed microRNAs in waterpipe smokers versus dual smokers
The 4 microRNAs differentially downregulated in WP compared to CS are presented with fold changes [hsa-miR-2355-5p (39.43), hsa-miR-1299 (21.29), hsa-miR-582-5p (21.25), hsa-miR-1-3p (7.29) and upregulated microRNAs are hsa-miR-139-5p (22.07), hsa-miR-424-3p (21.71) in Supplementary Table 8.
Overlap of microRNA expression in plasma exosomes
The overlaps of microRNAs expression between all four groups are presented as Venn diagram in Fig. 7 A-C. We have compared microRNAs expressed in all four groups: NS vs. CS, NS vs. WP, NS vs. E-cig and NS vs. DS. These groups have common expression of7 microRNAs that are hsa-let-7a-5p, hsa-miR-21-5p, hsa-miR-29b-3p, hsa-let-7f-5p, hsa-miR-143-3p, hsa-miR-30a-5pand hsa-let-7i-5p. The E-cig group has expressed 5 microRNAs (hsa-miR-224-5p, hsa-miR-193b-3p, hsa-miR-30e-5p, hsa-miR-423-3p, and hsa-miR-365a-3p|hsa-miR-365b-3p) that are specific for this group, not expressed in other three groups (Fig. 7 A). The comparison of up-regulated microRNAs (NS vs. other smoking/vaping groups) in all four groups revealed expression 6 microRNA (hsa-let-7a-5p, hsa-miR-21-5p, hsa-let-7i-5p, hsa-let-7f-5p, hsa-miR-143-3p and hsa-miR-30a-5p) common to all groups. However, there are 5 microRNA (hsa-miR-224-5p, hsa-miR-423-3p, hsa-let-7c-5p andhsa-miR-365a-3p|hsa-miR-365b-3p, and hsa-miR-193b-3p) expressed specifically to E-cig group (Fig. 7 B). When NS vs. other groups were compared for down-regulated microRNAs, hsa-mir-29b-3p was found to be common in all the groups. The microRNA expressed specifically in E-cig was hsa-mir-30e-5p and was downregulated. (Fig. 7 C).
Receiver operating characteristic (ROC) curve for evaluation of diagnostic utility
A ROC curve was built using one of the identified significant marker miRNA hsa-let-7a-5p to distinguish normal subjects from tobacco users (Fig. 8). The ROC curve showed that both the sensitivity and specificity of this miRNA marker are greater than 80% according to the optimal Youden Index point, which is an indicator of the performance of miRNA hsa-let-7a-5p to distinguish normal subjects from tobacco users.
Gene enrichment analysis of differentially expressed microRNAs
The FunRich enrichment analysis of differentially expressed microRNAs was performed to explore the potential target genes in NS vs. E-cig, NS vs. CS, NS vs. WP and NS vs DS pairwise comparisons (Fig. 9 A-F). The top six enriched functions with the lowest p values were biological pathway, biological process, molecular function, cellular component, site of expression and transcription factor in all groups. The top 3 biological pathway with the lowest p values were beta1 integrin cell surface interactions, integrin family cell surface interactions, and TRAIL signaling pathway common in NS vs. CS, NS vs. WP, NS vs. E-cig and NS vs. DS. The proteoglycan-mediated signaling events changed significantly in all three groups except NS vs. E-cig. In addition, endothelin biological pathways with lowest p values were in NS vs. E-cig and NS vs. DS.
The biological process in regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolism with lowest p value were present in all four groups. The two molecular functions with highly significant values related to transcription factor activity and extracellular matrix structural constituent were common in all four groups. The top two cellular components related to nucleus and cytoplasm were common in all four groups. The three site of expression of microRNAs with lowest p values were kidney, placenta and skeletal muscle. However, the other sites of expression related to lung with significant p values were in NS vs. CS, NS vs. E-cig and NS vs. DS. The transcription factor related EGR1, SP1, SP4 and POU2F1 were highly significant in all four groups while ZFP161 in only NS vs. E-cig and NS vs. DS. Further, pairwise comparisons were done in between E-cig vs. CS, WP vs. CS, DS vs. CS and WP vs. DS groups were presented in Supplementary Figure 3 A-F.
Target genes of plasma exosomal microRNA changes
The list of target genes of plasma exosomal microRNA changes observed in NS vs. E-cig, NS vs. CS, NS vs. WP, NS vs DS and CS vs. E-cig, CS vs. WP, CS vs. DS, WP vs. DS are presented in Supplementary Table 9.
Differential tRNA fragments identified in plasma exosomes
The changes in tRNAs were calculated based on the trimmed mean of M values (TMM) method using normalized tRNA counts in all NS, CS, WP, E-cig and DS groups. The read counts normalized and subjected to differential expression analysis by DESeq2, to get the change in tRNA expression between different groups by pairwise comparisons. The read raw count data analysis showed 25 different types of tRNAs in plasma exosomes of all groups NS vs. CS, NS vs. WP, NS vs. E-cig and NS vs. DS. The pairwise comparison data revealed significant changes in 7 tRNAs in all NS vs. CS, NS vs. WP and NS vs. E-cig groups. However, NS vs. DS group showed changes in eight tRNAs. All four groups have significant increase in six tRNAs (tRNAVal, tRNAGlu, tRNAAsp, tRNAGly, tRNAArg and tRNAHis) and decrease in tRNACys. In addition, NS vs. DS group also showed significant increase in tRNAIle (Supplementary Table 10-13). There was no significant change observed in tRNAs counts (total 24 raw counts in each) in all CS vs. DS, CS vs. WP, CS vs. E-cig and WP vs. DS.
The MDS plots were generated using the tRNAs that exhibited largest distance across all the samples from subjects of the NS, CS, WP, E-cig and DU groups. Majority of the samples from all the groups (NS, CS, WP, DS and E-cig) did not cluster together (Suppl. Fig. 4).
The differentially expressed tRNAs from NS vs. CS, NS vs. WP, NS vs. E-cig users, and NS vs. DS (Suppl. Fig. 5 A-D) are presented as volcano plot. The p-values on the y-axis, and the fold change between two groups on the x-axis are plotted showing the up- and down regulation appearance similar distance from the center. The tRNAs values plotted show two regions with highest magnitude of fold change and high statistical significance. The highlighted spots (red color) with tRNAs names are with greatest difference in expression and statistically significant (P<0.05) after correction for multiple measurement.
The hierarchical cluster analyses of differentially expressed tRNAs were done in all different groups. The heat maps were generated using normalized values from individual sample from NS vs. CS, NS vs. WP, NS vs. E-cig users, NS vs. DS (Fig. 6 A-D). In the heat map, each row represents individual tRNA and each column individual sample. The tRNA clustering on the left indicates hierarchical clustering of significant tRNA. The color scale at the right side indicates the relative expression level of tRNA in all samples. The red color indicates lower level than the mean and green a level higher than the mean.
The overlap of tRNAs expression between all four groups are presented as Venn diagram in Supp. Fig. 7. We have compared tRNAs expressed in all four groups: NS vs. CS, NS vs. WP, NS vs. E-cig and NS vs. DS. All the four groups have common changes in 7 tRNAs (tRNAVal, tRNAGlu, tRNAAsp, tRNAGly, tRNAArg and tRNAHis, and tRNACys). However, change in tRNAIle was only observed in NS vs. DS.
Differential expression of piRNA in plasma exosomes
Read counts of piRNAs from NS, CS, WP, E-cig, and DS groups were TMM-normalized and normalized counts were used to generate a MDS plot (Supp. Fig. 8). There was no close clustering observed between individual samples of these groups. Normalized counts were also were processed for differential expression analysis by DESeq2 to get the data of different group pairwise comparisons.
The hierarchical cluster analyses of differentially expressed piRNAs were done in all different groups. The heatmaps were generated using normalized values from individual sample from NS vs. CS, NS vs. WP, NS vs. E-cig users, NS vs. DS (Supp. Fig. 9 A-D). In the heat map, each row represents individual piRNA and each column individual sample. The piRNA clustering on the left indicates hierarchical clustering of significant piRNA. The color scale at the right side in indicates the relative expression level of piRNA in all samples. The red color indicates lower level than the mean and green a level higher than the mean. The significant piRNAs across all the pairwise comparisons were given in the supplementary Table 14.