Mice were housed in the vivarium under standard conditions (12-hr light/12-hr dark cycle at 21°C) with food and water ad libitum. The project was designed according to the guidelines of the NIH, the European Communities Council (2010/63/EU) and the Italian Law for Care and Use of Experimental Animals (DL26/2014). It was also approved by the Italian Ministry of Health (authorization 510/2020-PR to EB) and the Bioethical Committee of the University of Turin. The study was conducted according to the ARRIVE guidelines.

To induce chronic EAE, 8 week-old female C57BL/6J mice (Charles River, Calco, Italy) were immunized by 2 subcutaneous injections of 200 μg myelin oligodendrocyte glycoprotein 35–55 peptide (MOG_35-55; Espikem, Florence, Italy) in incomplete Freund’s adjuvant (IFA; Sigma-Aldrich, Milan, Italy) containing 8 mg/ml Mycobacterium tuberculosis (strain H37Ra; Difco Laboratories Inc., Franklin Lakes, NJ, USA), followed by 2 intravenous injections of 500 ng of Pertussis toxin (Duotech, Milan, Italy) on the immunization day and 48 h later (36).

We exposed healthy (Ctrl) and EAE mice to a commercially available urban outdoor PM₁₀ (NIST Standard Reference Material SRM1648a; Sigma-Aldrich), which is a chemically standardized compound (https://tsapps.nist.gov/srmext/certificates/1648a.pdf), widely used in toxicological studies (37–39). Stock suspensions of NIST SRM1648 PM were prepared in sterile ultrapure water and aliquots were thoroughly mixed under sonication for 1 hour prior to each experiment. Mice were randomly divided into control (saline, i.e. NaCl 0.9% solution) and treatment (PM) groups. The treatment group was treated by intratracheal instillations (Intubation stand, Kent Scientific Corporation) of a PM₁₀ suspension (10 μg in 50 μl of saline), and the control group was treated with saline (50 μl), as in our previous study (40). We opted for the intratracheal administration of PM to assure low variability of absorption and exclude a direct nose-to-CNS transit of PM via the olfactory mucosa.

PM₁₀ dose was calculated considering the daily respiratory volume of mice (0.04 m³) and the daily peaks of PM₁₀ concentration in East Europe polluted areas (>250 μg/m³; 41). Therefore, in order to use a PM₁₀ dose relevant for human exposure, the dose used in this study was set as 10 μg (0.04 m³/day × 250 μg/m³). PM₁₀ was administered 4 days after EAE immunization, i.e. during the presymptomatic phase of the EAE disease course, where debilitating symptoms are yet to be displayed by the animals (i.e. clinical score=0, according to the conventional EAE clinical score assessment (36)). Peripheral blood of each mouse was collected in EDTA Vacutainer tubes (Becton Dickinson, New Jersey, USA), 6 hours after PM₁₀/saline exposure, and processed within 3 hours from sampling. Blood was centrifuged at 1300xg for 15 min to separate plasma.

For EV studies, we followed the MISEV2018 guidelines (42; Supplementary Table S1 , Supplementary Figure S1A, B ). EVs have been isolated from plasma pools from six mice of each condition: healthy mice which received saline (Ctrl); healthy mice which received PM₁₀ (Ctrl PM); EAE mice which received saline (EAE); EAE mice which received PM₁₀ (EAE PM). The resulting pooled samples were centrifuged three times at increasing speeds (1000xg, 2000xg, 3000xg) for 15 min at 4°C. After every centrifugation, the supernatant was decanted into a new tube and centrifuged again to remove cell debris and aggregates. Shortly after centrifugation, plasma samples were transferred into ultracentrifuge tubes (Polycarbonate Centrifuge Bottles, Beckman Coulter) and filled up to 10.5 ml with NaCl 0.9% solution. Ultracentrifugation was performed at 110,000xg for 75 min at 4°C to allow sedimentation of EVs. After ultracentrifugation, the supernatant was discarded and tubes were allowed to air dry for 2 min.

To confirm that this isolation protocol resulted in preparations enriched in EVs, western blotting (WB) was used to assess the expression of EV markers (i.e. the tetraspanins CD63/CD9) and absence of contamination from other cellular components (i.e. calnexin, a marker for endoplasmic reticulum; Supplementary Figure S1A ). Briefly, plasma EVs lysates were obtained in RIPA buffer (150 mM NaCl, 1% Triton X-100, 0.5% Sodium deoxycholate, 0.1% sodium dodecyl sulfate, 50 mM Tris Base, 5mM EDTA, 1mM EGTA) containing protease and phosphatase inhibitors (100mM DTT, 1M NaF, 200 mM sodium orthovanadate, 100 mM PMSF), boiled in SDS buffer and separated using a 10% SDS-PAGE gel. Proteins were transferred to PVDF membranes, blocked in 5% non-fat milk in TBS (20 mM Tris-HCl, pH 7.5, 150 mM NaCl) + 0.1% Tween-20 (TBS-T) for 1h and then incubated with primary antibodies (anti-rabbit CD63, 1:500, Immunological Sciences, Italy; anti-mouse CD9, 1:250, BD Biosciences; anti-rabbit calnexin 1:250, Cell signaling, USA) overnight at 4°C. After washes with TBS-T, membranes were incubated with appropriate secondary antibodies (anti-rabbit or anti-mouse, 1:3000, Sigma, Italy) for 1h at room temperature. Protein extract from a mouse cerebral cortex was used as a positive control. The chemiluminescent signal was visualized using Westar Antares ECL Blotting Substrates (Cyanagen; Italy), acquired with Bio-Rad ChemiDocTM Imagers (Bio-Rad; Italy).

Nanoparticle tracking analysis (NTA) and flow cytometry were used to characterize EV number/size ( Supplementary Figure S1B ) and cellular source, respectively (as in 43). To determine EV cellular origins, immunophenotyping was achieved with the MACSQuant Analyser flow cytometer (Miltenyi Biotec, Bergisch Gladbach, DE) following the manufacturer’s protocol. The Fluoresbrite Carboxylate Size Range Kit I (0.2, 0.5, 0.75, and 1 µm) was used to set the calibration gate on the MACSQuant Analyser system. To evaluate the integrity and to highlight EV subsets, 60-µL sample aliquots were stained with 0.02 µM 5(6)-carboxyfluorescein diacetate N-succinimidyl ester (CFSE) at 37°C for 20 min in the dark. CFSE is a vital non-fluorescent dye that enters EVs, where intracellular esterase enzymes remove the acetate group and convert the molecule into the fluorescent ester form. To characterize and count EVs, the following panel of antibodies was used: monoclonal APC-anti-CD14 (Clone REA934; dilution 1:10) to identify EVs from macrophages and/or monocytes; APC-antiCD41 (Clone: REA1194; dilution 1:10) for EVs from platelets; APC-anti-CD25 (clone: REA568; dilution 1:10) for EVs from regulatory T (Treg) lymphocytes. All antibodies were purchased from Miltenyi Biotec. Quantitative multiparameter analysis of flow cytometry data was conducted using FlowJo software (Tree Star, Inc., Ashland, OR, USA) and antibody gating strategies were performed as previously described (43).

miRNA isolation from EV pellets in tubes was performed using the miRNeasy Mini Kit (Qiagen, Venlo, NL) according to the manufacturer’s protocol, with the addition of the “Rneasy Mini Elute Clean-up” kit (Qiagen). To assess the quality of miRNA purification, samples were analyzed by 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) using Agilent RNA 6000 Pico Kit ( Supplementary Figure S2 ) and stored at -80°C.

