Exposure to Fine Particulate Matter Air Pollution Alters mRNA and miRNA Expression in Bone Marrow-Derived Endothelial Progenitor Cells from Mice

Exposure to fine particulate matter (PM2.5) air pollution is associated with quantitative deficits of circulating endothelial progenitor cells (EPCs) in humans. Related exposures of mice to concentrated ambient PM2.5 (CAP) likewise reduces levels of circulating EPCs and induces defects in their proliferation and angiogenic potential as well. These changes in EPC number or function are predictive of larger cardiovascular dysfunction. To identify global, PM2.5-dependent mRNA and miRNA expression changes that may contribute to these defects, we performed a transcriptomic analysis of cells isolated from exposed mice. Compared with control samples, we identified 122 upregulated genes and 44 downregulated genes in EPCs derived from CAP-exposed animals. Functions most impacted by these gene expression changes included regulation of cell movement, cell and tissue development, and cellular assembly and organization. With respect to miRNA changes, we found that 55 were upregulated while 53 were downregulated in EPCs from CAP-exposed mice. The top functions impacted by these miRNA changes included cell movement, cell death and survival, cellular development, and cell growth and proliferation. A subset of these mRNA and miRNA changes were confirmed by qRT-PCR, including some reciprocal relationships. These results suggest that PM2.5-induced changes in gene expression may contribute to EPC dysfunction and that such changes may contribute to the adverse cardiovascular outcomes of air pollution exposure.


Introduction
Analysis of abundant epidemiological data has identified robust associations between exposure to fine particulate matter (PM 2.5 ) air pollution and adverse health outcomes in general [1], and cardiovascular disease in particular [2,3]. These studies have linked acute exposure to PM 2.5 with an increased incidence of myocardial infarction, stroke, cardiac arrhythmias, and sudden cardiac death [3][4][5], while chronic exposures appear to accelerate the progression of atherosclerotic disease [6,7]. Physiologically, these adverse cardiovascular outcomes are likely a consequence of a damaged or dysfunctional endothelium. Not only does the endothelial layer provide a physical separation between blood and Male C57BL/6J mice were obtained at 11 weeks of age from Jackson Laboratories and acclimated for 1 week in pathogen-free housing conditions. The animals were kept on a 12 h light/dark schedule and were supplied with food and water ad libitum. Mice (n = 3; 1 cage) were randomly selected and exposed to CAP generated by a Versatile Aerosol Concentration and Enrichment System (VACES) in the Inhalation Facility of the University of Louisville, as previously described [16,17]. Mice exposed to HEPA-filtered air (n = 3; 1 cage) were used as controls. Exposures to CAP or filtered air were 6 h per day (7 a.m. to 1 p.m.) for 30 consecutive days, and over this time course, the average CAP level was 67 µg/m 3 . These animals were monitored daily by veterinarian staff for well-being. The animals were euthanized (pentobarbitol injection (50 mg/kg, i.p., and exsanguination)) immediately after the final exposure (CAP-and air-exposed animals intermingled) and bone marrow cells isolated and EPCs cultured as previously described [18]. After 10 days in culture, the cells were harvested, and total RNA was isolated using the miRNeasy isolation kit (Qiagen). RNA quality was confirmed using an Agilent Bioanalyzer. A schematic of the experimental protocol illustrated in Figure 1.

Library Preparation
Libraries were prepared using the TruSeq Stranded mRNA Library Prep kit (20020594, Illumina, San Diego, CA, USA). In brief, poly-A enrichment of RNA samples (250 ng) was accomplished in two purification steps using oligo dT beads. Purified mRNA was chemically fragmented through a heated digestion with divalent cations. Eluted samples were then used for first-and second-strand cDNA synthesis, and the resulting cDNA was isolated using Agencourt AMPure XP Beads (A63881, Beckman Coulter Life Sciences, Indianapolis, IN, USA). After adenylation of the 3′ end to prevent blunt end fragments from ligating to each other, dual-index adapters were ligated, and the samples were again purified with the same beads. Adapter ligated DNA fragments were then enriched and amplified with PCR to create the final cDNA libraries. The concentrations of these libraries were measured by Qubit dsDNA HS Assay Kit (Q32851, Invitrogen, Waltham, MA, USA), and their quality was assessed using the Agilent Bioanalyzer.

