Microarray analysis of mRNA expression profiles in liver of ob/ob mice with real-time atmospheric PM2.5 exposure

Epidemiological studies have demonstrated the association between exposure to fine particulate matter (PM2.5) and the onset of non-alcoholic fatty liver disease (NAFLD). However, the potential biological mechanism is largely unknown. Our study was aimed to explore the impact of PM2.5 on the transcriptome level in the liver of ob/ob mice by atmosphere PM2.5 whole-body dynamic exposure system, and meanwhile preliminarily investigated the effects of metformin intervention in this process. More than three thousand differentially expressed genes (DEGs) was screened out by microarray analysis (p < 0.05, |FC|> 1.5). KEGG pathway enrichment analysis showed that these DEGs were mainly enriched in cancers, infectious diseases, and signal transduction, and the most significant pathways were thyroid hormone signaling pathway, chronic myeloid leukemia and metabolic pathways. Then, 12 hub genes were gained through weighted gene correlation network analysis (WGCNA) and verified by qRT-PCR. The expression of 5 genes in darkslateblue module (cd53, fcer1g, cd68, ctss, laptm5) increased after PM2.5 exposure and decreased after metformin intervention. They were related to insulin resistance, glucose and lipid metabolism and other liver metabolism, and also neurodegenerative diseases. This study provided valuable clues and possible protective measures to the liver damage in ob/ob mice caused by PM2.5 exposure, and further research is needed to explore the related mechanism in detail.