Reverse Transcriptase reaction was performed using Megaplex™ RT Primers Pool A and B Rodent (Life Technologies, Foster City, CA) and TaqMan ^® MicroRNA Reverse Transcriptase Kit (Life Technologies). Briefly, two distinct reactions (A and B) were performed, to cover the reverse transcription of 758 target miRNAs, including 0.75 µl of Megaplex RT Pool, 0.15 µl of dNTPs (100 mM), 0.75 µl 10X RT Buffer, 0.90 µl of MgCl₂ (25 mM), 0.1 µl of RNases Inhibitor (20 U/µl) and 1.5 µl of MultiScribe™ Reverse Transcriptase (50 U/µl). The reverse transcription thermal protocol consisted of 40 cycles at 16°C for 2 min, 42°C for 1 min, and 50°C for 1 s, plus one cycle at 85°C for 5 min and final stage at 4°C. The cDNA samples were pre-amplified using “Low Sample Input Protocol (LSI) for profiling human miRNA using OpenArray^® Platform” (Application Note 2011- Life Technologies) with a specific Megaplex™ Preamp Primer Mix for Pool A and Pool B (Rodent). The pre-amplified samples were diluted 1:20 with nuclease-free water and were analyzed by QuantStudio™ 12K Flex Real-Time PCR System with OpenArray^® Platform (Applied Biosystems) according to manufacturer’s instructions.