Sequencing
Similar molar amounts were pooled, and their quantity and quality were assessed using the MiSeq Reagent Nano Kit V2 (MS-103-1001, Illumina). The libraries and PhIX control (FC-110-3001, Illumina) were denatured and diluted using the manufacturer's directions and sequenced on an Illumina MiSeq. On the basis of the sequencing behavior of the libraries compared to PhIX, we re-pooled, denatured, and diluted equal amounts of the libraries, and 1.8 pM was loaded on the NextSeq with 1% PhIX spike-in. After the first sequencing run, the libraries were re-pooled for a second run with the goal of having an equal number of reads per sample after results from both runs were combined, which gave an average of 35 million reads per sample. All sequencing was performed at the University of Louisville Genomics Facility on the Illumina NextSeq 500 using the NextSeq 500/550 75 cycle High Output Kit v2 (FC-404-2002, Illumina).

Data Analysis
The 48 fastq single-end raw sequencing files representing the air and CAP conditions with three biological replicates and eight lances per replicate were downloaded from Illumina's BaseSpace (https://basespace.illumina.com/, accessed on 13 December 2019). The eight single-end raw fastq files across 8 sequencing lanes for each replicate designated L001 through L008 in the file names from the sequencing run were concatenated into one single-end fastq file using the unix cat command. Quality control (QC) of the raw sequence

Library Preparation
Libraries were prepared using the TruSeq Stranded mRNA Library Prep kit (20020594, Illumina, San Diego, CA, USA). In brief, poly-A enrichment of RNA samples (250 ng) was accomplished in two purification steps using oligo dT beads. Purified mRNA was chemically fragmented through a heated digestion with divalent cations. Eluted samples were then used for first-and second-strand cDNA synthesis, and the resulting cDNA was isolated using Agencourt AMPure XP Beads (A63881, Beckman Coulter Life Sciences, Indianapolis, IN, USA). After adenylation of the 3 end to prevent blunt end fragments from ligating to each other, dual-index adapters were ligated, and the samples were again purified with the same beads. Adapter ligated DNA fragments were then enriched and amplified with PCR to create the final cDNA libraries. The concentrations of these libraries were measured by Qubit dsDNA HS Assay Kit (Q32851, Invitrogen, Waltham, MA, USA), and their quality was assessed using the Agilent Bioanalyzer.

Sequencing
Similar molar amounts were pooled, and their quantity and quality were assessed using the MiSeq Reagent Nano Kit V2 (MS-103-1001, Illumina). The libraries and PhIX control (FC-110-3001, Illumina) were denatured and diluted using the manufacturer's directions and sequenced on an Illumina MiSeq. On the basis of the sequencing behavior of the libraries compared to PhIX, we re-pooled, denatured, and diluted equal amounts of the libraries, and 1.8 pM was loaded on the NextSeq with 1% PhIX spike-in. After the first sequencing run, the libraries were re-pooled for a second run with the goal of having an equal number of reads per sample after results from both runs were combined, which gave an average of 35 million reads per sample. All sequencing was performed at the University of Louisville Genomics Facility on the Illumina NextSeq 500 using the NextSeq 500/550 75 cycle High Output Kit v2 (FC-404-2002, Illumina).

Data Analysis
The 48 fastq single-end raw sequencing files representing the air and CAP conditions with three biological replicates and eight lances per replicate were downloaded from Illumina's BaseSpace (https://basespace.illumina.com/, accessed on 13 December 2019). The eight single-end raw fastq files across 8 sequencing lanes for each replicate designated L001 through L008 in the file names from the sequencing run were concatenated into one single-end fastq file using the unix cat command. Quality control (QC) of the raw sequence data was performed using FastQC (version 0.10.1). The FastQC results showed the sequences for all the samples had a Phred score above 30, which indicated 99.9% accuracy in base calling and were considered to be good quality. Therefore, sequence trimming was not performed. The concatenated sequence samples were directly aligned to the Mus musculus reference genome assembly (mm10.fa) using STAR (version 2.6) [23]. The alignment rate ranged from 99.12 to 99.98 percent across the samples. For the identification of differentially expressed genes (DEG) using DESeq2 [24], the raw reads counts were extracted from the STAR-aligned bam format files using HTSeq version 0.10.0 and the Mus musculus gtf file (mm10.gtf). Supplementary Table S1 indicates the number of raw reads and number of reads successfully aligned for each of the six samples. The raw read counts for six samples being normalized using the relative log expression (RLE) method and then filtered to exclude genes with fewer than 10 counts across the samples with the aid of DESeq2. DEGs (CAP versus air) were identified and a Benjemini-Hochberg adjusted p-value (q-value) cutoff ≤0.05 with log 2 fold change (FC) ≥0.6 was used to define differential expression. RNAseq data are available at GEO accession number GSE153038.