Introduction
Fine particulate matter (PM 2.5 ) refers to particles with an aerodynamic diameter (AED) of less than 2.5 μm in the ambient atmosphere. PM 2.5 generally comes from natural sources and man-made production, while the sources in cities are mainly industrial emissions and traffic-related exhaust particles (Mukherjee and Agrawal 2018). It was estimated that more than 92% of people in the world living in places with PM 2.5 concentrations higher than the 2005 World Health Organization (WHO) air quality guideline of 10 μg/m 3 . Lately, according to the findings in recent 15 years, WHO released the updated air quality guideline in which the reference value of PM 2.5 was limited as 5 μg/m 3 (WHO 2021b). An investigation of 47 representative cities from six continents revealed that PM 2.5 concentrations in only 2 cities complied with the new guideline, while in places from underdeveloped countries, PM 2.5 concentrations needed more than 90% decrease, which reflected that the majority of global population would live in the areas with PM 2.5 concentration above WHO guideline reference value in a rather long term in the future, especially Responsible Editor: Ludek Blaha Lisen Lin and Li Tian contributed equally to the work. in the low-income and middle-income countries (Carvalho 2021). Due to the special physical and chemical properties, a variety of chemical substances are easily absorbed to the particle surface of PM 2.5 , including metals, polycyclic aromatic hydrocarbons (PAHs), and endotoxins, which are typical toxicant that could give rise to gene mutations, DNA damage, and epigenetic changes (Wu et al. 2017). Strong and accumulative evidence demonstrated the causal relationships between PM 2.5 exposure and all-cause mortality, pulmonary, and cardiovascular diseases (WHO 2021a). The major exposure route of PM 2.5 in human is through inhalation, but Liang et al. observed the deposition of PM 2.5 not only in the alveolar region of lung as generally recognition but also in extra-pulmonary organs such as liver and kidney detected by a fluorescent imaging method in vivo . As the hazardous health effects of PM 2.5 have been raised global concern for nearly 30 years, the detrimental effects of PM 2.5 beyond respiratory and cardiovascular system have drawn more and more attention in recent years.
Non-alcoholic fatty liver disease (NAFLD) is a spectrum of liver diseases from simple non-alcoholic fatty liver, nonalcoholic steatohepatitis (NASH) to irreversible cirrhosis; it was the most prevalent chronic liver disease which would gradually progress to extrahepatic cancer, hepatocellular carcinoma, cirrhosis, or cardiovascular disease and was associated with overall and cause-specific mortality (Simon et al. 2021). Lately, a 16-year prospective cohort study reported that long-term PM 2.5 exposure was associated with higher risk of NAFLD in 58,026 participants when the concentrations of PM 2.5 exceeded 23.5 μg/m 3 ; they found that each 1 μg/m 3 elevation in PM 2.5 concentration was related with an HR of 1.06 for NAFLD (Sun et al. 2021). This study provided the first epidemiology evidence to illustrate that PM 2.5 exposure could be an important environmental risk factor for NAFLD. Zheng et al. observed that 10-month exposure of real-world PM 2.5 induced NASH-like phenotype in mice with calculated mean daily exposure concentration at 11.6 μg/m 3 (Zheng et al. 2013). The subsequent animal studies demonstrated that PM 2.5 exposure decreased hepatic glycolysis, the Krebs cycle and GSH synthesis, increased hepatic lipogenesis, thus disrupting redox balance in the liver, gradually causing inflammation and lipid steatosis in liver . These studies revealed that oxidative stress, inflammation, insulin resistance, and circadian rhythm played important roles in contributing the imbalance of liver metabolism. Even so, the related mechanism and pathways involved in this process were largely unknown.
Obesity population, as susceptible population, may be at a greater risk of air pollution. A large number of epidemiological studies have shown that obese population are at higher risk of cardiovascular disease and Alzheimer's disease when exposed to particulate matter (Weichenthal et al. 2014). Meanwhile, obesity is one of the most important driving factors for liver diseases. It has been shown that the prevalence rate of steatohepatitis was approximately 3% in non-obese persons, 20% in obesity population, and 40% in extremely obese patients (Fabbrini et al. 2010). The risk of liver cancer in obese patients was 1.4-4.1 times of that in healthy people (Huang et al. 2021). Hiesh et al. observed that traffic-related air pollution was associated with serum level of cytokeratin-18, the indicator of NASH risk, in 74 overweight and obese adolescents (Hsieh et al. 2018). The ambient PM exposure combined with high-fat diet treatment could cause a synergistic effect on the changes of lipid accumulation, oxidative stress, and inflammation in the mouse liver (Ghassabian et al.). However, currently known information about the effects of PM 2.5 on the liver in obesity is limited, and the regulation mechanism is not yet clear. Transcriptomics based on gene chips is a reliable approach to provide comprehensive understanding on the changes in mRNA level caused by PM 2.5 exposure. The screened differentially expressed mRNA will provide clues and basis for follow-up research.
Metformin is a universal first-line medication for treating type 2 diabetes, and it was reported that metformin could restrain gluconeogenesis in liver in a redox-dependent manner (Madiraju et al. 2018). Haberzettl et al. found that metformin could prevent PM 2.5 -induced vascular insulin resistance and activation of NF-κB and inflammasomes thus preventing endothelial progenitor cells (EPCs) mobilization and restoring EPCs levels, so as to maintain EPCs homeostasis (Haberzettl et al. 2016). A number of studies have concluded that metformin can reduce the risk of liver cancer by about 50%, indicating that it had hepatic protection effect while the mechanism was not clear (Huang et al. 2021). It was found that metformin could effectively alleviate hyperglycemia in obese mice by stimulating CBP (CREB-binding protein) phosphorylation to block insulin signaling pathway and suppress hepatic gluconeogenesis (He et al. 2009). However, the effect of metformin on PM 2.5 -induced liver toxicity in obesity has not been reported so far. The health interventions, like metformin, should be investigated to reduce the undesirable environmental pollution effect.
Therefore, the aim of this study is to preliminarily explore the influence of PM 2.5 on the transcriptome level in the liver in ob/ob mice and to discover whether metformin has hepatic protection under exposure of PM 2.5 . The study could provide a basis for the screening of biomarkers and the mechanisms after atmospheric PM 2.5 exposure.

Animal treatment and real-time whole-body PM 2.5 exposure
Twenty male ob/ob mice (C57BL/6 J background, Huafukang Bio-Technique Co., Ltd, Beijing, China) were purchased from the Animal Experimental Center of Capital Medical University and acclimated for a week before the experiment. Twenty mice were divided into four groups randomly and the interaction experiment was designed with two interaction factors (PM 2.5 and metformin). The control group (Con group) was exposed to filtered air and drank pure water, the PM 2.5 exposure group (PM 2.5 group) was exposed to concentrated PM 2.5 , the drug group (Met group) drank prepared metformin solution, and the intervention group (PM 2.5 + Met group) has both factors. The metformin concentration (250 mg/kg/day) in water was determined on a per-mouse basis and adjusted daily, based on measured daily water intake and body weight (Luo et al. 2016).
Mice were housed four to five per cage on corncob bedding with ad-lib access to food and water. The humidity was 50% and the temperature was 22-26 °C with a 12 h light/ dark cycle. The PM 2.5 exposure was carried out in Hinnerstype stainless-steel whole-body inhalation chambers and concentrated by the ambient particulate matter whole-body dynamic exposure system; the cleaned air was filtered by animal cage air filter. The exposure period was 6 h per day, 6 days per week from November 14th, 2019 to December 11th, 2019 at a total of 4 weeks, then the mice was sacrificed. The experimental protocol was approved by the Committee of the Ethics Animal Experiments of Capital Medical University (AEEI-2019-161) and carried out under the institutional guidelines for ethical animal use.