The Gene Expression Suite Software was used to analyze the expression of 754 unique miRNAs. Four small RNA species (RNU48, RNU44, U6, and ath-miR159a), included in the array as internal technical controls to assess assay performance and RNA integrity, were also present but not used for normalization purposes. Only miRNAs with Relative Threshold Cycle (C_rt) value< 27 or AmpScore> 1.24 were considered amplified. The expression data were normalized using the “Global Normalization Factor” which calculates the average Crt of all detected miRNAs in a given sample. This method is recommended for EV-derived RNA, where validated endogenous reference genes are lacking due to the heterogeneity of vesicle content. Relative miRNA expression (also referred as Fold Change - FC) was determined using the relative quantification RQ= 2^(-Δcrt) formula, with -ΔC_rt = measured C_rt -medium normalization C_rt (44).

Statistical analyses were carried out with GraphPad Prism 9 (GraphPad software, Inc, RRID: SCR_002798). As for EV characterization (i.e. EV concentration, mean EV size and percentage of marker-positive EVs; Figure 1 ), the Shapiro-Wilk test was first applied to test for a normal distribution of the data. Then, as data were normally distributed, a Two-ways ANOVA test followed by Bonferroni’s post-hoc analysis was used for multiple group comparisons. P<0.05 was considered as statistically significant. Statistical differences were indicated with * P<0.05, **P<0.01, ***P<0.001, ****P<0.0001. The list of the applied tests and number of animals in each case are included in Supplementary Table S2 . Differential expression analysis of 754 unique miRNAs between groups was assessed using unpaired t-tests followed by the Benjamini-Hochberg FDR correction for multiple testing. A threshold of 0.20 was applied to the FDR p-value significance level to identify the set of top miRNAs.

