MX1 and UBE2L6 are potential metaflammation gene targets in both diabetes and atherosclerosis

Animal models and interventions

The animal study protocol was approved by the Ethics Committee of the Third Medical Center of Chinese PLA General Hospital (KY2021-008; 2021/02/02). Pigs have been widely used as models in studies on metabolism and cardiovascular disease because they share similarities with humans in anatomy, physiology, and genetics (Zhao et al., 2018). As described previously (Zhao et al., 2022), 24 genetically wild-type male Ba-Ma mini pigs (5–6 months, 12.6–15.2 kg) were purchased from Beijing Strong Century Minipigs Breeding Base and randomly divided into five groups. Pigs in the STZ+HC and STZ+HL groups were injected with STZ first and then fed the HC or HL diet for 9 months. In contrast, pigs in the HC+STZ and HL+STZ groups were fed the HC or HL diet for 9 months and injected with STZ at 3 months. The control group was only fed a regular diet for 9 months. The proportions of the main components of tallow/cholesterol/bile salt/basic feed were 15/0/0.5/84.5% and 15/2/0.5/82.5% in the HC and HL experimental diets, respectively. The animals were housed in separate cages (20–26 °C, relative humidity of 40–70%, and a light/dark cycle of 12/12 h). After deletion of the missing data, 20 cases were included in the study: control group (n = 4), HC+STZ group (n = 4), STZ+HC group (n = 3), HL+STZ group (n = 4), and STZ+HL group (n = 5). Anaesthesia for all animals was induced by intramuscular injection of xylazine hydrochloride (25 mg/kg) and midazolam (0.25 mg/kg) and maintained by intravenous injection of pentobarbital sodium (6 mg/kg), after which the femoral artery was bled until the animals died. The blood glucose levels and the size of abdominal aorta atherosclerotic plaque are the basis to judge the success of DM and AS modelling. The plaques and their sizes were observed and estimated through oil red O staining. There were no surviving animals after the experiment.

Transcriptome sequencing and data preparation

Venous blood from pigs was placed in a tube containing EDTA, red blood cell lysate was added, and the sediment was taken after repeated blowing and centrifugation. We extracted RNA as described in a previous report (Hua et al., 2018) and assessed its integrity using the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) with the RNA Nano 6000 Assay Kit as previously described (Nawaz et al., 2018). As input material for RNA sample preparation, 1 µg of RNA was used per sample. According to the manufacturer’s instructions, sequencing libraries were generated using the NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA). For each sample, index codes were added to the attribute sequences. PCR products were purified using an AMPure XP system, and the library quality was assessed using the Agilent Bioanalyzer 2100 system. Next, the index-coded samples were clustered using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA) on a cBot Cluster Generation System. A 150-bp pair-end read was generated from the libraries prepared on an Illumina NovaSeq platform after cluster generation. The biological replicate numbers of the control, HC+STZ, STZ+HC, HL+STZ, and STZ+HL groups were 4, 4, 3, 4, and 5, respectively. The technical replicate number of all samples was 2.

The statistical power of this experimental design, calculated in RNASeqPower, is 0.87. We processed in-house Perl scripts to process quality control raw data in fastq format. Clean data were obtained by removing adapters, poly-N sequences, and low-quality reads from raw data, which were the basis of all downstream analyses. Additionally, the Q20, Q30, and GC contents of the clean data were calculated.

Reads mapping to the reference genome and quantification of gene expression levels

