Exosomal microRNAs are novel circulating biomarkers in cigarette, waterpipe smokers, E-cigarette users and dual smokers

Electronic cigarettes (e-cigs) vaping, cigarette smoke, and waterpipe tobacco smoking are associated with various cardiopulmonary diseases. microRNAs are present in higher concentration in exosomes that play an important role in various physiological and pathological functions. We hypothesized that the non-coding RNAs transcript may serve as susceptibility to disease biomarkers by smoking and vaping. Plasma exosomes/EVs from cigarette smokers, waterpipe smokers and dual smokers (cigarette and waterpipe) were characterized for their size, morphology and TEM, Nanosight and immunoblot analysis. Exosomal RNA was used for small RNA library preparation and the library was quantified using the High Sensitivity DNA Analysis on the Agilent 2100 Bioanalyzer system and sequenced using the Illumina NextSeq 500 and were converted to fastq format for mapping genes. Enrichment of various non-coding RNAs that include microRNAs, tRNAs, piRNAs, snoRNAs, snRNAs, Mt-tRNAs, and other biotypes are shown in exosomes. A comprehensive differential expression analysis of miRNAs, tRNAs and piRNAs showed significant changes across different pairwise comparisons. The seven microRNAs that were common and differentially expressed of when all the smoking and vaping groups were compared with non-smokers (NS) are hsa-let-7a-5p, hsa-miR-21-5p, hsa-miR-29b-3p, hsa-let-7f-5p, hsa-miR-143-3p, hsa-miR-30a-5p and hsa-let-7i-5p. The e-cig vs. NS group has differentially expressed 5 microRNAs (hsa-miR-224-5p, hsa-miR-193b-3p, hsa-miR-30e-5p, hsa-miR-423-3p, hsa-miR-365a-3p, and hsa-miR-365b-3p), which are not expressed in other three groups. Gene set enrichment analysis of microRNAs showed significant changes in the top six enriched functions that consisted of biological pathway, biological process, molecular function, cellular component, site of expression and transcription factor in all the groups. Further, the pairwise comparison of tRNAs and piRNA in all these groups revealed significant changes in their expressions. Plasma exosomes of cigarette smokers, waterpipe smokers, e-cig users and dual smokers have common differential expression of microRNAs which may serve to distinguish smoking and vaping subjects from NS. Among them has-let-7a-5p has high sensitivity and specificity to distinguish NS with the rest of the users, using ROC curve analysis. These findings will pave the way for the utilizing the potential of exosomes/miRNAs as a novel theranostic agents in lung injury and disease caused by tobacco smoking and vaping.