The interactions between deregulated miRNAs and their target gene transcripts in each experimental condition were visualized using the web-based platform miRNet 2.0 (https://www.mirnet.ca/miRNet/home.xhtml 45). Significantly upregulated or downregulated miRNAs from each of the four groups, identified by their miRBase IDs, were uploaded into the “miRNA” module of miRNet 2.0. The target list was specified as “Genes (miRTarBase v9.0)” to link the input miRNAs to their predicted target gene transcripts, generating interaction tables and network visualizations. In the resulting networks, miRNAs may appear as specific strands (e.g., miR-9-5p from the input miR-9) or as closely related variants (e.g., miR-451a from miR-451). For the sake of clarity, only the most highly connected miRNAs in each network were labeled in the visualizations. KEGG pathway enrichment analysis was performed on the target gene lists from each network using the “Function Explorer” toolbox within miRNet, which employed a hypergeometric test to identify enriched (p-value < 0.05) pathways by comparing the input genes against the full set of mouse genes annotated in the KEGG database. The resulting p-values were adjusted for multiple testing using the Benjamini-Hochberg procedure to control the false discovery rate (FDR).

The list of target genes identified by the miRNA network analysis was used as input for a gene-to-disease association analysis using the DisGeNET Curated dataset. Genes associated with a disease class were displayed using the DisGeNET application in Cytoscape (46), operated through a DisGeNET specific R script. Edge thickness in the network graphical representation is proportional to the DisGeNET score for association robustness (S_GDA, score of Gene to Disease Association), which takes into consideration the number and type of sources and number of publications that support the association, following the subsequent formula:

To assess whether individuals primed toward the development of autoimmunity against CNS myelin (i.e. immunized to develop EAE) react differently to PM₁₀ exposure, we characterized plasmatic EVs (number/size, source and miRNAs cargo) in healthy (Ctrl) vs. presymptomatic (i.e. 4 days after immunization) EAE mice exposed to PM₁₀, compared to Ctrl and presymptomatic EAE mice which received saline ( Figure 1A ). Although not resulting in an overt motor impairment, nor in demyelination, in the presymptomatic stage of EAE the biological processes leading to the pathology have already started, as exemplified by cytokine and chemokine alterations in the nervous tissue, cerebrospinal fluid (CSF), and plasma (47), changes in the expression levels of molecules involved in neuronal-microglial communication (48), as well as by mouse cognitive impairment (49), thus recapitulating the early stages of MS.

EAE mice basally showed a slightly lower concentration of plasmatic EVs compared to Ctrl. PM₁₀ exposure induced a further significant decrease in plasmatic EV release in both Ctrl and EAE mice, compared to unexposed individuals ( Figure 1B ). Although reduced, about 1.2 x10⁹ EVs/ml were still found in the plasma of PM₁₀-exposed EAE mice ( Figure 1B ), showing a moderate increase in size, compared to those of unexposed EAE and PM₁₀-exposed and unexposed Ctrl mice ( Figure 1C ). To characterize EV cellular source, we exploited a panel of markers specific for cell types reported to be at the origin of plasma EVs upon PM exposure in humans, i.e. monocytes/macrophages and platelets (8, 15), and important players in MS pathophysiology, i.e. T regulatory cells (Tregs; 50). While in Ctrl mice 20-25% of EVs were produced by monocytes/macrophages and by Tregs, only 10% of EVs have these cellular origins in EAE mice, as assessed by immunopositivity for CD14 and CD25 surface antigens respectively ( Figure 1D ). Although the fraction of monocyte/macrophage-derived EVs appeared slightly higher in PM₁₀-exposed EAE mice compared to unexposed EAE, overall PM₁₀ exposure did not alter the percentage of marker-positive EVs in both Ctrl and EAE mice ( Figure 1D ).

Taken together, these analyses showed that – in terms of plasmatic EV number and source - mice primed to develop autoimmunity do not respond differently to PM₁₀ exposure, compared to healthy Ctrl mice.

To assess whether the response to PM₁₀ of presymptomatic EAE mice differed in terms of EV miRNA cargo, miRNAs have been extracted and quantified from plasmatic EVs. This analysis provided a list of deregulated miRNAs for each of the following comparisons: PM-exposed Ctrl (i.e. Ctrl PM) vs. unexposed Ctrl (i.e. Ctrl); unexposed EAE (i.e. EAE) vs. unexposed Ctrl mice; PM-exposed EAE (i.e. EAE PM) vs. unexposed EAE; PM-exposed EAE vs. PM-exposed Ctrl mice ( Figure 2 ). Deregulated miRNAs are listed in Supplementary Table S3 . In Ctrl mice, PM₁₀ exposure induced changes in the expression of 59 miRNAs, among which the most downregulated ones were miR-2134, miR-712-5p and miR-29a-5p, while the most upregulated were miR-215-5p, miR-199a-3p and miR-146b-5p (7 miRNA species with FDR<0.05; Figures 2A, E ). Of note, miR-2134 was also one of the most downregulated miRNAs among the 45 differentially expressed miRNAs in EAE vs. Ctrl mice, together with miR-1839-3p, miR-200b-3p and miR-451a. Top upregulated miRNAs in EAE vs. Ctrl mice instead included miR-223-3p, miR-1195 and miR-685 (6 miRNA species with FDR<0.05; Figures 2C, E ).

When comparing EAE mice exposed to PM₁₀ to unexposed EAE mice, we observed 26 differentially expressed miRNAs ( Figures 2B, E ). In particular, miR-192-5p, miR-451a and miR-122-5p were the most downregulated, whereas miR-9-5p and miR-680 were the most up-regulated. Finally, 94 miRNAs were differentially expressed between EAE and healthy (Ctrl) mice exposed to PM₁₀, with miR-200c, miR-451a and miR-25-3p among the most downregulated, and miR-9-5p, miR-1195 and miR142-3p among the most upregulated miRNAs (7 miRNA species with FDR<0.05; Figures 3D, E ).

Notably, a certain degree of overlap existed between the responses of healthy mice to immunization and to PM₁₀, especially in terms of downregulated miRNAs (58.84% of the miRNAs downregulated in Ctrl PM vs. Ctrl mice overlapped with those downregulated in EAE vs. Ctrl mice; Figure 3 ; Supplementary Table S4 ), suggesting a similar biological effect. In contrast, the responses to PM₁₀ of EAE and healthy mice were largely divergent, with only 22% of upregulated miRNAs and 12% of miRNAs downregulated in EAE PM vs. EAE overlapping with those downregulated in Ctrl PM vs. Ctrl mice ( Figure 3 ; Supplementary Table S4 ).

Overall, these analyses show that the response to PM₁₀ – in terms of deregulated EV-associated miRNAs – is qualitatively different in mice primed to develop autoimmunity compared to healthy Ctrl mice.

To further explore the potential functional impact of miRNA deregulations triggered by PM₁₀ in healthy and EAE mice, we identified the target gene transcripts of miRNAs that were differentially expressed in the 4 comparisons ( Supplementary Table S5 ). Notably, several targets were shared among multiple deregulated miRNAs, suggesting potential regulatory hubs. In line with this idea, network-based analyses performed on upregulated ( Figure 4 ) and downregulated ( Figure 5 ) miRNAs for each experimental comparison highlighted a significant number of interactions between the deregulated miRNAs (yellow squares) and their predicted target gene transcripts (red dots for upregulated miRNAs in Figure 4 and blue dots for downregulated miRNAs in Figure 5 ). Overall, the networks of upregulated miRNAs ( Figure 4 ) were remarkably more complex than those of the downregulated miRNAs ( Figure 5 ), which revealed relatively few target gene connections (especially for Ctrl PM vs Ctrl and EAE vs Ctrl, with fewer than 50 interactions; Figures 5A, C ). In contrast, the networks of downregulated miRNAs in EAE PM vs both EAE ( Figure 5B ) and Ctrl PM ( Figure 5D ) exhibited more intricate interactions, with a few miRNAs having over 300 connections each.

Examining the shared miRNAs within these networks allowed us to provide more insightful information about the actual biological consequences of miRNA deregulation, accounting for both the miRNA-target gene transcript interactions and their collective regulatory impact, rather than simply considering the shared miRNAs across conditions. Interestingly, among the upregulated miRNAs displaying the highest number of connections, a partial overlap appeared again between the response of healthy mice to immunization (EAE vs Ctrl) and to PM₁₀ (Ctrl PM vs Ctrl; e.g. miR-24-3p, and miR-223-3p among the upregulated miRNAs and miR-221-3p among the downregulated miRNAs; Figures 4A, C , 5A, C ), suggesting that these miRNAs may lead to similar cascades in response to both environmental and autoimmune triggers. Except for miR-1195, miR223-3p and miR-451a which were shared across EAE vs Ctrl and EAE PM vs Ctrl PM ( Figures 4B–D ), “hub” miRNAs deregulated in EAE mice responding to PM₁₀ were exclusive. Among them, miR-9-5p ( Figures 4B, D ) and miR-15a-5p ( Figures 5B, D ) were the most connected upregulated and downregulated miRNAs, respectively.

To further investigate the biological pathways potentially influenced by the detected miRNA deregulation, KEGG pathway enrichment analysis was performed on the predicted target genes from the upregulated and downregulated miRNA networks across the different comparisons. Of relevance in the context of an autoimmune CNS disease like MS, this analysis revealed a significant enrichment of terms related to synapse function, nervous system, and immune response, for both up- and downregulated miRNAs ( Figures 6A, B , respectively, Supplementary Table S6 ). Terms such as Glutamatergic synapse, GABAergic synapse, Dopaminergic synapse, Cholinergic synapse suggest that miRNA deregulation may play a role in modulating synaptic activity and plasticity, as further supported by the enrichment of Long-term potentiation and Long-term depression pathways. Moreover, several terms associated with addiction - such as Morphine addiction, Cocaine addiction, Amphetamine addiction, and Nicotine addiction - emerged across the comparisons, suggesting a possible influence of deregulated miRNAs on neurotransmitter systems. Besides nervous system-related pathways, immune- and inflammation-related pathways were also significantly enriched. These pathways were significantly enriched in the networks of upregulated miRNAs across the different comparisons ( Figure 6A ; Supplementary Table S5 ). In contrast, for downregulated miRNAs, nervous system- and immune-related pathways emerge almost exclusively in the comparisons involving EAE mice (i.e. EAE vs Ctrl, EAE PM vs Ctrl PM and EAE PM vs EAE comparisons). Although this might be due only to the very small size of the network of downregulated miRNAs in Ctrl PM vs Ctrl comparison ( Figure 5A ; Supplementary Table S5 ), this pointed again to the special character of the response of EAE mice to PM₁₀. In line with this idea, while the pathways described so far were largely shared across the four comparisons, several terms associated with fatty acid (FA) metabolism (Biosynthesis of unsaturated fatty acids, Fatty acid metabolism, Fatty acid elongation) were instead specifically enriched for the upregulated miRNA networks of EAE mice responding to PM (i.e. PM vs Ctrl PM and EAE PM vs EAE comparisons; Figure 6A ; Supplementary Table S6 ).

To link EV-packaged miRNA dysregulation with possible pathological alterations elicited by PM exposure, we used the list of predicted target gene transcripts to perform a “gene-to-disease association” analysis. Again of relevance in the context of MS, and consistent with the KEGG pathway analysis, many of the targets of the miRNAs selectively deregulated in PM-exposed EAE vs PM-exposed healthy mice (i.e. EAE PM vs Ctrl PM; see Supplementary Table S5 ) are involved in nervous systems diseases and immune system diseases, and particularly in autoimmune diseases including MS ( Figures 7A–D ; Supplementary Table S7 ). Consistently, when assessing the fraction of gene transcripts targeted by miRNAs selectively deregulated in PM-exposed EAE mice over the entire number of genes listed in each of the disease classes included in the DisGeNET database, “Nervous System Diseases” was the most represented class, followed by Mental Disorders, Hemic and Lymphatic Diseases, Behavior and Behavior Mechanisms and Immune System Diseases ( Figures 7E, F ).

Overall, this analysis showed that gene transcripts targeted by miRNAs selectively altered in EAE mice exposed to PM₁₀ were associated with nervous system disorders and with pathologies closely related to MS, suggesting that the plasmatic EV cargo may mediate the transfer of PM-induced damage to the CNS and synergize with the immune setting of the individual, possibly contributing to immune dysregulation in the CNS via yet unknown mechanisms.

Overall, our findings support the view that, in immunologically primed individuals - such as people with MS - exposure to PM₁₀ can trigger a unique response through the delivery of EV-associated miRNAs, which can contribute to a worse manifestation or an unfavorable evolution of the pathology. Yet, our study has limitations that must be disclosed to correctly interpret our findings and to identify persisting open issues. Here, we combined exposure to PM₁₀ with EAE, which is the closest animal model that approximates human MS pathogenesis by reflecting its autoimmune component (104). Yet, as our results point to a likely participation of exposure to PM₁₀ in immune dysregulation, it is relevant to acknowledge that some differences exist in the innate and adaptive immune system of mice and humans (105), potentially limiting the translatability of our findings and calling for future investigations in MS patients. Moreover, besides the relatively small size of individuals used in this kind of in vivo experiments, here we used only female mice, as in most studies focusing on EAE. This choice relies on the well-known sexual dimorphism existing in MS, with women showing a higher susceptibility to the disease (106). Yet, disease progression appears worse in men (106) and mechanisms leading to EAE partially differ in female and male mice (107). Inclusion of male mice in future studies will validate the generalizability of our findings or, alternatively, unveil sex-specific aspects in the response to PM₁₀ of immunologically primed individuals. Moreover, as regards the analysis of EV-packaged miRNAs, this study did not include EV subtype-specific miRNA profiling due to limited EV yield, and sample pooling - while necessary due to sample constraints and ethical (i.e. 3Rs) considerations - limited the assessment of inter-individual variability. Finally, our study investigated an acute (i.e. occurring within 6 hours) biological response triggered by a single exposure to PM₁₀. While this was instrumental to investigate the short-term effect of the peaks of airborne PM₁₀ which were associated with higher rates of hospital admissions for MS and MS relapses (22, 29), results obtained in this experimental setting might not reflect the biological responses elicited upon long-term or chronic exposure to PM₁₀. Similarly, in our study PM₁₀ dose was calculated to mimic daily peaks of PM₁₀ concentration in the most polluted European areas (>250 μg/m³; 41). Whether exposure to higher PM₁₀ concentrations - as detected in other geographical regions (108) - or to other risk factors epidemiologically associated with MS risk may converge on the same biological mechanisms unveiled by our study deserves further investigation.