. Library Preparation
Libraries were prepared with 100 ng from each of the RNA samples using the QIAseq miRNA Library Kit (331502, Qiagen, Beverly, MA, USA) following the manufacturer's instructions. Briefly, specially designed 3 and 5 adapters were ligated to mature miRNAs. The ligated miRNAs were then reverse-transcribed to cDNA and purified with magnetic beads (QIAseq miRNA NGS beads). The eluted cDNA was amplified with a universal forward primer and a sample indexing reverse primer. The amplified libraries were again purified with magnetic beads to remove adapter dimers and unwanted RNA species. Library quantity and quality were analyzed on an Agilent Bioanalyzer using the Agilent high sensitivity DNA Kit (5067-4626, Agilent Technologies, Santa Clara, CA, USA).

Sequencing
Similar molar amounts of the libraries were pooled and run on MiSeq to test for quantity and quality. Then, the libraries and PhIX control (FC-110-3001, Illumina) were denatured and diluted following manufacturer's directions, and 300 µL of each were combined and sequenced on Illumina MiSeq. The concentration of the libraries was corrected and the libraries re-pooled on the basis of the results from the initial MiSeq run. After denaturation and dilution, sequencing was performed on the University of Louisville Genomics Facility's Illumina NextSeq 500 using the NextSeq 500/550 75 cycle High Output Kit v2 (FC-404-2002, Illumina).

Data Analysis
The 48 fastq single-end raw sequencing files representing the air and CAP conditions with three biological replicates and eight lances per replicate were downloaded from Illumina's BaseSpace (https://basespace.illumina.com/, accessed on 3 January 2020). The eight single-end raw sequencing files (fastq) for each replicate designated L001 through L008 from the sequencing run were concatenated into one single-end fastq file using the unix cat command. Quality control (QC) of the raw sequence data was performed using FastQC (version 0.10.1). The FastQC results indicated quality trimming was unnecessary since the minimum quality value of a Phred score for all samples was well above 30 (1 per 1000 error rate). Given this, preliminary adapter trimming was performed on each of the samples to remove the Qiagen 3 and 5 adapter sequence with Trimmomatric v0.33 [25]. For all of the samples, a unique peak around 22 bp was identified that indicated that little if any other non-coding RNAs were selected in the library processing. As expected, this peak represented the length of the mature miRNAs. The data were further analyzed using mirDeep2 [26] and DESeq2 [24]. Briefly, the trimmed sequence reads were directly mapped to the mouse mm10 reference genome assembly using mirDeep2 package (v0.07) implementing bowtie mapper (v1.1.1) [27]. Supplementary Table S2 indicates the number of raw reads, number of reads after trimming, and number of reads successfully aligned for each of the six samples. The aligned sequences were then used as inputs into mirDeep2 along with the latest release of the mirBase version 22 database (v22) containing mature miRNA and miRNA hairpin sequences. The result was a file containing the number of reads mapping to each of the 2061 mouse miRNAs. After quantification, the resulting counts for each miRNA in each sample were combined into a reads matrix. Using the counts table resulting from the previous step, we determined the differentially expressed miRNAs (q ≤ 0.05; upregulated: log 2 FC ≥ 0.6; downregulated: log 2 FC < −0.6) using DESeq2. miRNA-seq data are available at GEO accession GSE153038.

In Silico INGENUITY Network Analysis
Canonical pathway and biological function analysis of DEGs and differentially expressed miRNAs (DE miRNAs) was performed using Ingenuity Pathway Analysis (202 Qiagen version 57662101).

qRT-PCR
Quantitative RT-PCR (qRT-PCR) of selected mRNA species was performed by initial, first-strand cDNA synthesis (High Capacity cDNA Reverse Transcription Kit, Thermo-Fisher), followed by real-time PCR using specific primers as recommended by the manufacturer (Supplementary Table S3: TaqMan gene expression assays, Thermo-Fisher). Quantitative RT-PCR of selected miRNA species was performed using specific TaqMan MicroRNA assays (Supplementary Table S3; Thermo-Fisher). All amplifications were performed on a Quant Studio 5 thermal cycler (Applied Biosystems). The internal controls used were Hprt for mRNAs and sno202 for miRNAs. Fold changes were calculated using the Ct method [28] and were reported with a standard error.