Background
Electronic cigarettes (E-cigs) are battery operated various kind of device, also known as electronic nicotine delivery system. These devices produce aerosolized substances by heating a liquid, which generates a large number of chemicals due to the presence of propylene glycol, vegetable glycerin, nicotine in various concentrations, flavoring agents and other additive compounds [1]. Though E-cig aerosol contains less number of toxic chemicals than cigarette smoke, but it may produce adverse health effects. E-cig aerosol contains ultra-fine particles, heavy metals, volatile organic compounds, and numerous toxic chemicals including acetaldehyde, acrolein, toluene, and formaldehyde in lower concentration than cigarette smoke [2,3]. Recently, E-cig users are also using this device to deliver some other harmful substances, such as tetrahydrocannabinol, cannabidiol, and hash oil [4,5]. There are several pulmonary illnesses that have been reported in users who used nicotine or cannabis extract in E-cig [4,6]. E-cig users develop changes in lung function which manifest peripheral obstructive airway involvement [7]. E-cig use can cause oxidative stress and endothelial cell dysfunction [8,9], and compromised innate immune response in lungs [10]. E-cig users have increased levels of biomarkers of inflammation and oxidative stress, reduced pro-resolving anti-inflammatory mediators, and endothelial dysfunction [11]. Cigarette smoke exposure causes various diseases involving changes in miRNA expressions in the pathogenesis of atherosclerosis, chronic obstructive pulmonary disease (COPD) and lung cancer [12][13][14][15]. Waterpipe tobacco smoking also causes various acute and chronic health effects including cardiovascular disease, chronic bronchitis and cancer [16,17]. Similarly, other study has shown that waterpipe smokers may develop COPD after chronic exposure. Exosomes are extracellular vesicles membrane bound particles of 30-150 nm. These are lipid vesicles carrying various types of biological materials such as proteins, RNAs, DNAs, transcription factors, receptors, and lipids. microRNAs are non-coding molecules of 19-23 nucleotides (nt), play major role in genes regulation associated with various biological pathways [18][19][20]. MicroRNAs recognized as molecular markers for different diseases [21,22]. microRNAs are present in higher concentration in whole RNAs [23]. Exosomes are also present in various body fluids and has several pathophysiological roles [24]. The exosomes act as signaling messengers by delivering various molecules to the adjacent cells and also via systemic circulation to distant cell population.
Cigarette smoke induced release of exosomes may be involved in the progression of chronic diseases [25,26]. Also, the role of exosomes release by cigarette smoke extract have been implicated in altering cellular functions such as endothelial dysfunction and inflammatory conditions [26][27][28][29][30]. Transfer RNAs (transfer ribonucleic acid, tRNAs) are small molecule have multiple role in cellular functions including protein synthesis. The tRNAs play a significant role in various cellular functions such as cell signaling, proliferation, differentiation, apoptosis and metabolism, as well as in regulation of gene expression at transcription and translation levels. Additionally, they have been reported to involve in DNA damage response, viral infection, neurodegeneration and cancer [31][32][33][34]. PIWI-interacting RNA (piRNA) are 24-31 nt long and form RNA-induced silencing complex with PIWI family proteins and plays a role in stem cell division, apoptosis, epigenetic control of transposons, telomeres and translational control. The piRNA serve as gene expression regulators by inducing histone modification and DNA methylation. Given the role of tRNAs and piRNA in several important biological functions and diseases, it was thought to determine whether their levels are altered in plasma exosomes from E-cig users, waterpipe and dual smokers. To the best of our knowledge this is the first study to report the detailed and differential analysis of tRNA and piRNA in plasma exosomes from E-cig users, waterpipe, and dual smokers.
The current study is performed to identify, characterize and compare the plasma-derived exosomal microRNAs in Ecig users (E-Cig), waterpipe smokers (WPS), dual smokers (DS) and cigarette smokers (CS) with normal/non-smokers (NS). Understanding the changes in plasma-derived exosomes and their microRNA profile will help in developing biomarkers during disease progression in smokers. This study will also be useful in unraveling the molecular mechanisms of disease progression in smokers. Therefore, we did comprehensive analysis of exosomal microRNA expression profile and their biological functions in plasma samples from normal/non-smokers, E-cig users, cigarette smokers, waterpipe, and dual smokers groups.

Ethic statements
All protocols, procedures and subject recruitments were approved by the Institutional Review Board (IRB)/Research Subject Review Board (RSRB) committee of the University of Rochester Medical Center, Rochester, NY [11,35]. All subjects provided written informed consent before participating in the study. The samples from nonsmokers, E-cig users, waterpipe, cigarette smokers and dual smokers (cigarette and waterpipe smokers) used in this study are described previously [11,35].

Biochemical and molecular methods
The following methods were used by us in our recent publication and some of the methods used here are given in detail [29,30].

Collection of blood plasma
The venous blood (20-25 ml) was collected from cigarette smokers, waterpipe smokers, smokers (cigarette and waterpipe), Non-smokers and E-cig users.

Plasma exosome purification
Plasma was used for exosomal isolation as per the manufacturer's instructions (Norgen Biotek Corporation, Cat# 57400).

Transmission electron microscopy
Transmission electron microscopy (TEM) was used to visualize exosomes as reported earlier [30].
Nanoparticle tracking analysis (NanoSight NS300) for size and particle concentrations NanoSight Technology NS300 was used for particle size analyses and concentration of plasma derived EVs by nanoparticle tracking analysis (NTA) as described earlier [30,36].