The exposure equipment and parameters monitor
The real-time PM 2.5 concentrated exposure was conducted by small animal whole-body dynamic exposure system (HRH-300L, Beijing Huironghe Technology Co., Ltd. Beijing, China). The inhalation exposure chambers were outfitted with air quality monitor and aerosol generator to concentrate PM 2.5 , which ensure that the PM exposure is consistent with the changes of the external environment, and the concentration is not too low (concentrate 6-10 folds).
The condition inside the chambers was closely monitored to maintain a relatively constant 20-25 °C temperature, 40-60% humidity, and 18-20/h ventilation frequency. The ambient PM 2.5 concentration was monitored by using the Aerosol Detector Dusttrak Drx Aerosol Monitor 8533 (TSI Instrument, Shoreview, MN). The characteristics of particles were measured by an Aerodynamic Particle Sizer (APS) Spectrometer 3938 N (including particle number, particle surface, particle mass) and analyzed by Aerosol Instrument Manager Software Version 10.1 (TSI Instrument, Shoreview, MN).

Total RNA extraction and microarray analysis
For Affymetrix microarray profiling, the total RNA of mice liver was isolated by using TRIzol reagent (Invitrogen, Carlsbad, Canada), then purified with RNeasy Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol. Finally, the amount and quality of RNA were determined using a UV-Vis spectrophotometer (Thermo, Nan-oDrop 2000, USA) at an absorbance of 260 nm. Each group had three replicate samples tested, and for each sample, the experiment was performed in triplicate as technical replicates. The mRNA expression profile was measured using Clariom ™ S Assay (Affymetrix GeneChip, USA). Gene-Chips were washed and stained in the Affymetrix Fluidics Station 450. All arrays were scanned by using Affymetrix® GeneChip Command Console (AGCC) which was installed in GeneChip® Scanner 3000 7G. The microarray analysis was performed using Affymetrix Expression Console Software (version 1.2.1). The row data(.cel file) were normalized by the software TAC (Transcriptome Analysis Console; Vension:4.0.1) with Robust Multichip Analysis (RMA) algorithm using Affymetrix default analysis settings and global scaling as a normalization method. Values presented are log2 RMA signal intensity. Data from the microarray analysis (ECL files) discussed in this article were deposited in the National Center for Biotechnology Information (NCBI). All microarray data is MIAME compliant and the raw data has been deposited in NCBIs Gene Expression Omnibus (NCBIs GEO ID: GSE186900, https:// www. ncbi. nlm. nih. gov/ geo/ query/ acc. cgi? acc= GSE18 6900).

Differentially expressed genes (DEGs) analysis
In microarrays, the R package "limma" was used to filter the differentially expressed genes (DEGs). R package "limma" used moderated F-statistic to filter the multi-group differentially expressed genes. The p-values were corrected by empirical Bayes moderation. Benjamini-Hochberg was used for multiple tests correction (FDR was used to adjust the p-values for multiple comparisons). The threshold sets (fold change > 2.0, p-value < 0.05, and FDR < 0.05) were conducted to filter these up and downregulated genes.

Gene Ontology (GO) enrichment analysis and pathway enrichment analysis
Gene Ontology analysis is an internationally standardized system that classifies the gene function and provides a series of dynamically controlled vocabulary to comprehensively describe the attributes of genes and gene products in organisms. After the GO function classification annotation of DEGs given, the GO enrichment analysis of DEGs was conducted. GO has three ontologies: molecular function, cellular component, and biological process. The GO function enrichment analysis can determine the main biological functions performed by the DEGs.
In organisms, different genes coordinate with each other to perform their biological functions. Pathway analysis helps to further understand the biological functions of genes. KEGG is the main public pathway database. Pathway enrichment analysis takes KEGG pathway as the unit and applies hypergeometric test to find pathways that are significantly enriched in DEGs compared with the background of the entire genome.
Because the basic unit of GO is term, map the DEGs to each term of the GO database (http:// www. geneo ntolo gy. org/) and calculate the number of genes in each term. Then apply hypergeometric test to find GO terms that are significantly enriched in DEGs compared with the whole genome background. After the calculated p-value is corrected by FDR, GO terms with corrected p-value ≤ 0.05 are defined as GO terms that are significantly enriched in DEGs. Pathway enrichment analysis is the same as GO enrichment analysis, but needs to be replaced with KEGG pathway database.

Trend analysis and series test of cluster of Gene Ontology (STC-GO) analysis
Trend analysis was to cluster the gene expression pattern based on the characteristics of multiple continuous samples. The same gene expression pattern can be found through trend analysis. First of all, all terms enriched by GO/KEGG terms were identified and hierarchically clustered the terms statistically based on similarities among their gene memberships. Then a subset of representative terms was selected and converted into a network layout. More specifically, all analysis has been carried out with the following ontology sources: GO biological processes, KEGG pathway, Reactome gene sets, CORUM, TRRUST, PaGenBase, Wiki pathways, and PAN-THER pathway. All genes in the genome have been used as the enrichment background. Terms with a p-value < 0.01, a minimum count of 3, and enrichment factor > 1.5 are collected and grouped into clusters based on their membership similarities (similarity > 0.3 is considered a cluster).

Weighted gene correlation network analysis (WGCNA)
The genes evaluated for availability and the gene coexpression network was constructed by using the R package "WGCNA" in R (version 4.1.1). The strength of correlation between genes was tested by Pearson's correlation coefficient. The adjacency matrix was constructed to describe the correlation strength between the gene nodes, and further transformed into a topological overlap matrix (TOM). The formulas were as follows (Langfelder and Horvath 2008): Then, the TOM matrix is used to quantitatively describe the similarity in gene nodes by comparing the weighted correlation between two gene nodes and other gene nodes. WGCNA identifies gene modules using unsupervised clustering, then the modules of similar gene composition were identified, the characteristic genes were calculated, the modules were hierarchically clustered, and similar modules were merged.
The co-expression modules analyzed by WGCNA are gene clusters with high topological overlap similarity. Genes in the same module have a higher correlation and degree of co-expression. The module eigengene E (ME), as the first principal component of the module, is used to represent the expression pattern of the genes of the module in each sample.

Quantitative RT-PCR
The intramodule connectivity of a gene is equal to the sum of the degree of correlation between genes in that module. Therefore, the top 5 genes with the highest intramodule connectivity (the highest degree in each module) were selected as hub genes, and qRT-PCR were used to verify the expression. Notably, in the section of qRT-PCR validation, the same liver samples were used as those in the genechip, and each sample was repeated independently for 3 times to ensure random error. According to the protocol provided by the manufacturer, the Direct-zol RNA Mini-Prep kits (R2050, ZYMO) were used to extract total RNA. PrimeScrip RT reagent Kit (RR037A, Takara, Japan) and SYBR Premix Ex Taq II (Tli RNaseH Plus) (RR820B, Takara, Japan) are used for reverse transcription and amplification, respectively. Liver samples from 5 animals were used in each group. GAPDH was used as an internal reference, and the 2 −ΔΔCt value was normalized to its expression level. The sequence of qPCR primers used in this study is placed in Supplementary Table 2 (Table S2). The experiment was performed in triplicate.

Statistical analysis
The OmicShare online analysis tools (http:// www. omics hare. com/ tools), R, Cytoscape, and Metascape were performed in the part of bioinformation analysis (Zhou et al. 2019). GraphPad Prism 8.0 was used to test and present the results of qRT-PCR. p < 0.05 was considered statistically significant.

Concentration and characterization of PM 2.5 in exposure chamber
The experiment was carried out in November and December, and the PM 2.5 pollution in Beijing showed strong spatiotemporal variations. Day to day variation of PM 2.5 possessed a long-term trend of fluctuations, with 2-6 peaks each month (Huang et al. 2015). The in-chamber monitoring data showed that the PM 2.5 concentration had been fluctuating dramatically during the 28-day poisoning process (Fig. S1). The maximum value of PM 2.5 concentration was 293 μg/m 3 with the minimum value was 18 μg/m 3 , and the average concentration was 164 μg/m 3 (Table S1).
The particle cutter equipped with the system can block the entry of most large particles and cut the large particles into fine particles (AED < 2.5 μm, aerodynamic equivalent diameter), even the ultrafine particles (AED < 0.1 μm). After measurement, it can be seen that the almost all PMs are less than 1.0 μm, indicating that the system is operating reliably. As for the particle surface area, most of the particles are also concentrated below 1.0 μm, and there are almost no particles above 10.0 μm. The particle mass increased sharply between 2.5 and 10 μm, probably because although there is a small amount of PM 10 , the weight of PM 10 is quite considerable compared to PM 2.5 (Table S5).

Global differentially expressed genes expression in liver tissues
A total of 22,207 genes were detected using the chip, of which 3574 differentially expressed genes (DEGs) were screened out using the set conditions (p < 0.05, FC > 1.5). The cluster heat map shows the relative expression of differential genes after PM and metformin intervention (Fig. 1A). It can be seen that the expression levels of most genes in the control group and the PM 2.5 group are opposite, and some genes and expression levels in the MP intervention group have recovered to a certain extent. It should be noted that the two samples in the MP group were eliminated due to the large error and the poor clustering effect. Ranked by statistical significance, the top ten gene symbols are aacs, gm17530, synrg, rassf6, anks4b, pnpla3, arntl, ddi2, angptl8. The relative gene expression level of each group, fold change (FC), p-value, adjusted p-value, and gene description are listed in detail in the table (Table 1).

Gene Ontology enrichment analysis and pathway enrichment analysis of differentially expressed genes
The GO enrichment analysis showed that the differentially expressed genes were mainly enriched in protein binding, transferase activity, nucleotide binding, ATP binding, and RNA binding in terms of molecular function; nucleus, cytosol, cytoplasm, membrane, and nucleoplasm in terms of cellular component; while regulation of transcription, DNA-templated transcription, DNA-templated metabolic process, negative regulation of transcription from RNA polymerase II promoter, and protein transport in terms of biological process (Fig. 1B-D, Table 2).
The pathway enrichment analysis was also conducted by using KEGG database, and found that these differentially expressed genes were mainly involved in thyroid hormone signaling pathway, chronic myeloid leukemia, metabolic pathways, HTLV-I infection, and AMPK signaling pathway. Among them, thyroid hormone signaling pathway has the smallest p value, indicating that it is most likely to change the signal transduction of thyroid hormone; metabolic pathways have the smallest RichFactor with the most enriched genes, which indicated that metabolic pathways was the most important and meaningful pathway for differential genes ( Fig. 2A, Table 3).
The cluster analysis of the enriched pathways found that the pathways ranked from high to low belong to human diseases, metabolism, organismal systems, environmental information processing, and cellular process by the top pathway classification, while cancers, infectious diseases, signal transduction, endocrine system, and carbohydrate metabolism by the middle pathway classification, respectively ( Fig. 2B and C).

Trend analysis and STC-GO analysis of differentially expressed genes
Trend analysis showed that all differentially expressed genes were fitted to 26 trends, of which 9 trends highlighted by color were statistically different. Red represented a gene set with an overall upward trend, green represented a gene set with an overall downward trend, and black represented a gene set with an uncertain trend (Fig. 3A). The number of genes contained in each fitted trend and the specific expression trend of each gene are presented in the figure (Fig. 3B).
Trend analysis was used to find the trends of opposite expression between the PM 2.5 -exposed group and the metformin intervention group. Among them, trend 6 (contains 60 genes) and trend 15 (contains 39 genes) met the requirement. Therefore, these two trends were separately proposed to use STC-GO analysis. The results showed that the main functions of the trend 6 genes were pentose phosphate pathway oxidative phase glucose 6P to ribulose 5P, staphylococcus aureus infection, and monocarboxylic acid metabolic process (Fig. 4A, Table S3) and trend 15 was mainly on HTLV-I infection, antigen processing and presentation of exogenous peptide antigen via MHC class II, and cytokinemediated signaling pathway (Fig. 4B, Table S4).
Through STC-GO analysis, the relevance of these gene functions has been clearly displayed. In trend 6 STC-GO network, the largest and most important cluster was composed of carbon metabolism, monocarboxylic acid metabolic process, pentose biosynthetic process, and small molecule biosynthetic process. The relationships in trend 15 were also shown in clusters with non-alcoholic fatty liver disease (NAFLD), pertussis, and staphylococcus aureus infection (Fig. 4C and D).

Results of WGCNA
First of all, a gene correlation heatmap was constructed to measure the all DEGs co-expression network and to make use of interaction patterns among genes (Fig. 5A). Then the hierarchical clustering Dynamic Tree Cut was used to identify modules (Fig. 5B). Through analysis of the scale independence and mean connectivity, the soft threshold power was determined (Fig. 5C and D). After being raised to a suitable height, we got five modules with different expression trend (ME darkmagenta, ME red, ME darkslateblue, ME antiquewhite4, ME darkgrey). Module-trait relationships heatmap indicated the expression level and p-value of different modules in different groups (Fig. 5E) and the eigengene adjacency heatmap could intuitively reflect the relationship between the modules (Fig. 5F). Figure 6A showed all the DEGs in the five modules and their relationships. Five genes with top degree in each module was selected as the hub genes, and the bar graph showed the expression of these 25 hub genes. In order to find biomarkers that can monitor and predict diseases, we use the WGCNA analysis to find out those important genes (hub genes) in the entire gene co-expression network. The WGCNA clustering criterion has a great biological significance (Fig. 5). Due to the unique soft threshold algorithm of WGCNA, the gene expression network tended to be distributed with free-scale network, which made the results have higher reliability (Tian et al. 2020). This analysis divided all DEGs (p < 0.05, FC > 1.5) into 5 modules (Fig. 6A). However, after verification and screening, only 12 genes in 4 modules were met the criteria ( Fig. 6B and C).

Validation of gene chip results by quantitative RT-PCR
To verify the WGCNA analysis results, the expression of 25 hub genes in the liver of ob/ob mice was detected using qRT-PCR. Repeated independent experiments for three times, there are 12 hub genes with stable trends and consistent with the WGCNA analysis results (Fig. 6C, Table 4). The results of verification showed good consistency with WGCNA analysis. In darkslateblue cluster, PM 2.5 exposure upregulated the gene expressions of cd53, fcer1g, cd68, ctss, and laptm5, and the effects were alleviated by metformin treatment. The gene expressions of mup6, mup8 in the red cluster, and sub1, snrpd2, etohi, zfp931 in the darkmagenta cluster were reduced in response to PM 2.5 exposure, the gene expressions of sub1, etohi1, and zfp931 were significantly reduced (p < 0.01), while metformin did not restore the decrease of the gene expressions. The gene expression of egln1 in antiquewhite4 cluster was also downregulated after PM 2.5 exposure while metformin treatment totally regained the depression.

Discussion
Air pollution has become a severe environmental problem all over the world. Therefore, the study explored the changes in the liver of obese mice at the transcriptome level under the cross-intervention of PM 2.5 and metformin. Through the analysis of 12 hub genes, we discussed the  pathophysiological and functional changes that may occur in the liver, especially the physiological processes and diseases associated with the liver, such as thyroid function, insulin resistance, and lipid metabolism; NAFLD, HCC, even neurodegenerative diseases. The study could also provide some biomarkers for early screening of diseases caused by PM 2.5 . GO function analysis (Fig. 1) showed that the most obvious term in molecular function is "protein binding." Impairment of liver function may not only disturb liver metabolism but also affect plasma protein binding, which in turn affects the distribution and removal of metabolites in the body (Verbeeck 2008). What needs attention in cellular component was nucleus, cytosol, and cytoplasm, which may be related to cytosol-nucleus traffic and colocalization in hepatocytes (Romanque et al. 2011). Transcription and metabolism were more important in biological processes. This may be because the liver uses a series of liver transcription factors to regulate the expression of genes involved in all aspects of lipid metabolism (including catabolism, transportation, and synthesis) (Karagianni and Talianidis 2015).
According to the KEGG pathway analysis, PM 2.5 could affect the thyroid hormones signal pathway with most statistically different, even thyroid cancer (Fig. 2). Thyroid hormones (THs) had a significant effect on the anabolism of fatty acids and cholesterol in the liver, and direct regulate de novo lipogenesis, tricarboxylic acid cycle (TCA), fatty acid β-oxidation, OXPHOS, lipolysis, and lipophagy pathway, which mainly involves genes, such as acc1, me, fasn, thrsp, cpt1a, pdk4, mcad, ucp2, hmgcl, atgl (Sinha et al. 2018). It was reported that low thyroid function in population is associated with increased likelihood of chronic fibrotic diseases of the liver (Bano et al. 2020). Population studies found that prenatal exposure to PM 2.5 can damage neonatal thyroid function (Ghassabian et al. 2019). Our research found that PM 2.5 upregulates the thyroid signaling pathway in the liver, which was also proved from the level of metabolic organs (Kim et al. 2020). The risk of NAFLD was inversely correlated with free thyroxine levels (Ritter et al. 2020). THs modulated the homeostasis of hepatic lipid metabolism by regulating lipoprotein, triglyceride (TAG) storage, and cholesterol levels, which had a key effect on liver-related diseases, such as NAFLD and hypercholesterolemia (Martínez-Sánchez et al. 2017) and THs may modulate co-activators and co-repressors through the hypothalamic-pituitary-thyroid axis, thereby altering cholesterol metabolism in the liver (Ritter et al. 2020). In addition, our results showed that the pathway with the most DEGs enrichment was the "metabolic pathways" (Fig. 2). The liver played an important role in glucose homeostasis by controlling various pathways of glucose and lipid metabolism, including oxidation, gluconeogenesis, and adipogenesis (Han et al. 2016). It also regulated other important metabolisms, including purines and pyrimidine synthesis, histidine catabolism, methionine recycling, and formic acid utilization (Zaitsev et al. 2019). The results showed that the enriched related pathways were AMPK signaling pathway, PPAR signaling pathway, insulin resistance, and hepatocellular carcinoma. These pathways were involved in the pathophysiological process of oxidative stress, inflammation, abnormal metabolism leading the accumulation of glucose, and lipid in the liver ). The excessive production and accumulation of hepatic lipid might induce liver fibrosis in further, which was in line with the evolution of NAFLD. Animal studies have confirmed that PM 2.5 can induce excessive extracellular matrix accumulation in liver tissues and eventually lead to liver fibrosis, which was a foreshadow to liver cancer (Zheng et al. 2015). A prospective epidemiological study in the USA shown that environmental PM 2.5 exposure may be a risk factor for HCC (VoPham et al. 2018). So, the disturbance of metabolic pathways might be the first step of PM 2.5 -induced liver injury and the long-term hazardous hepatic effects of PM 2.5 exposure could be overwhelming.
In order to understand the changes in DEGs, we conducted a trend analysis and STC-GO analysis and obtained two reasonably interpretable trends (Figs. 3 and 4). As can be seen from the plot, the genes in trend 6 increased after being exposure and decreased with metformin ( Fig. 4A and C). The term "GO 0019322: pentose biosynthetic process" has the most significant difference. Study observed an increasement in the oxidative branch of the pentose phosphate pathway and 13 C incorporations suggestive of enhanced capacity for the de novo synthesis of fatty acids,

Fig. 3
Total differential mRNA expression trend by p-value (A) and total differential mRNA expression trend profile line (B) in the STC analysis. When the trend is statistically significant (p < 0.05), the leg-end is colored and the red legend indicates an upward trend, the green legend indicates a downward trend, and the black legend indicates an amorphous trend which indicates an increase in insulin resistance (Reyes-Caballero et al. 2019). On the other hand, metformin can relieve the pentose phosphate pathway, inhibit gluconeogenesis, and promote glycogen retention to reduce insulin resistance (Atangwho et al. 2014). Another term "KO 04932: Non-alcoholic fatty liver disease (NAFLD)" also involved gluconeogenesis, glycogen, and insulin resistance, which was in accordance with the above research results. We also found other interesting pathways, such as "KO 05150: staphylococcus aureus infection," "GO 0046942: carboxylic acid transport," "GO 0007188: adenylate cyclase-modulating G protein-coupled receptor signaling pathway" and there were also many articles confirmed that these terms are related with liver or PM 2.5 exposure. For example, monocarboxylate transporter 1 (MCT1) expression was downregulated in adipocytes of diabetic rats thus to impair the ability to transport lactic acid (Hajduch et al. 2000); metastasis and glycolysis could be induced by upregulation of MCT1 expression and subsequently activating Wnt/β-catenin signaling pathway in HCC (Fan et al. 2018). The expression of trend 15 genes only increased significantly after PM 2.5 exposure, and metformin intervention does not change the rising trend ( Fig. 4B and D). The terms "WP 447: Adipogenesis genes" and "GO 0046890:  Results of weighted correlation network analysis. A Network heatmap plot. Branches in the hierarchical clustering dendrograms correspond to modules. Color-coded module membership is displayed in the color bars below and to the right of the dendrograms. In the heatmap, high co-expression interconnectedness is indicated by progressively more saturated yellow and red colors. Modules correspond to blocks of highly interconnected genes. Genes with high intramodular connectivity are located at the tip of the module branches since they display the highest interconnectedness with the rest of the genes in the module. B Gene dendrogram obtained by average linkage hier-archical clustering. The color row underneath the dendrogram shows the module assignment determined by the Dynamic Tree Cut. Determination of soft-threshold power in the WGCNA by C analysis of the scale-free index for various soft-threshold powers (β) and D analysis of the mean connectivity for various soft-threshold powers. E and F Heatmap of the relationship between modules and correlation of the five module eigengenes. In the heatmap, green color represents low adjacency (negative correlation), while red represents high adjacency (positive correlation) regulation of lipid biosynthetic process" are noteworthy and need further study. PM 2.5 exposure could increase the expression of genes related to lipid synthesis through different mechanism. For example, srebp1 was involved in regulating the expression of fasn, acc, and scd1; exposure to PM 2.5 resulted in increased expression of bmal1, rev-erbα, and pparα, affecting circadian rhythm, liver triglyceride, free fatty acid levels, or fatty acid transport (Yan et al. 2020). It can be seen from the network diagram that these interesting gene groups we discussed are not adjacent. The other genes linking them are worthy of attention and further research.
In WGCNA analysis, the red module (mup6, mup8) is mainly related to lipid metabolism, oxidative stress, and inflammation; the darkslateblue module (cd53, fcer1g, cd68, ctss, laptm5) is mainly related to cell activation, innate immune system, and atherosclerosis; the darkmagenta module (sub1, snrpd2, zfp931, etohi1) is mainly related to transcriptional regulation, mRNA splicing, and gene expression; the antiquewhite4 module (egln1) is mainly related to cellular oxygen sensor that catalyzes, under normoxic conditions (Fig. 6). More specifically, the cd68 was a surface marker for M1 macrophages and it was involved in liver damage such as inflammation, liver fibrosis, and HCC (Liu et al. 2019). Shi et al. reported that PM 2.5 upregulated the expression of CD68 both in cell model and in lung tissues (Shi et al. 2019), indicating PM 2.5 promoted the pro-inflammatory transformation of macrophage thus inducing tissue damage. The included cd68 and cd53 are related to liver inflammation and insulin sensitivity (Ehses et al. 2009). Inhibition of the expression of the cd family may be a therapeutic target for HCC. Additionally, PM 2.5 exposure could increase the glucocorticoids in plasma by reducing the expression of glucocorticoid receptors in the hippocampus, thereby activating the inflammatory response and inducing neurotoxicity (Jia et al. 2018) cd68, ctss, laptm5, fcgr3a, and cd53 were related to the regulation of microglia polarization and can detect out neuropathic pain early (Yu et al. 2020). On the other hand, PM 2.5 has been confirmed in the population to cause neurodegenerative diseases such as Parkinson's syndrome, even if the concentration is lower than the current American national standard (Liu et al. 2016); animal experiments have shown that PM 2.5 may aggravates Parkinson's disease via inhibition of autophagy and mitophagy pathway (Wang et al. 2021). It has also showed that ctss, cd53, igsf6, ptprc, and laptm5 may be potential pathological target gene for the Parkinson's syndrome, which is highly similar to our darkslateblue module (Cui et al. 2015). We can infer that darkslateblue module can be used as biomarker for neurodegenerative diseases such as Parkinson's disease. It is important that the genes in darkslateblue module increase after PM 2.5 exposure, but decrease with metformin. These potential biomarkers might be helpful for the prediction and early screening of these related diseases.

Conclusion
This study investigated the disturbance of transcriptome level in ob/ob mice liver induced by concentrated PM 2.5 exposure through small animal whole-body dynamic exposure system, and meanwhile preliminarily explored the effects of metformin intervention in this process. The results showed that PM 2.5 could affect thyroid function, insulin resistance, glucose, and lipid metabolism in obese fatty liver, which may be related to the mechanism of PM 2.5 -induced liver diseases such as NAFLD and HCC. The screened 12 hub genes might be used as potential biomarkers for air pollution health risk assessment and the expression of 5 genes in darkslateblue module (cd53, fcer1g, cd68, ctss, laptm5) increased after PM 2.5 exposure and decreased after metformin intervention, which could provide clues for the related mechanism and the protective effect of metformin in the detrimental effects in the obesity fatty liver caused by PM 2.5 . But further research is still needed to explore the unequivocal mechanism involved in above-mentioned process.
Abbreviations KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: Gene Onotology; STC-GO: Series test of cluster of gene ontology; DEGs: Differentially expressed genes; WGCNA: Weighted correlation network analysis; qRT-PCR: Real-time quantitative polymerase chain reaction Author contribution LL: writing-first draft, software, visualization, and investigation. LT: writing-first draft, conceptualization, and methodology. TL: visualization and investigation. MS: data curation and investigation. JD: investigation and writing-reviewing and editing. YY: supervision, writing-reviewing and editing. ZS: writingreviewing and editing.
Funding This work was supported by the National Natural Science Foundation of China (91943301, 92043301, 81602876). The authors would like to thank Weiping Tang of Cnkingbio biotechnology Co. Ltd. for bioinformatics assistance.

Declarations
Ethics approval This work has received approval for research ethics from the Animal Care and Use Committee of Capital Medical University, which ethical approval number is AEEI-2019-161.

Consent to participate Not applicable.
Consent for publication All listed authors have approved the manuscript before submission, including the names and order of authors.

Competing interests
The authors declare no competing interests.