mRNA Analysis
To evaluate similarities and differences in mRNA expression between the air and CAP groups, we initially performed a principal component analysis (PCA) on the relative log expression (RLE)-normalized counts of the six samples for all genes. In this analysis (Figure 2A), the two principal components (PC1 and PC2) represent two-dimensional data with 71.7% and 11.9% variation in the six mRNA samples, respectively. The PC1 axis (the first principal direction with all six mRNA samples) showed the largest variation, which explained most of the variance in the original data. The PC2 axis is the second most important direction and orthogonal to the PC1 axis. For these samples, the CAP group (blue) was partially separated from the air group (red), and some variation within the groups was observed.
To identify transcriptional alterations driven by CAP exposure, we performed mRNA expression profiling. Hierarchical clustering of the expression of 68 differentially expressed genes (DEGs) revealed a distinct transcriptomic profiling pattern between the samples as illustrated in a heatmap ( Figure 2B) from low expression (green) to high expression (red). Given a q-value cutoff ≤ 0.05, we identified 166 genes that were differentially expressed. This included 122 genes that were upregulated and 44 that were downregulated in the CAP vs. air groups. The most upregulated or downregulated DEGs are located in a volcano plot (p-value vs. FC) in Figure 2C. We further analyzed these DEGs according to a fold change. Given a q-value ≤ 0.05, we identified 104 that were upregulated by >1.5 fold (log 2 FC ≥ 0.6). The top 20 of these genes are listed in Table 1, while the complete list is provided in Supplementary Table S4. Given a q-value ≤ 0.05, we identified 34 that were downregulated by <−1.5 fold (log 2 FC ≤ −0.6). The top 20 of these genes are listed in Table 2, while the complete list is provided in Supplementary Table S5.

miRNA Analysis
Using the same samples, we performed a similar analysis of miRNA expression. PCA analysis of these samples ( Figure 3A) revealed 69.0% and 19.2% variation of the six miRNA samples. This analysis reveals clear separations between the CAP (blue) and air (red) samples. Hierarchical clustering of the expression of 38 differentially expressed miRNAs (DE miRNAs) revealed a unique miRNA profiling pattern, as illustrated in a heatmap ( Figure 3B). Given a q-value cutoff ≤ 0.05, we identified 108 unique miRNAs that were differentially expressed. This included 55 miRNAs that were upregulated and 53 that were downregulated in the CAP vs. air groups. The most upregulated and downregulated miRNAs are located in a volcano plot ( Figure 3C). We further analyzed these DE miRNAs according to fold change. Given a q-value ≤ 0.05, we identified 19 that were upregulated by >1.5 fold (log 2 FC ≥ 0.6) and 19 that were downregulated by >1.5 fold (log 2 FC ≤ −0.6). These miRNAs are listed in Table 3.
heatmap ( Figure 3B). Given a q-value cutoff ≤ 0.05, we identified 108 unique miRNAs that were differentially expressed. This included 55 miRNAs that were upregulated and 53 that were downregulated in the CAP vs. air groups. The most upregulated and downregulated miRNAs are located in a volcano plot ( Figure 3C). We further analyzed these DE miRNAs according to fold change. Given a q-value ≤ 0.05, we identified 19 that were upregulated by >1.5 fold (log2 FC ≥ 0.6) and 19 that were downregulated by >1.5 fold (log2FC ≤ -0.6). These miRNAs are listed in Table 3.

Biological Function Analysis
To identify biological functions involving the DEGs and DE miRNAs, we used ingenuity pathway analysis (IPA). For the 166 DEGs, we found the most impacted, biological functions involved cellular movement, tissue development, cardiovascular system development and function, cellular assembly and organization, and cellular function and maintenance (Table 4). Other functions implicated to a lesser extent were those involving metabolic processes and signal transduction. Similarly, many of the 108 DE miRNAs were predicted to be involved in functions involved in cell movement, survival, development, growth, cardiovascular system development and function, signaling, and organization ( Table 5). Listed are the top biological functions impacted by the mRNA changes we detected, their significance, the number of DEGs that impact these functions, and representative DEGs as listed in Supplementary Tables S4 and S5. Selected genes are boxed in Figure 2B. Listed are the top 10 biological functions impacted by the miRNA changes we detected, their significance, the number of DE miRNAs that impact these functions, and representative miRNAs as listed in Table 3. Selected miRNAs are boxed in Figure 3B.