Small RNA library construction and sequencing
Exosomal RNA isolated was used for small RNA library using the small RNA library preparation kit (Small RNA library prep kit, cat# 63600, Norgen Biotek Corp., Canada) as reported previously [37], and quantified using High Sensitivity DNA Analysis kit on the Agilent 2100 Bioanalyzer system (Agilent Technologies, USA). Libraries were sequenced on the Illumina NextSeq 500 sequencing platform using the NextSeq 500/550 High Output kit v2 at the Norgen Biotek Corp.

Sequence read mapping and RNA annotation
The sequence data (Illumina NextSeq 500) derived fastq format was used in the Genboree Workbench's exceRpt small RNA-seq pipeline for mapping with human genome [38].

Data and statistical analysis
The RNA-Seq fastq files were analyzed [29], and the trimmed mean of M values (TMM) method was used to normalize the summarized count data for further miRNA analysis [42]. The DESeq package of R/Bioconductor was used for miRNA differential analysis across different groups. Pairwise miRNA differences were conducted through linear contrasts within the DESeq2 model framework [43]. The false discovery rate (FDR) was accounted for at 5% using the Benjamini-Hochberg method for each pairwise comparison. A heatmap of significant miRNA selected from the pairwise comparisons was generated using the pheatmap function with the ward. D2 clustering method in statistical analysis software R [29]. A receiver operating characteristic (ROC) curve was generated using the ROCit package using software R (R Core Team, 2019). Pathway and gene enrichment analyses were conducted through the FunRich software [44].
Functional over-representation analysis and gene set enrichment analysis of microRNA Enrichment analysis and annotation tool, such as miEAA (http://ccb-compute2.cs.uni-saarland-de/mieaa_tool/; accessed in August 2019) was used to identify various microRNA expression levels and functions [45]. Gene-Trail was used for miRNA set enrichment (FDR adjustment) analyses and gene ontology, pathways, and disease and over-representation analysis [29].

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. Immunoblot analysis showed the presence of exosomal markers CD63 and CD81 and negative for calnexin and any endoplasmic reticulum contamination (Fig. 1a). 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).

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 (14,031,359 ± 2,019,014), CS (10,688,617 ± 917,325), WP (8,765,589 ± 768,710), E-cig users (9,570,944 ± 455,585) and DS (8,988,900 ± 798, 445). 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 and 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 Representative TEM images of plasma-derived EVs/Exosomes. c Particle size depicted as mean and mode, and particle concentration were estimated using NanoSight NS300 (n = 3/ group). Equal amount of protein was loaded as assessed by Ponceau staining. Figure 1a full length western blot images are available in Suppl. Fig. 1A and labeled as "original western blots". Positive control: human lung homogenate Fig. 2 The small RNA analysis quality control results of individual samples in non-smokers, cigarette smokers, waterpipe smokers, E-cig users and dual smokers 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. 2 A-D) are presented as volcano plot. The volcano plots were made by plotting the -log 10 of adjusted pvalues on the y-axis, and the log 2 fold change between two groups on the x-axis showing the up-and down regulation appearance similar distance from the center. Some micro-RNAs 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. 3 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 dual smokers versus non-smokers
The DS vs NS groups have shown total 20 differential expressed microRNAs that are shown in Suppl. 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 (Suppl. Table 5). Among the top fold changed microRNAs, upregulated in E-cig are
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. 9a-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  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 (tRNA Val , tRNA Glu , tRNA Asp , tRNA Gly , tRNA Arg and tRNA His ) and decrease in tRNA Cys . In addition, NS vs. DS group also showed significant increase in tRNA Ile (Suppl. Tables 10, 11, 12, and 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. 5).
The differentially expressed tRNAs from NS vs. CS, NS vs. WP, NS vs. E-cig users, and NS vs. DS (Suppl. Fig. 6 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 (Suppl. Fig. 7). 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 Suppl. Fig. 8. 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 (tRNA Val , tRNA Glu , tRNA Asp , tRNA Gly , tRNA Arg and

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 (Suppl. Fig. 9). 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 (Suppl. 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 (Suppl. Fig. 10E-H). The color scale at the right side in indicates the relative expression level of piRNA in all samples. Venn diagram showing piRNA for different groups is given in Suppl. Fig.  11). 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 are given in Suppl. Table 14.
Overall sequencing and metadata The above described sequence data in fastq format with associated metadata generated from current study as described above have been deposited into the NCBI Sequence Read Archive (SRA) public repository. The deposited BioProject database with associated URL links (See Availability of Data and Materials) are described in Suppl. Tables 15, 16, 17.

Discussion
Exosomes are secreted into biological fluids, and communicate with cells by transferring exosomal content including microRNAs to the recipient cells, thereby mediate several cellular and biological processes [37,46,47]. In addition, exosome releases micoRNAs during airways lung injury that involve in the progression of lung diseases [29,47].
In this study, we have investigated the small RNA profile of plasma samples of NS, E-cig users, WP and DS using high-throughput sequencing. We have presented a comprehensive microRNAs profile of plasma exosomes in all the four groups in comparison to the NS. Our results showed that some microRNAs are common in all tobacco/E-cig users. Interestingly, we also found group specific microRNAs in different groups.
To perform a comprehensive analysis of small RNA species, we analyzed our data for several key abundant RNA biotypes. The input reads were lower in all groups, but the significant lower were seen in WP and DS. Several biotypes of known biological importance were identified, such as piRNA, tRNA, snoRNA, snRNA, lincRNA, Mt_tRNA, Mt_rRNA and microRNAs. The microRNA was most diverse RNA biotype found in comparison to all biotypes in samples from all groups. Biotypes Mt_tRNA and snoRNA were significantly lower in all groups in comparison to NS, suggesting nicotinespecific effect. Recent studies of small RNA sequence analysis of exosomes in various biofluids from healthy donors, smokers, and COPD patients reported showing changes in relative distribution of biotypes, such as microRNA, tRNA, Mt_tRNA, and snoRNA [29,37]. We have performed further analysis of tRNA and piRNA in all NS, CS, E-cig, WP and DS groups, showing a differential regulation of these transcripts.
There are several studies related to expression of coding and non-coding RNA in CS exposed cellular models and COPD patients [48][49][50]. The studies related to transcriptome profiling of gene expression in E-cig vapor and CS exposed human bronchial epithelial cells and nasal epithelial cells from E-cig users [29,51,52]. Recently, microRNA expression profiling of E-cig exposed human airway epithelial cells has been reported [53]. The waterpipe smokers have differential gene expression in small airway epithelial cells and alveolar macrophages in comparison to non-smokers [54]. Further, waterpipe smokers small airway epithelial cells have epigenetic related changes in gene expression [55]. However, the studies related to plasma exosomal microRNA in E-cig users and WP are still unknown.
To our surprise, when we analyzed the overlap of microRNAs in all four groups, there were 6 downregulated miRNAs (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 the groups. The microRNAs hsa-let-7a-5p and hsa-let-7f-5p are involved in NSCLC and typical and atypical carcinoid tumors of the lung, respectively [66,67]. A ROC curve evaluation of diagnostic utility using miRNA hsa-let-7a-5p to distinguish normal subjects from tobacco users showed high sensitivity and specificity (greater than 80%). In addition, miRNA hsalet-7a-5p was also found to be an identifier for the cigarette smoke treatments in lung cells [29]. This supports the validity and the usage of this miRNA as a possible theranostic marker in the future. Further, exosomal microRNA-21 derived from bronchial epithelial cells play a role in myofibroblast differentiation in COPD by cigarette smoking [68]. Further, when groups compared for upregulated microRNAs, the only micro-RNA has-miR-29b-3p was common in all groups (Fig. 7 c) which is involved in epithelial-mesenchymal transition [57]. As changes in these microRNAs expression are common in all the groups and may reflect tobaccospecific effect. There are four downregulated microRNAs (hsa-miR-224-5p, hsa-miR-423-3p, hsa-miR-365a-3p|hsa-miR-365b-3p, and hsa-miR-193b-3p), and one upregulated hsa-miR-30e-5p expressed specifically to Ecig group (Fig. 7 b, c). These microRNAs may be specific biomarkers for E-cig users.
Next, we performed FunRich functional enrichment analysis of differentially expressed microRNAs to explore the potential target genes in NS vs. E-cig, NS vs. CS, NS vs. WP and NS vs DS pairwise comparisons. 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, 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 of regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolism with lowest p value were 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 was 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. Our analysis revealed that many vital biological functions associated genes are potential target of selective plasma exosomal microRNAs in all four groups, which may be associated with lung pathophysiological conditions.
To determine the target genes of differentially expressed microRNAs in different groups, we identified 3244 (NS vs.  Table 9. The large numbers of genes are the target of differentially expressed microRNAs in all the groups suggesting several compromised biological functions in these groups. Studies have shown that tRNA and piRNA are associated with several human diseases. Our differential expression 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 8 tRNAs. All four groups have significant increase in 6 tRNAs (tRNA Val , tRNA Glu , tRNA Asp , tRNA Gly , tRNA Arg and tRNA His ) and decrease in tRNA Cys . In addition, NS vs. DS group also showed significant increase in tRNA Ile (Suppl. Tables 10,11,12,13). There is no information available regarding the tRNA expression in plasma exosomes of E-cig users, WP, and DS. A recent study from our lab has shown that five tRNAs differentially expressed (tRNA Lys , tRNA Gly , tRNA Tyr , tRNA Glu and tRNA) in nonsmokers vs COPD and five tRNAs (tRNA Gly , tRNA Tyr , tRNA Leu , tRNA and tRNA Met ) in smoker vs. COPD pairwise comparison among the 25 enriched tRNAs [29]. tRNAs plays important role in tumor development and associated with the pathological characteristics of lung adenocarcinoma and cancer-specific survival [69]. The significant differential expression of piRNA in plasma exosomes of NS vs. CS (piR-004153, pir-019825), NS vs. WP (piR-004153, piR-019825, piR-000552, piR-014620, and piR-020450), NS vs. E-cig (piR-016658, piR-016659, piR-019825, piR-000552, and piR-017591), and NS vs. DS (piR-020365, piR-000552, and piR-017591) by pairwise comparisons were found. Further, pairwise comparison of piRNA expression in CS vs. WP (piR-000552 and piR-020450), CS vs. E-cig (piR-000552), CS vs. DS (piR-000552), and WP vs. DS (piR-019825, piR-014620, and piR-020450) also showed significant changes (Suppl. Table 14). piRNA form RNA-induced silencing complex with PIWI family proteins and plays a role in stem cell division, apoptosis, epigenetic control of transposons, telomeres and translational control. There is no study as such to study the differential expression of piR-NAs in plasma exosomes of WP, E-cig users, and DS. A differential expression of piRNA-004153 was observed in non-smokers vs. smokers, similar to our finding in NS vs. CS [29].
We mapped input reads to human genome and classified to the various small RNA biotypes. To our surprise, significantly lower counts of Mt-tRNA and small nucleolar RNA were observed in CS, WP, E-cig and DS (P < 0.05) in comparison to NS (Fig. 3 a, b). Small nucleolar RNAs are non-coding RNAs consist of 60-300 nt long which accumulate mostly in nucleoli. snoRNAs are involved in various pathophysiological processes. The aberrant expression of snoRNA may induce cell transformation, tumorigenesis, and metastasis. Mt-tRNA plays a role in mitochondrial protein synthesis and oxidative phosphorylation (OXPHOS) to produce energy (ATP) and reactive oxygen species (ROS). Any alteration in the level of Mt-tRNA may have impact on mitochondrial function. The mitochondrial dysfunction has been associated with several lung diseases, such as asbestos-related lung fibrosis and COPD [70,71]. The Mt-tRNA gene mutations were reported in idiopathic pulmonary fibrosis, lung cancer and hypertension [72][73][74]. However, we did not perform detailed analysis of Mt-tNRA and snoRNA genes in our samples. In our future studies, we are interested to perform a detailed analysis of Mt-tRNA, snoRNA, snRNA and other differentially expressed small RNA enriched biotypes.

Conclusions
In conclusion, we have performed a comprehensive plasma exosomes small RNA-sequence analysis that includes exosome isolation, purification, RNA extraction, library preparation, RNA sequencing, and RNA annotation in samples from NS, CS, WP, E-cig users and DS. The isolated RNAs from plasma exosomes were used for small RNA-sequencing analysis. Our data show the enrichment of various non-coding RNAs that include microRNAs, tRNAs, piRNAs, snoRNA, snRNAs, Mt-tRNAs, and other biotypes. Further, the detailed differential expression analysis of microRNAs, tRNAs and piRNA showed significant changes between pairwise comparisons of different groups. Gene set enrichment analysis showed significant changes in the top six enriched functions that were biological pathway, biological process, molecular function, cellular component, site of expression and transcription factor in all groups. The pairwise comparison of tRNAs and piRNA also revealed significant changes in differential expression in all groups. We mapped input reads to human genome and classified to the various small RNA biotypes that showed significant lower counts of Mt-tRNA and snoRNA in all groups. Our comprehensive differential transcriptome analysis done in this study will help to better understand molecular mechanisms of plasma exosome non-coding RNAs, and in developing biomarkers that may be helpful in diagnosis and therapy (theranostics) of lung injury and disease mechanisms in smokers and vapers.