From the ensembl genomes website, we downloaded the reference genome (https://ftp.ensembl.org/pub/release-98/fasta/sus_scrofa_largewhite/dna/) and annotation files (https://ftp.ensembl.org/pub/release-98/gtf/sus_scrofa_largewhite/) for the Sus scrofa. HISAT2 v2.0.5 was used to build the index of the reference genome, and paired-end clean reads were aligned with it. The read numbers mapped to each gene were counted using FeatureCounts v1.5.0-p3. Based on the gene length and read counts mapped to this gene, the expected number of fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) of each gene was calculated.

Identification of differentially expressed genes in DM and AS

The gene expression dataset was normalized by log₂ transformation with Benjamini–Hochberg correction to control the false discovery rate. A p value < 0.05 and absolute log₂ fold change (FC) ≥ 1 were regarded as threshold criteria for significant differentially expressed genes (DEGs) of interest. The DESeq2 R package (1.16.1) was used in this study to analyse differential expression between any two groups based on a negative binomial model. Benjamini and Hochberg’s approach to controlling false discovery rates was used to adjust the resulting p values. In DESeq2, DEGs with an adjusted p value of 0.05 were considered.

Gene ontology and KEGG pathway enrichment analyses

In the clusterProfiler R package, the gene length bias was corrected, and a corrected p value of 0.05 was considered significant enrichment for DEGs. The statistical enrichment of DEGs was also tested among the KEGG (http://www.genome.jp/kegg/) pathways. An adjusted p value < 0.05 was considered significant for all enrichment analyses.

Protein–protein interaction analysis

The STRING database (Szklarczyk et al., 2021) was used to construct protein–protein interaction (PPI) networks of the proteins encoded by the identified DEGs. Network analysis and hub protein identification were performed using the cytoHubba modules in Cytoscape 3.7.0 (Shannon et al., 2003).

Statistical analysis

The statistical software, statistical tests and significance thresholds applied for DEGs, Gene Ontology (GO) and KEGG analysis have been described above separately. Correlation analysis of blood glucose levels and gene expression was performed. Comparisons of multiple groups of measurement data were performed with one-way ANOVA, and the subsequent pairwise comparison adopted the Bonferroni method. A p value < 0.05 was considered significant for the analysis.

Identification of DEGs common to HC+STZ vs control, STZ+HC vs control, HL+STZ vs control, and STZ+HL vs control

The blood glucose levels of pigs after modelling were 3.20 ± 0.91, 23.23 ± 3.28, 2.93 ± 0.92, 24.50 ± 2.68, and 3.80 ± 0.71 mmol/L in the control, HC+STZ, STZ+HC, HL+STZ, and STZ+HL groups, respectively. Abdominal aorta atherosclerotic plaques were found in all pigs from the HC+STZ group. In other models, only three cases from the STZ+HL group showed plaques (n = 1) or lipid streaks (n = 2). The significant DEGs from Ba-Ma mini pigs were identified using applied combinatorial statistical methods. A total of 103 DEGs were identified as common between the HC+STZ and STZ+HC groups, which were different from the HL+STZ and STZ+HL groups (Fig. 1). The numbers of consistently upregulated and downregulated genes in both groups were 44 and 42, respectively. Moreover, 17 genes were upregulated in the HC+STZ group and downregulated in the STZ+HC group. As shown in Fig. 2, the distribution of 103 DEGs was different not only between the groups but also within the groups, suggesting the potential application value of DEGs in distinguishing and refining the molecular mechanism of DM and AS. According to the hyperglycaemia and plaque results, it is speculated that these 103 DEGs may be involved in the common mechanism of DM and AS, and these genes were included in the further follow-up analysis.

Significant GO terms and pathways

The enrichment analysis of identified genes was used to assess the biological processes (BPs), cellular components (CCs), and molecular functions (MFs) in which the enriched genes participated. The top three enriched BPs of DEGs in HC+STZ vs control or STZ+HC vs control were regulation of catabolic process, cellular catabolic process, and organic substance catabolic process (p < 0.05, Fig. 3A) or antigen processing/presentation, carbohydrate derivative metabolic process, and immune response (p < 0.05, Fig. 4A), respectively. The top three significant CCs of DEGs in HC+STZ vs control or STZ+HC vs control were proteasome core complex, proteasome complex, and endopeptidase complex (p < 0.05, Fig. 3B) or microtubule organizing centre, preribosome, and chromosome (p < 0.05, Fig. 4B), respectively. Carboxy-lyase activity, purine nucleoside binding, and GTP binding (p < 0.05, Fig. 3C) or phosphoric diester hydrolase activity, sulfuric ester hydrolase activity, and hydrolase activity (p < 0.05, Fig. 4C) were identified as significant MFs. The most significantly enriched pathways analysed through KEGG were endocytosis, regulation of actin cytoskeleton, fructose, mannose metabolism (p < 0.05, Fig. 3D) and DNA replication and the toll-like receptor signalling pathway (p < 0.05, Fig. 4D). Collectively, these data suggest that the DEGs in the HC+STZ vs control and STZ+HC vs control groups were mainly enriched in catabolic pathways, metabolic processes, immune responses and DNA replication.

The top three enriched BPs were regulation of response to biotic stimulus, immune response, and immune system process (p < 0.05, Fig. 5A). The top three significant CCs of the 103 DEGs common to HC+STZ and STZ+HC were extracellular region, proteinaceous extracellular matrix, and viral capsid (p < 0.05, Fig. 5B). Chemokine activity, chemokine receptor binding, G-protein coupled receptor binding and others were identified as significant MFs (p < 0.05, Fig. 5C). The most significantly enriched pathways were influenza A, hepatitis C, and measles (p < 0.05, Fig. 5D). These results demonstrated that 103 common DEGs were mainly enriched in immune response and regulation.

Identification of hub proteins

The PPI network was first constructed by retrieving the global and internal interactions of the 103 common DEGs between the HC+STZ and STZ+HC groups from the STRING database (Fig. 6A). The PPI network consisted of 648 nodes (72 nodes from the common DEGs) and 2,106 edges, revealing the top 10 hub proteins, namely, ISG15, IRG6, IRF7, IFIT3, MX1, UBE2L6, DDX58, IFIT2, USP18, and IFI44L (Fig. 6B). The interactions within 103 common DEGs (14 nodes and 31 edges) were further displayed using the PPI network shown in Fig. 6C. Unexpectedly, MX1 and UBE2L6 were the intersection of internal and global PPI networks of 103 common DEGs, suggesting their key role in the comorbidity mechanism of DM and AS (Figs. 6A and 6C). Our results revealed that MX1 and UBE2L6 have different roles, in addition to the traditional anti-virus activity and ubiquitination degradation of proteins; furthermore, the data indicate that the above biological processes may be also involved in the pathophysiology of DM and AS.

Differential expression of MX1 and UBE2L6 and their contributions to DM and AS

The expression of MX1 and UBE2L6 was 507.22 ± 342.56 and 96.99 ± 49.92 in the HC+STZ group, respectively, which was significantly higher than that in the control group (p = 0.0178 and p = 0.0188, respectively; Figs. 7A and 7B). In contrast, the expression of MX1 and UBE2L6 in the other group with stable hyperglycaemia was only 31.10 ± 19.21 and 28.21 ± 14.30, similar to that in the control group (p < 0.05 and p < 0.05, respectively; Figs. 7A and 7B). The Sankey diagram traced the glucose levels and gene expression of models with or without AS. Obviously, in the HC+STZ group, there were two cases with atherosclerotic plaques exceeding 50% of the lumen and two other cases with plaques that were less than half of the lumen, with serious hyperglycaemia. The coexistence of DM and AS was related to the overexpression of MX1 and UBE2L6 in the HC+STZ group (Figs. 7C and 7D).

MX1- and UBE2L6-related PPI networks for calcium and micronutrients

A total of 13,802 and 382 potential candidate genes were identified in the GeneCards databases. Calcium- and micronutrient-related PPI networks were constructed by retrieving the interactions of the top 1,998 candidates (if applicable) with MX1 and UBE2L6 from the STRING database. The calcium-related PPI network consisted of 1,744 nodes and 33,764 edges. The micronutrient-related PPI network consisted of 342 nodes and 2,825 edges. The calcium-related PPI network in which MX1 and UBE2L6 participate contained 58 nodes and 424 edges (Fig. 8A). The micronutrient-related PPI network in which MX1 and UBE2L6 participate contained 18 nodes and 58 edges (Fig. 8B). The 13 common proteins were MX1, UBE2L6, IFNG, CD4, TNF, VEGFA, CXCL10, PRDX5, F2, IL10, CCL2, IL6, and TP53. The most significantly enriched KEGG pathways were glutathione metabolism, pyrimidine metabolism, purine metabolism, and metabolic pathways. The data revealed that MX1- and UBE2L6 have a close relation to calcium and micronutrient metabolic pathways.