qRT-PCR Validation
To validate results from our transcriptomics analysis, we performed qRT-PCR for a subset of differentially expressed mRNAs ( Figure 4, Table 6) that are involved in those functions identified in the IPA analysis (Table 4). We also performed a similar qRT-PCR analysis of selected DE miRNAs (Figure 4, Table 6). Our results both confirmed the direction of change (upregulated vs. downregulated) and had a high level of consistency with regards to the extent of the change (fold change). Since there is a reciprocal relationship between miRNA expression and the expression of their target mRNAs, we sought to confirm some of these relationships among our DEGs and DE-miRNAs. Consistent with prior reports [29][30][31][32], our RNA-seq data and our qRT-PCR results confirmed reciprocal relationships between mir-214-3p and Ccl5, mir-450a-5p and Dusp-10, and mir-92a-3p and Tgfb2 (Figure 4, Table 7). Finally, we used a bioinformatics database [30] to identify additional, documented, reciprocal relationships between the 166 DEGs and the 108 DE miRNAs ( Figure 5). Our detected changes in the levels of mir-7b-5p, mir-466k, mir-499-5p, mir-11-3p, mir-3110-5p, mir-34b-5p, and mir-450a-5p (Table 3) were consistent with the inverse expression levels of some target mRNAs, further supporting the validity of our data.

Discussion
Exposure to PM 2.5 is associated with cardiovascular morbidity and mortality in humans [2,3] and promotes pre-clinical vascular disease in animal models of exposure [33][34][35]. Some studies propose that quantitative [16,17,36,37] or qualitative [18] defects in EPCs may underlie these outcomes. Yet the mechanistic basis for PM 2.5 -mediated EPC depletion or dysfunction is not completely understood. To this end, in this study, we used a transcriptomic approach to identify gene expression changes in EPCs isolated from the bone marrow of mice exposed to CAP. Compared with control cells isolated from mice breathing filtered air, we identified multiple mRNAs and miRNAs that were either upregulated or downregulated in exposed mice. The top predicted biological functions impacted by these changes are those involving cell movement, cell and tissue development, cell death and survival, and cellular assembly and organization. These predictions are consistent with the role of EPCs in vascular maintenance and repair, and with previously identified, CAP-induced impairments in proliferation, in vitro tube formation, and angiogenesis [18]. Thus, these PM 2.5 -induced changes in EPC gene expression may form the basis of EPC defects and be early indices of incipient cardiovascular disease.
Multiple mechanisms may contribute to the gene expression changes we observed. Epigenetic modifications are implicated, given that the RNA samples used for transcriptomics analysis were isolated from EPCs after 10 days of culture following bone marrow collection from exposed mice. One possible epigenetic mechanism is the post-transcriptional regulation of mRNA levels through targeting miRNAs. Consistent with this idea, we did identify reciprocal relationships between representative dysregulated genes and at least one of their documented, targeting miRNAs (Table 7). Supporting this mechanism of PM 2.5 -induced gene expression changes, other studies also suggested that particle inhalation [38,39] or exposure to particles in cell culture models [40,41] induces the upregulation or downregulation of specific miRNAs. Other epigenetic mechanisms may involve the modification of nuclear material. Methylation of cytosine residues in DNA promoter regions generally inhibits gene expression [42]. In support of this mechanism, some prior studies have determined that exposure to PM 2.5 induces site-specific DNA methylation, thereby limiting mRNA and protein expression [43][44][45]. DNA methylation is catalyzed by the family of DNA methyltransferases (DNMTs) and, consistently, prior evidence also suggested that PM 2.5 exposure can activate DNMTs [46]. Furthermore, it has also been shown that reactive oxygen species (ROS) can upregulate DNMT expression [47]. As ROS production is an early and robust outcome of PM 2.5 inhalation [48], such exposures can thus effectively regulate DNMT expression or activity in promoting methylation. The regulation of gene expression through promoter DNA methylation may be a common outcome of exposure to other environmental toxins as well, and has also been reported following exposure to cigarette smoke [49], volatiles [50], or metals [51], which can also be a constituent of air pollution particles. Finally, in addition to DNA modifications, histone modifications (e.g., methylation, acetylation) can alter chromatin structure, thereby facilitating or impairing transcriptional processes [42]. Recent evidence likewise suggests PM 2.5 may promote histone acetylation and that this may be accomplished in a dose-and time-dependent manner [52,53].
Exposure to PM 2.5 is known to negatively impact vascular function and/or the ability of EPCs to maintain vascular homeostasis [18,33,34]. This is supported by our function analysis in the current study, which suggested the influence of PM 2.5 on genes that regulate cardiovascular system development and function. Furthermore, some of the specific gene expression changes we observed are consistent with PM 2.5 -induced impairments in EPC function. EPC-expressed Cxcr3 is a receptor for multiple CXC chemokines that collectively promote the mobilization, migration, and homing of these cells [54,55]. We found that Cxcr3 is downregulated in EPCs from CAP-exposed mice (Tables 2 and 6), thus potentially limiting the availability of these cells for the repair of distal tissue damage. Other gene expression changes we observed are consistent with an impairment in EPC differentiation. For instance, signaling through the angiopoietin receptor Tek (Tie2), a commonly used marker of mature EPCs, is important in controlling the differentiation of these cells [56]. Its downregulation (Tables 2 and 6) suggests inefficient EPC maturation in exposed animals. Among its effects, the proto-oncogene Myc appears to support cellular de-differentiation and, indeed, was one of the original factors used in the generation of induced pluripotent stem cells [57]. Its upregulation in our study (Supplementary Table S4) further suggests a mechanism whereby PM 2.5 exposure may impede EPC differentiation. Finally, Ccl5 (Rantes) is a pro-inflammatory cytokine released by early EPCs but to a much lesser extent by late EPCs [58]. The upregulation of Ccl5 in EPCs derived from CAP-exposed mice (Tables 1 and 7) is consistent with the idea that these cells remain in an early, immature state, being inefficient in promoting vascular repair. Other DEGs we identified and confirmed by qRT-PCR (Itga3, Marcks, Dock9, Dusp-10; Tables 6 and 7) play general roles in cell adhesion, cytoskeletal organization, growth, and differentiation, and may likewise play essential roles in EPC mobilization, homing, migration, and tissue repair. Some of the DEGs and DE-miRNAs identified in our analysis (Tables 1-3) have an uncertain role in EPC function, and a delineation of their roles requires more sophisticated approaches in future studies. Nevertheless, the current work suggests PM 2.5 -induced changes in gene expression may underlie widespread functional defects in EPCs that could lead to impairments in vascular repair and promotion of CVD.
Our findings are significant because they show that PM 2.5 -induced gene expression changes can occur in a tissue (bone marrow) that is distal from the tissue of initial exposure (lungs). These long-range PM 2.5 effects are believed to initiate largely from redox-catalyzed lipid peroxidation and the generation of toxic intermediates such as aldehydes in the lungs. These intermediates are then distributed systemically through peripheral blood circulation to impact the function of other tissues [59]. There is strong evidence that lipid peroxidation products and aldehydes can influence DNA methylation and histone acetylation patterns [60], as well as miRNA expression levels [61,62]. Thus, PM 2.5 -induced oxidative stress may be causative in the gene expression changes we observed. Supporting this idea, measures taken to limit oxidative stress [18,37], or to neutralize aldehydes [63], appear to limit the downstream, toxic effects of PM 2.5 exposure and can reverse some of gene expression changes we previously observed [18].
One limitation of this study is that it is not clear if our results are applicable to all exposure scenarios, as PM 2.5 concentration and composition may vary according to locale and in a seasonal and temporal manner. However, the location of our exposure facility in downtown Louisville assures that particles are derived from diverse and voluminous vehicle exhausts and are generally reflective of traffic-derived PM 2.5 in any locale. Furthermore, with a 30-day exposure period, we expect the effects of daily variations in composition would be minimized. Consistent with this idea, the predicted biological functions impacted by the gene and miRNA expression changes we observed are consistent with air pollution-induced quantitative and qualitative defects in EPCs that have been previously reported [18,36,37].

Conclusions
In summary, we identified multiple changes in the expression of distinct mRNAs and miRNAs in bone marrow-derived EPCs isolated from mice inhaling CAP. Some of the predicted biological functions impacted by these changes are those regulating EPC availability or vascular repair potential. These results shed further mechanistic insight with regards to the origin of vascular pathologies associated with air pollution exposures in humans and add to the growing literature suggesting that environmental exposure to volatiles, metals, and chemicals impact health by inducing gene expression changes in diverse cell and tissue types [64][65][66][67][68].