Transcriptomics of type 2 diabetic and healthy human neutrophils.

OBJECTIVES: Chronic inflammatory diseases, including diabetes and cardiovascular disease, are heterogeneous and often co-morbid, with increasing global prevalence. Uncontrolled type 2 diabetes (T2D) can result in severe inflammatory complications. As neutrophils are essential to normal and aberrant inflammation, we conducted RNA-seq transcriptomic analyses to investigate the association between neutrophil gene expression and T2D phenotype. As specialized pro-resolving lipid mediators (SPM) act to resolve inflammation, we further surveyed the impact of neutrophil receptor binding SPM resolvin E1 (RvE1) on isolated diabetic and healthy neutrophils.
METHODS: Cell isolation and RNA-seq analysis of neutrophils from N = 11 T2D and N = 7 healthy individuals with available clinical data was conducted. Additionally, cultured neutrophils (N = 3 T2D, N = 3 healthy) were perturbed with increasing RvE1 doses (0 nM, 1 nM, 10 nM, or 100 nM) prior to RNA-seq. Data was evaluated through a bioinformatics pipeline including pathway analysis and post hoc false discovery rate (FDR)-correction.
RESULTS: We observed significant differential expression of 50 genes between T2D and healthy neutrophils (p < 0.05), including decreased T2D gene expression in inflammatory- and lipid-related genes SLC9A4, NECTIN2, and PLPP3 (p < 0.003). RvE1 treatment induced dose-dependent differential gene expression (uncorrected p < 0.05) across groups, including 59 healthy and 216 T2D neutrophil genes. Comparing T2D to healthy neutrophils, 1097 genes were differentially expressed across RvE1 doses, including two significant genes, LILRB5 and AKR1C1, involved in inflammation (p < 0.05).
CONCLUSIONS: The neutrophil transcriptomic database revealed novel chronic inflammatory- and lipid-related genes that were differentially expressed between T2D cells when compared to controls, and cells responded to RvE1 dose-dependently by gene expression changes. Unraveling the mechanisms regulating abnormalities in diabetic neutrophil responses could lead to better diagnostics and therapeutics targeting inflammation and inflammation resolution.

Background

The increasing global prevalence of chronic inflammatory diseases, such as diabetes, is of critical concern to human health. The global prevalence of diagnosed adult diabetes was ~ 415 million in 2015, and is estimated to rise to ~ 642 million by 2040 [1]. Type 2 diabetes (T2D) is modulated by defective metabolic and exacerbated immune responses that can lead to chronic inflammatory complications, including kidney, nerve, cardiovascular, eye, and periodontal diseases [2], with up to 60% of diabetic individuals having moderate to severe periodontitis [3]. Nevertheless, the endogenous pathways, genes, and cells that trigger and sustain unresolved low-grade inflammation in T2D are not completely understood. While acute inflammation is part of a protective host response, chronic diseases result from uncontrolled inflammation and failure of immune cells to restore homeostasis. In particular, the critical role of the most abundant immune cells, neutrophils, in the initiation and termination (resolution) of inflammation has become apparent. Elucidating how neutrophil abnormalities are a common link across chronic inflammatory diseases, including diabetes [4-8], can lead to better understanding of the disease mechanistically. We now know that such neutrophil abnormalities in chronic, unresolved inflammation include impaired neutrophil adhesion, migration, chemotaxis, cytokine signaling, and phagocytosis, as well as increased neutrophil infiltration, oxidative stress and degranulation [4, 5, 9]. Inflammatory impairment in diabetes has shown improvements after clinical treatments, demonstrating that cells are not permanently damaged. This is consistent with reduction of glucose therapy and periodontal inflammation [10]; yet no specific markers or treatments are available to target inflammatory cells. Appropriate resolution of inflammation programs aim to restore tissue homeostasis following acute inflammation and neutrophil migration [5, 11]. This active process involves key endogenous lipid ligand mediators, such as the specialized pro-resolving mediator resolvin E1 (RvE1), which can bind and transduce agonist signals through innate immune receptors on neutrophils to resolve inflammation [4, 12].

Despite considerable exploration of T2D genetics, primarily through genome-wide association studies (GWAS), and the identification of over 200 implicated genomic regions, to date there has been limited ability to translate genetic information into clinically actionable subtypes of this highly heterogeneous condition [13, 14]. Further, identified genomic variants tend to have very small effect sizes and these genetic studies have been heavily biased by inclusion of primarily European populations. While genetics unequivocally plays an important role in T2D, with heritability estimates ranging from a low of 25% to a large twin cohort that estimated 72% heritability [15], the role of gene expression has been less studied, particularly in the context of disease-relevant cell types. Indeed, as T2D is associated with both increased adipose tissue and circulating inflammatory mediators, which induce insulin resistance and develop a state of chronic and unresolved inflammation [16], a global survey of transcriptomics provides an unbiased method to investigate functional mechanisms of T2D. Non-targeted, unbiased transcriptomics thus has the potential to elucidate novel mechanisms important for disease biomarkers and drug efficacy, especially by revealing unknown pathways of myeloid immune cells and providing dynamic molecular information important to understanding neutrophil heterogeneity in disease.

As neutrophils are essential to both initiation and resolution of inflammation and their gene expression in diabetes had not been previously studied, we conducted RNA sequencing (RNA-seq) transcriptomic analyses to investigate gene expression in human neutrophils from T2D and healthy subjects as a baseline investigation. In this study, we found that T2D individuals had higher blood glucose and more advanced periodontal disease compared to healthy individuals, consistent with worsening overall systemic health. When we investigated global neutrophil gene expression, we observed significant downregulation of several inflammatory- or lipid metabolism-related genes in T2D neutrophils. We further conducted a perturbation study of ex vivo cells to investigate early neutrophil markers, as well as the impact of exogenous RvE1 treatment on T2D and healthy neutrophil gene expression and the inflammatory cytokine signaling response. Our results demonstrate that diabetic neutrophils were not permanently damaged and that cell plasticity was modulated by RvE1. In fact, distinct trends in gene expression and specific cytokine level differences between T2D and healthy neutrophils were modified by RvE1 treatment. To our knowledge, this study represented the first investigation into developing a transcriptome database of diabetic neutrophils, revealing novel signals from type 2 diabetics. By investigating gene patterns and cellular pathways that control aberrant inflammation, novel hypotheses can be generated to improve neutrophil biology and improve precision medicine.

Results

Clinical characteristics of study subjects

Our study population (including both the main transcriptomics and cell culture perturbation analyses) consisted of 13 T2D cases and 8 healthy controls with clinical and demographic information on: age, sex, self-reported ethnicity, self-reported smoking status, body-mass index (BMI), blood total cholesterol, blood glucose, percent hemoglobin A1C (HbA1c), periodontal condition, neutrophil cell counts, and monocyte cell counts (Table 1). The T2D individuals in the cohort had higher blood glucose and more advanced periodontal disease compared to healthy individuals. Diabetes diagnosis was confirmed by elevated blood glucose and HbA1c levels, with blood glucose > 126 mg/dl and HbA1c >  6.5% indicating diabetes, consistent with clinical criteria and the literature [17]; HbA1c was not measured for clinically diagnosed healthy control individuals in this cohort. T2D individuals tended to have higher cholesterol, BMI, neutrophil counts, and monocyte counts, but these factors did not reach statistical significance, while T2D individuals were significantly older than healthy individuals (p = 0.007). Principal component analysis (PCA) was conducted for all demographic and clinical characteristics, and the primary differential clustering was by disease, not demographic or clinical characteristics (Supplementary Fig. S1).

Table 1

Clinical demographics of study subjects

		Healthy (N = 8)	Type 2 Diabetic (N = 13)	Total Cohort (N = 21*)	p-value
Age (mean in years ± SD)		39 (± 9.63)	58 (± 11.05)	51 (± 14.04)	0.007
Sex (no., %)	Male	4 (50%)	7 (53.85%)	11 (52.38%)
	Female	4 (50%)	6 (46.15%)	10 (47.62%)	–
Ethnicity	Caucasian	5 (62.50%)	4 (30.77%)	9 (42.86%)
	Hispanic	0	1 (7.69%)	1 (4.76%)
	African-American	3 (37.50%)	8 (61.54%)	11 (52.38%)	–
Smoking Status (no., %)	Smokers	1 (12.50%)	1 (7.69%)	2 (9.52%)
	Former Smokers	0	2 (15.38%)	2 (9.52%)
	Never Smokers	7 (87.50%)	10 (76.92%)	17 (80.95%)	–
BMI (kg/m2; mean ± SD)		28.56 (± 4.47)	30.46 (± 6.27)	29.74 (± 5.61)	0.47
Blood Cholesterol (mg/dl; mean ± SD)		144.25 (± 20.61)	183.77 (± 66.97)	168.71 (± 56.80)	0.12
Blood Glucose (mg/dl; mean ± SD)	Normal (<  100 mg/dl)	96.50 (± 3.54)	–	96.50 (± 3.54)
	Pre-diabetes (100-125 mg/dl)	107 (± 3.81)	–	107 (± 3.81)
	Diabetes (> 126 mg/dl)	209	267.77 (± 123.38)	263.57 (± 119.58)	0.003
HbA1c % (mean ± SD)	Normal (<  5.7% HbA1c)	–	–	–
	Pre-diabetes (5.7–6.4% HbA1c)	–	6.15 (± 0.17)	6.15 (± 0.17)
	Diabetes (>  6.5% HbA1c)	–	9.64 (± 2.63)	9.64 (± 2.63)	–
Periodontal Condition (no., %)	Healthy	8 (100%)	0	8 (38.09%)
	Mild	0	3 (23.08%)	3 (14.29%)
	Moderate	0	6 (46.15%)	6 (28.57%)
	Severe	0	2 (15.38%)	2 (9.52%)
	Gingivitis	0	2 (15.38%)	2 (9.52%)	–
Neutrophil Count (mil. of cells, mean ± SD)		104.88 (± 57.23)	123.69 (± 65.68)	116.52 (± 61.83)	0.51
Monocyte Count (mil. of cells, mean ± SD)		72.50 (± 34.13)	99.59 (± 73.72)	89.27 (± 62.05)	0.34

* N = 2 healthy and N = 1 diabetic subjects were included in both the serum and cell culture perturbation analyses. p-values calculated by unpaired t-tests (two-sided p-values; italicized p < 0.05 significant). SD standard deviation, BMI body-mass index

Global gene expression levels of neutrophils in health and disease

We observed statistically significant differential expression of 50 genes (FDR-corrected p < 0.05, Fig. 1A, Table 2, Supplementary Fig. S2) between T2D and healthy neutrophils. Differentially expressed genes were primarily inflammatory- or lipid-associated and were almost universally downregulated in T2D individuals relative to healthy individuals, with the exception of GTSCR1 (Gilles de la Tourette syndrome chromosome region, candidate 1), a non-coding RNA gene that was upregulated in T2D neutrophils.

Fig. 1

Neutrophils from type 2 diabetic subjects show down-regulated expression of inflammatory and lipid genes relative to healthy neutrophils. RNA-seq transcript expression of T2D (N = 11) vs. healthy (N = 7) subject neutrophils. A Volcano plot of ln-transformed gene expression (x-axis) and the FDR-adjusted negative log-transformed p-value (y-axis). The dotted lines represent selected significant differentially expressed gene thresholds: 0.6 for fold-change, 0.05 for significance. Blue and red indicate significantly down- and up-regulated genes, respectively. B Bean plots of the ln-transformed relative abundance of gene expression (line at mean) of the top 3 most significant down-regulated genes and the one significant up-regulated gene (GTSCR1) plotted for T2D vs. healthy individuals

Table 2

All statistically significant differentially expressed neutrophil genes comparing T2D (N = 11) vs. healthy (N = 7) subjects

Gene	Gene Name	Annotated Role	Log(FC)	p-value	FDR-corrected p-value
SLC9A4	solute carrier family 9 member A4	pH Regulation / Signaling	−5.42	7.46E-08	0.001
NECTIN2	Nectin cell adhesion molecule 2	Immune	−3.77	2.23E-07	0.002
PLPP3	phospholipid phosphatase 3	Lipids / Inflammation	−5.37	3.38E-07	0.002
C3	complement C3	Inflammation / Immune	−3.61	7.92E-07	0.003
IRAK2	interleukin-1 receptor-associated kinase 2	Inflammation / Immune	−4.01	8.93E-07	0.003
KCNH4	potassium voltage-gated channel subfamily H member 4	Lipids / Insulin	−4.41	9.15E-07	0.003
C1orf61	chromosome 1 open reading frame 61	Transcription	−4.71	2.08E-06	0.005
CCRL2	C-C motif chemokine receptor like 2	Inflammation / Immune	−3.37	2.61E-06	0.005
GNPDA1	glucosamine-6-phosphate deaminase 1	Glucose Metabolism	−2.95	2.70E-06	0.005
CLIC4	chloride intracellular channel 4	Cell Processes	−3.54	3.76E-06	0.007
GPR84	G protein-coupled receptor 84	Inflammation / Immune	−3.41	4.09E-06	0.007
PLAU	plasminogen activator, urokinase	Immune / Blood	−3.83	7.31E-06	0.01
SLC43A3	solute carrier family 43 member 3	Transmembrane Transport	−3.34	1.04E-05	0.01
RHOV	ras homolog family member V	Cell Processes / Signaling	−4.11	1.14E-05	0.01
SGPP2	sphingosine-1-phosphate phosphatase 2	Inflammation / Immune	−3.42	1.19E-05	0.01
TNFAIP6	TNF alpha induced protein 6	Inflammation / Immune	−3.01	1.69E-05	0.02
C15orf48	chromosome 15 open reading frame 48	Immune	−4.01	1.84E-05	0.02
TLCD2	TLC domain containing 2	Lipids	−1.60	2.28E-05	0.02
GTSCR1	Gilles de la Tourette syndrome chromosome region, candidate 1	-*	3.71	2.39E-05	0.02
SLC25A13	solute carrier family 25 member 13	Mitochondrial / Cell Processes	−2.32	2.57E-05	0.02
RP11-1193F23.1	–	–	−2.68	2.68E-05	0.02
SPHK1	sphingosine kinase 1	Lipids / Inflammation	−3.33	2.73E-05	0.02
DNAH17	dynein axonemal heavy chain 17	Motility (Sperm)	−2.94	2.79E-05	0.02
SPATC1	spermatogenesis and centriole associated 1	–	−3.09	2.96E-05	0.02
RALGDS	ral guanine nucleotide dissociation stimulator	Cell Processes	−2.69	3.01E-05	0.02
VEGFA	vascular endothelial growth factor A	Cell Growth / Angiogenesis / Lipids	−2.01	3.05E-05	0.02
HTRA3	HtrA serine peptidase 3	Signaling / Binding	−3.81	3.37E-05	0.02
CTB-114C7.4	–	–	−2.36	3.48E-05	0.02
MFSD2A	major facilitator superfamily domain containing 2A	Lipids / Transmembrane Transport (Brain)	−2.96	3.88E-05	0.02
AGAP3	ArfGAP with GTPase domain, ankyrin repeat and PH domain 3	Signaling (Brain)	−2.61	3.90E-05	0.02
PNPLA1	patatin like phospholipase domain containing 1	Lipids	−2.53	4.15E-05	0.02
ABCG1	ATP binding cassette subfamily G member 1	Inflammation / Immune / Lipids	−1.45	4.28E-05	0.02
PTGES	prostaglandin E synthase	Inflammation / Lipids	−3.13	4.46E-05	0.02
TBC1D30	TBC1 domain family member 30	Transport / GTPase	−2.14	4.55E-05	0.02
GPR4	G protein-coupled receptor 4	Angiogenesis / Lipids	−3.48	5.35E-05	0.03
CD300LD	CD300 molecule like family member d	Immune	−2.41	5.97E-05	0.03
KLHDC8B	kelch domain containing 8B	Cell Processes / Mitosis	−2.25	7.90E-05	0.04
ANKRD33B	ankyrin repeat domain 33B	–	−2.05	8.21E-05	0.04
LPCAT1	lysophosphatidylcholine acyltransferase 1	Lipids	−1.86	8.31E-05	0.04
RP1-239B22.5	–	–	−3.87	9.27E-05	0.04
CTA-126B4.7	uncharacterized LOC101927344	–	−2.70	1.06E-04	0.05
TRAF3	TNF receptor associated factor 3	Immune	−2.01	1.22E-04	0.05
LAMB3	laminin subunit beta 3	Cell Processes	−3.02	1.23E-04	0.05
TTPAL	alpha tocopherol transfer protein like	–	−2.13	1.25E-04	0.05
TMEM231	transmembrane protein 231	Ciliogenesis	−2.28	1.28E-04	0.05
AZIN1	antizyme inhibitor 1	Cell Growth / Homeostasis	−1.90	1.30E-04	0.05
MARCH3	membrane associated ring-CH-type finger 3	Endosomal Transport	−2.31	1.38E-04	0.05
LIMS1	LIM zinc finger domain containing 1	Cytoskeleton / Signaling	−2.14	1.38E-04	0.05
NCR3LG1	natural killer cell cytotoxicity receptor 3 ligand 1	Immune	−3.55	1.41E-04	0.05
SQSTM1	sequestosome 1	Immune / Signaling	−1.78	1.42E-04	0.05

* - indicates no information available. FDR-corrected p-value < 0.05 (before rounding) considered significant. FC fold-change; FDR false discovery rate

The top 3 differentially expressed genes were downregulated in T2D neutrophils (Fig. 1B): inflammatory signaling gene SLC9A4 (solute carrier family 9 member A4; log(FC) = − 5.42, FDR-corrected p = 0.001), immune regulating gene NECTIN2 (nectin cell adhesion molecule 2; log(FC) = − 3.77, FDR-corrected p = 0.002), and anti-inflammatory gene PLPP3 (phospholipid phosphatase 3; log(FC) = − 5.37, FDR-corrected p = 0.002); while only GTSCR1 was upregulated in T2D neutrophils (log(FC) = 3.71, FDR-corrected p = 0.02, Fig. 1B). GTSCR1 has been identified as a regulatory gene in cardiac inflammation via human CD4+ T cell transcriptomics cultured in Th17-promoting conditions [18].

Of the significant differentially expressed neutrophil genes, 34% (17 / 50) had known immune or inflammation related roles, with 24% (12 / 50) of genes linked to lipid or glucose metabolism (Table 2). In addition to NECTIN2 and PLPP3, 15 other annotated immune or inflammation associated genes had significantly lower gene expression in T2D neutrophils (Table 2), while KCNH4, GNPDA1, TLCD2, VEGFA, MFSD2A, PNPLA1, GPR4, and LPCAT1 showed lower expression in T2D neutrophils and had previously been linked to lipid or glucose metabolism. Indeed, several significant differentially expressed genes, including PLPP3, SPHK1, ABCG1, and PTGES, have previously been indicated to have both lipid- and inflammatory-related roles.

This interaction between lipid and inflammatory pathways was further evident when we investigated neutrophil gene expression by biological pathways and observed that most of the significant biological pathways were lipid-related, despite the abundance of immune- and inflammatory-related genes among our significant results (Table 3, Supplementary Figs. S3-S4).

Table 3

Top over-represented pathways included in pathway analyses

Pathway	p-value	p-value (FDR)	p-value (Bonferroni)
KEGG Pathways
Sphingolipid metabolism	2.52E-04	0.02	0.02
Ether lipid metabolism	0.004	0.14	0.28
Phospholipase D signaling pathway	0.02	0.24	1
Fc gamma R-mediated phagocytosis	0.02	0.24	1
Glycerophospholipid metabolism	0.02	0.24	1
GO Biological Processes
Ammonium ion metabolic process	9.30e-5	0.10	0.15
Surfactant homeostasis	1.20e-4	0.10	0.20
GO Molecular Functions
Sphingosine-1-phosphate-phosphatase activity	1.20e-4	0.03	0.03
GO Cellular Components
Plasma membrane	2.00e-4	0.03	0.04
Integral component of plasma membrane	4.80e-4	0.03	0.11

Uncorrected p-value < 0.05 considered over-represented; Bonferroni or FDR-corrected p-value < 0.05 considered significant. FDR false discovery rate; KEGG Kyoto Encyclopedia of Genes and Genomes; GO Gene Ontology

Overall, 13 pathways and 19 diseases were over-represented (uncorrected p < 0.05) in our cohort by gene expression differences. The top 5 over-represented KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways were heavily lipid-related: sphingolipid metabolism, ether lipid metabolism, phospholipase D signaling pathway, Fc gamma R-mediated phagocytosis, and glycerophospholipid metabolism; though only sphingolipid metabolism was statistically significant after Bonferroni or FDR correction (corrected p = 0.02, Supplementary Table S1). Of the 19 over-represented diseases, 15 remained significant by FDR correction, but none remained significant after Bonferroni correction, and all 19 had only a single differentially expressed gene represented in that particular disease pathway. The top 5 disease pathways were: citrullinemia, Paget’s disease, endometriosis, junctional epidermolysis bullosa, and avascular necrosis/osteonecrosis of femoral head (Supplementary Table S2). Additionally, top Gene Ontology (GO) biological process, molecular function, and cellular component pathways were interrogated [19, 20]. Following FDR correction, the top GO pathways also included lipid-related and plasma membrane pathways (Table 3).

When we investigated the top KEGG and GO biological pathways (Table 3) by individual gene expression across the top 50 genes, we observed hierarchical clustering of pathway and gene expression by groups of individuals (Fig. 2). These clusters primarily separated out diabetic and healthy individuals and included disease heterogeneity in gene levels.

Fig. 2

Neutrophil gene expression in top biological pathways is primarily clustered by diabetic vs. health status. RNA-seq transcript expression of T2D (D; N = 11) vs. healthy (H; N = 7) subject neutrophils. Heatmap of T2D and healthy gene expression by over-represented KEGG and Gene Ontology (GO) pathways and top 50 genes (FDR-corrected p < 0.05) in those pathways. Blue indicates decreased and red indicates increased gene expression

Neutrophil gene expression by lipid ligand Resolvin E1 dose-response

Because endogenous inflammatory pathways are deficient in chronic inflammatory diseases like T2D, we aimed to understand whether treating isolated neutrophils with an exogenous lipid ligand resolution mediator, RvE1, would impact transcription phenotypes.

We conducted an ex vivo perturbation experiment across a range of clinically relevant RvE1 doses (0–100 nM) and observed that, as expected, dosing of the RvE1 ligand impacted T2D neutrophils differently than healthy neutrophils. We investigated this dose-response phenotype between three comparisons: 1) T2D neutrophils only, 2) healthy neutrophils only, and 3) T2D vs. healthy neutrophils.

RvE1 treatment induced differential gene expression (uncorrected p < 0.05) across doses (Fig. 3), including 59 genes in healthy and 216 genes in T2D neutrophils. Comparing T2D to healthy neutrophils, 1097 genes were differentially expressed across treatment doses, including two statistically significant (FDR-corrected p < 0.05) inflammatory genes, LILRB5 (leukocyte immunoglobulin like receptor B5) and AKR1C1 (aldo-keto reductase family 1 member C1). Interestingly, the gene CD177, a beta-2 integrin associated heterodimer with ITGAM/CD11b, was also significantly increased in neutrophils from diabetics (p < 0.002). While this gene did not maintain its significance after FDR correction, this is a key regulator of TNF-α primed neutrophils through superoxide production regulation. Neutrophils from T2D subjects had a stronger response to RvE1 treatment dose-dependently (more differentially expressed genes), particularly at the clinically relevant dose, 100 nM [12] (117 differentially expressed genes between 100 nM and 10 nM RvE1), showing modified gene expression by RvE1 treatment. In contrast, healthy neutrophils were not as perturbed by RvE1 dosing (few differentially expressed genes between doses), though some effect was observed at the 10 nM dose (25 differentially expressed genes between 10 nM and 0 nM RvE1). Treatment of T2D neutrophils with 100 nM RvE1 resulted in a reduction in the number of differentially expressed genes in T2D neutrophils compared to healthy neutrophils (98 genes when T2D neutrophils were treated with 100 nM RvE1 versus 169 genes without RvE1 treatment, 0 nM).

Fig. 3

Resolvin E1 treatment has a stronger impact on gene expression in T2D relative to healthy neutrophils. Number of genes differentially expressed (uncorrected p < 0.05) by resolvin E1 (RvE1) treatment dose (0 nM, 1 nM, 10 nM, or 100 nM) across groups: A healthy only, B T2D only, C T2D (higher dose) vs. healthy (lower dose) in cultured neutrophils (N = 3 T2D, N = 3 healthy). D Number of genes differentially expressed (uncorrected p < 0.05) in healthy, T2D and T2D vs. healthy groups, comparing low (1 nM vs. 0 nM) to high (100 nM vs. 0 nM) RvE1 dose. E Bean plots of the relative abundance of gene expression (line at mean) of the significant differentially expressed genes (LILRB5 and AKR1C1; FDR-corrected p < 0.05) and top genes (RN7SL2, LILRB2, and NECTIN2). Samples were plotted for T2D vs. healthy individuals. T2D, type 2 diabetes

An interesting observation was that, in general, the differentially expressed neutrophil genes across RvE1 doses tended to be unique to each comparison: 1) T2D, 2) healthy, or 3) T2D vs. healthy groups. Within each comparison the genes differentially expressed remained mostly consistent across RvE1 doses. However, the genes differentially expressed between T2D and healthy neutrophils in the cell culture perturbation model without RvE1 treatment were generally not the same as the significant differentially expressed serum neutrophil genes observed in the main analysis, with the exception of NECTIN2, HTRA3, and ABCG1, which were differentially expressed in both serum and cell culture. Interestingly, when we perturbed T2D neutrophils with 100 nM RvE1 and did not perturb healthy neutrophils, only HTRA3 remained strongly differentially expressed (p < 0.05), while NECTIN2 and ABCG1 showed less differential expression. Overall, there were distinct trends in neutrophil responses to RvE1 perturbation between diabetic and healthy neutrophils in cell culture (Fig. 3).

Diabetic and healthy neutrophil cytokines respond differently to Resolvin E1 perturbation

To confirm that gene expression changes influenced neutrophil protein expression, we analyzed cytokine production of neutrophil cell culture supernatants following the same RvE1 dose-response studies. Results showed trends of differential inflammatory cytokine concentrations between T2D and healthy neutrophils and across RvE1 doses. Out of a panel of 20 human cytokines, the concentration of 8 cytokines (MIP-1α, IL-4, IL-8, MIP-1β, P-Selectin, sICAM-1, TNF-α, and IL-1α) were within the range of quantification and were included in the analysis (Fig. 4, Supplementary Tables S3-S4). Without any RvE1 perturbation there were differences in levels of these cytokines between T2D and healthy cell cultured neutrophils, with T2D neutrophils having higher TNF-α and P-Selectin levels but lower MIP-1β, IL-8, and sICAM-1 levels. Across all RvE1 doses, P-Selectin levels remained consistently higher, while IL-8 levels remained lower in T2D compared to healthy neutrophils. However, at 1 nM RvE1 treatment, MIP-1β levels in T2D neutrophils rose to the baseline levels of healthy neutrophils, while healthy neutrophil levels stayed consistent. At higher RvE1 doses, including 10 nM and 100 nM, MIP-1β levels rose in both T2D and healthy neutrophils; there was a stronger, dose-dependent increase among healthy neutrophils. After treatment with 10 nM or 100 nM RvE1, TNF-α levels were elevated in healthy neutrophils, reaching levels seen in T2D neutrophils; a similar effect was observed at 100 nM RvE1 for P-Selectin, though healthy levels (1511 pg/mL) did not fully rise to match T2D levels (1751 pg/mL). We also observed some distinct neutrophil cytokine profiles in individual subjects, consistent with known inter-individual cytokine variation. The cytokine level differences did not reach statistical significance at this sample size. Overall, we observed trends of several potentially interesting differences in cytokines between T2D and healthy neutrophils, indicating that both diabetic and healthy neutrophils may have a distinct functional signaling response to RvE1 treatment in a dose-dependent manner.

Fig. 4

Type 2 diabetic and healthy neutrophils showed differential inflammatory cytokine concentrations across resolvin E1 doses. Human inflammatory cytokine levels (cytokines: MIP-1α, IL-4, IL-8, MIP-1β, P-Selectin, sICAM-1, TNF-α, IL-1α; average pg/mL) in response to resolvin E1 (RvE1) dose, comparing T2D (N = 3) to healthy (N = 3) cultured neutrophils with: A vehicle only (0 nM), B 1 nM, C 10 nM, or D 100 nM RvE1 treatment. For each cytokine, values at or below the lower limit of quantification (LLOQ; based on a 7-standard curve at 1:4 dilution) were reported as the LLOQ value

Discussion

We conducted the first study of viable neutrophil gene expression in the chronic inflammatory disease type 2 diabetes. We observed differential expression between T2D and healthy neutrophils of several inflammatory- or lipid-associated genes, indicating the potential for distinct dysbiotic transcriptomic profiles between these groups among a key and novel immune cell type. Further, as many of the significant differentially expressed genes had direct lipid- or glucogenesis-related roles, neutrophil transcriptomics may provide a window directly into impairment of the host response among diabetics and disease pathogenesis beyond their role in immune or inflammatory genes.

The top genes significantly down-regulated in T2D neutrophils compared to healthy neutrophils tended to be biologically relevant across inflammatory- and lipid-associated gene pathways and included SLC9A4, NECTIN2, and PLPP3 (Fig. 1).

SLC9A4 is a plasma membrane solute carrier protein that acts in homeostatic pH regulation and proton transport, including to eliminate metabolism-generated acids and absorb sodium, with an important role in signal transduction [21, 22]. SLC9A4 is in a genomic region with genes involved in the IL-1 receptor (IL-1R) and IL-18 pathways [23, 24]. Ground-breaking studies have elucidated the function and structure of a network of proteins (termed the “inflammasome”) responsible for processing and releasing key pro-inflammatory cytokines (e.g. IL-1β and IL-18) and controlling cell death and resolution of inflammation [25, 26].

This SLC9A4/IL-1R/IL-18 locus has important roles in cytokine signaling and inflammatory response [23, 24]. The IL-18 pathway induces synthesis of the cytokine IFN- γ in T cells, which has been shown to act in mucosal inflammation in celiac disease [23, 27]. Celiac disease is a chronic inflammatory and autoimmune disease that shares several genetic risk loci with type 1 diabetes [23]. A celiac disease GWAS identified associations with the SLC9A4/IL-1R/IL-18 locus, and a variant in this locus further showed an allelic dosage effect that reduced IL18RAP mRNA expression in celiac patients, indicating a direct link between this locus’ genotype and biologically relevant gene expression in a chronic inflammatory disease [23]. This reduced IL-18 cytokine expression in celiac disease is consistent with the overall reduced expression of SLC9A4 observed among diabetic neutrophils in our study. Another GWAS linked SLC9A4 with regulation of IL-33/ST2 (suppression of tumorigenicity 2, part of the IL-1R family) signaling, which has previously been implicated in both immune and inflammatory diseases, including cardiovascular disease [24]. In addition, a variant in SLC9A4 was associated with eosinophil, though not neutrophil, cell counts in a GWAS of white blood cell subtypes utilizing subjects from the BioBank Japan Project, further indicating a direct role for SLC9A4 in immune cell activity [28].

NECTIN2 is a cholesterol responsive, adherens junction plasma membrane glycoprotein [29, 30]. NECTIN2 is also involved in extravasation and angiogenesis, and may act at vascular inflammation sites to regulate transendothelial leukocyte migration [29, 30]. This study is the first exploration of NECTIN2 expression in neutrophils. We observed downregulation of NECTIN2 among T2D subjects, which is consistent with a continuous cell migration phenotype. Further, we observed higher cholesterol levels among our T2D subjects, consistent with both the role of NECTIN2 and the known propensity for high cholesterol in diabetic patients. Previous evidence has shown that NECTIN2 also regulates immune system responses, including modulating T cell signaling and viral pathogenesis [31]. Further, NECTIN2 has shown potential associations with several diseases, including chronic inflammatory diseases.

PLPP3 is a plasma membrane adherens junction phospholipid phosphatase that converts phosphatidic acid to diacylglycerol and has roles in lipid metabolism and Wnt signaling [21, 32, 33]. The Wnt signaling pathway is involved in a menagerie of cellular functions and development, and its dysfunction has been linked to diseases including cancer, metabolic syndrome, diabetes, and diabetic neuropathy [34]. It has also previously been demonstrated that silencing PLPP3 in an endothelial cell line enhanced inflammatory cytokine secretion, leukocyte adhesion, cell survival, and migration, while its overexpression reversed those effects and instead induced apoptosis [35]. This is consistent with our observation of reduced PLPP3 expression among T2D neutrophils, as PLPP3 normally acts in a protective role as a negative regulator of inflammatory cytokines and leukocyte adhesion. In addition, GWAS have previously associated variants in PLPP3 with several chronic inflammatory diseases, including coronary artery disease [36-39] and eosinophilic esophagitis [40]. Thus, the lipid regulatory and anti-inflammatory roles of PLPP3 in chronic inflammatory diseases lend further evidence for a metabolic-immune axis in diabetes [8].

It is worth note that our study of neutrophil transcriptomics found that most significant gene expression was downregulated in T2D neutrophils compared to healthy neutrophils. Out of several significant genes, the only upregulated FDR-corrected gene in T2D neutrophils was GTSCR1, a non-coding RNA with unknown biology. While non-coding RNA may act to alter gene expression, it is impossible to speculate on a role for GTSCR1 in T2D given the overall paucity of knowledge about its function. However, initial investigations into GTSCR1 showed it may have importance in cardiac inflammation [18]. In Th17-promoting conditions, this was an important gene in regulating T-cell regulation and plasticity [18]. Interestingly, in our study, promoter gene RN7SL2 demonstrated a dose-response to RvE1, showing that neutrophils were viable in cell culture and not terminally differentiated [41].

It has traditionally been thought that gene expression in diabetic individuals would be increased, particularly among inflammatory and immune genes, and indeed gene expression studies of pancreatic islet cells showed increased gene expression in diabetes [42]. However, as neutrophils experience severe dysregulation in individuals with chronic inflammatory disease [4, 5, 9], we do not believe it is surprising to see a downregulation of many key immune and inflammatory genes in neutrophils. Diabetic neutrophils are known to exhibit impaired immune responses and cell migration [4, 5, 9], and we observed a marked reduction in neutrophil chemotaxis in T2D compared to healthy neutrophils (data not shown), consistent with the downregulated gene expression observed in this study. Further, a recent study in adipose tissue showed that apolipoprotein M expression was decreased in T2D and obese individuals compared to lean individuals [43], further indicating a role for downregulation of genes in diabetes pathogenesis. A recent study investigating neutrophil gene expression in a small number of newly diagnosed diabetics compared to controls also observed neutrophil dysregulation, with > 15% of neutrophil genes differentially expressed in recently diagnosed T2D, ~ 60% of which were upregulated [44]. These vast differences in gene expression may be explained by neutrophil heterogeneity in a small number of samples, failure to fully exclude monocytes, or differences in neutrophil gene expression at onset of T2D compared to established disease.

When we further explored neutrophil gene expression by biological pathways, we observed strong links between lipid- and inflammatory-associated genes and pathways, further emphasizing the role of the metabolic-immune axis in diabetes [8]. Plasma membrane cellular pathways were also over-represented, consistent with the roles of our top differentially expressed genes. Hierarchical clustering of over-represented and biologically relevant pathways by individual expression of significant differentially expressed genes tended to sort individuals by disease status (Fig. 2). This is in agreement with our PCA clustering for disease status (Supplementary Fig. S1), where type 2 diabetic and healthy individuals tended to cluster in different regions, consistent with hyperglycemia among diabetics. However, disease heterogeneity could be observed in the clustering, with some T2D individuals in a healthy cluster and vice versa, indicating the potential to identify subsets of disease severity through neutrophil transcriptomics.

To investigate the functional activity of T2D and healthy neutrophils and whether diabetic neutrophil dysregulation could be restored to a healthy, rather than diseased, phenotype, we treated cell cultured neutrophils from T2D and healthy individuals with the small lipid ligand resolution mediator RvE1. We previously investigated the impact of RvE1 perturbation in a murine model of T2D, where RvE1 treatment was able to rescue impaired neutrophil phagocytosis [45]. In this study, we observed that RvE1 dose-dependently modified T2D neutrophil gene expression in cell culture, suggesting a tendency for RvE1 treatment to shift neutrophils toward a healthy phenotype, whereas healthy neutrophil gene expression required a minimum amount of the RvE1 ligand for perturbation, consistent with their lack of chronic inflammation for RvE1 to resolve and consistent with our previous receptor-ligand findings [12]. Diabetic genes were mostly affected at 10–100 nM of RvE1, demonstrating that diabetic neutrophils are not permanently impacted by the chronic disease, and inflammation could be modulated in diabetes. While we showed that multiple genes were impacted by RvE1 treatment, after FDR correction only two genes were significantly differentially expressed between T2D and healthy neutrophils: LILRB5 and AKR1C1.

LILRB5 is an immune system gene that can bind major histocompatibility complex (MHC) class I molecules on antigen-presenting cells, including neutrophils, inhibiting stimulation of an immune response [46, 47]. The increased LILRB5 expression in T2D neutrophils may represent a functional neutrophil change, possibly reflecting the overall neutrophil dysregulation observed in T2D [47]. AKR1C1 catalyzes the reaction of progesterone to its inactive form but may also regulate inflammatory cytokine signaling pathways [48], consistent with some altered cytokine concentrations observed in this study. As LILRB5 and AKR1C1 remained significantly differentially expressed regardless of RvE1 treatment dose, they may reflect molecular changes in neutrophil inflammatory profiles.

In general, the genes differentially expressed between untreated T2D and healthy neutrophils differed between the human serum and cell culture models. This may be due to in vivo vs. ex vivo neutrophil differences, inter-individual variation in neutrophil responses, as not all individuals whose neutrophils were cultured were also included in the baseline serum neutrophil analysis, or they may represent different genes targeting similar pathways. Indeed, when we investigated the top KEGG and GO pathways associated with differential gene expression between T2D and healthy cell cultured neutrophils with a low (1 nM) dose of RvE1, the top pathways were primarily immune related, though there was less of a lipid pathway role (data not shown).

Interestingly, three genes were differentially expressed in both neutrophil cell culture without RvE1 treatment and the baseline serum analysis, including NECTIN2, HTRA3, and ABCG1. Importantly, treatment of the T2D neutrophils in cell culture with a clinically relevant dose, 100 nM, of exogenous RvE1 was able to reduce both the overall number of genes differentially expressed relative to healthy neutrophils and the level of differential expression of NECTIN2 and ABCG1, suggesting a partial normalization of these neutrophils with RvE1 treatment. The top 50 significant genes from neutrophil gene expression in type 2 diabetic versus healthy subjects are listed in Supplementary Fig. S2.

We then used cytokine profiling to further interrogate the functional effect of RvE1 treatment on T2D and healthy neutrophils. We observed trends in differences in cytokine profiles between T2D and healthy neutrophils both with and without RvE1 treatment. In general, T2D neutrophils tended to have a stronger response to RvE1 treatment, particularly at higher, clinically relevant, RvE1 doses (Fig. 3). Without any RvE1 treatment, we observed higher levels of pro-inflammatory cytokines P-Selectin and TNF-α among T2D compared to healthy neutrophils, consistent with recent evidence for enriched cytokine expression in T2D neutrophils [44], though levels of pro-inflammatory MIP-1β, IL-8, and sICAM-1 cytokines were decreased in T2D neutrophils (Fig. 4). As cytokines have complex and often overlapping roles in pro- and anti-inflammatory pathways, the relative contribution of each cytokine to neutrophil dysfunction in T2D will require further investigation. We also observed differences in cytokine profiles by subject, consistent with known inter-individual variation [49]. Our study aimed to understand unprimed and baseline levels of neutrophils in order to allow investigators to map exactly how neutrophils from diabetic subjects behave at baseline levels. We have previously shown that baseline levels of intracellular neutrophil signaling was detectable as soon as 15 mins, demonstrating the feasibility of our premise [12]. In the current study, we further aimed to investigate early stages of neutrophil cytokine production, where we only detected limited amounts of cytokines; a more detailed subgroup analysis was not feasible due to limitations of sample size. However, even at this sample size, this novel cell type has shown distinct differences in transcriptomic and functional cytokine profiling between T2D and healthy individuals, including following RvE1 treatment, indicating that larger follow-up studies are warranted. While intracellular changes were mostly detected after RvE1 treatment, cytokine secretion was limited in the early stage settings. In the future, it will be possible to clarify the differences in cytokine secretion between T2D and healthy volunteers in the presence of RvE1 by addition of stimulus or incubation for a longer time period. Importantly, our cytokine profiling experiments showed that both T2D and healthy neutrophils were viable and capable of both producing cytokines and responding to RvE1 treatment, demonstrating the importance of concentration-based therapeutics for exogenous lipid ligands. Notably, production of TNF-α was clearly increased in T2D neutrophils and modulated by RvE1 treatment. This cytokine trend was corroborated by gene regulation of CD177, demonstrating that that an intracellular regulator of TNF-α was also differentially expressed in diabetic neutrophils when compared to healthy controls (p < 0.05). Absolute values of cytokine expression are listed in Supplementary Tables S3 and S4.

Key strengths of our study include the novel exploration of a biologically relevant cell type (neutrophils) in chronic inflammatory diseases, investigations of both gene level transcriptomics and pathway analyses, and the further exploration of how a key molecule in inflammation resolution (RvE1) impacts gene expression and functional cytokine levels. However, our study does have some limitations, including limited sample size, which made stratification on demographic and clinical characteristics unfeasible. While there were some differences in average demographics between T2D and healthy individuals in our cohorts that align with known risk factors for T2D, PCA of demographic and clinical characteristics (including age, BMI, and ethnicity) did not show any evidence that these differences explained the neutrophil gene expression results, though some residual confounding remains possible. Our T2D subjects were significantly older than our healthy subjects on average. While older individuals experience immune system changes that may include cell deterioration and DNA damage, the similar live neutrophil and monocyte cell counts for T2D and healthy individuals indicate we were able to capture healthy neutrophils from both groups despite any age-related effects on the immune system. A higher proportion of our diabetic subjects were African-American compared to our healthy subjects, who were primarily Caucasian, reflecting known differences in T2D risk across ethnicity. While racial and ethnic population stratification can bias DNA sequencing study results, this was an exploratory study of neutrophil gene expression and PCA did not show a significant impact of ethnicity, though some residual confounding remains possible.

Overall, we showed that neutrophils may act differently in individuals with chronic inflammatory diseases, specifically diabetes, compared to healthy individuals. Further validation of proteomic differences between diabetic and healthy individuals could elucidate important inflammation mechanisms and potential pathways related to resolvin treatment effects, including in subgroups of disease with clinically relevant diabetes comorbidities, such as those with cardiovascular disease or periodontitis. Unraveling these mechanisms of neutrophil dysregulation in chronic inflammatory diseases could ultimately elucidate inflammation and resolution targets for better diagnostic and therapeutic treatment options.

Materials and methods

Study subjects

Subject recruitment has been described previously [12]. Briefly, T2D and healthy subjects were recruited from the patient cohorts of the Center for Clinical and Translational Research at the Forsyth Institute under Forsyth Institute Institutional Review Board-approved protocols (Protocol #11–03 and #13–07). All experiments were performed in accordance with the relevant human subjects ethical guidelines and regulations. All subjects gave signed informed consent prior to study evaluations. Clinical periodontal data and peripheral venous blood were collected. Additional information was collected on subject demographics (age, sex, self-reported ethnicity, and self-reported smoking status), BMI (kg/m2), blood total cholesterol, point-of-care blood glucose, and percent HbA1c [50]. T2D was diagnosed by the subject’s primary care physician according to American Association of Diabetes guidelines [17]. For diabetic subjects, HbA1c was used to determine the level of glycemic control. Neutrophil and monocyte cell counts were determined by lab assay (described below) [12]. Individuals were excluded if they had taken insulin sensitizers, nonsteroidal anti-inflammatory drugs, or antimicrobials within 3 months of study initiation. For this transcriptomics study, a total of 13 T2D and 8 healthy subjects were included for analysis, all of whom were unrelated and over 18 years of age (range: 29–70 years of age). This included 11 T2D and 7 healthy subjects in the main transcriptomics analysis of serum neutrophils, as well as 3 T2D and 3 healthy subjects whose neutrophils were cultured for perturbation experiments (1 T2D and 2 healthy subjects were included in both the main and perturbation analyses). All cell culture and transcriptomics analysis were performed at the Forsyth Institute. Data analysis and cytokine experiments were completed at the J. Craig Venter Institute.

Human neutrophil isolation and cell culture

Isolation of human neutrophils from whole blood was conducted by Ficoll-Histopaque density-gradient centrifugation (Histopaque-1077 and Histopaque-1119; Sigma-Aldrich), as has been described previously [12]. Briefly, neutrophils were isolated after isotonic lysis of red blood cells and counted with a hemocytometer; neutrophil viability (> 95%) was assessed using trypan blue live/dead staining [51]. Monocyte exclusion was completed by plating. A separate monocyte fraction was isolated for monocyte counting. For culture experiments, isolated neutrophils were incubated with RPMI 1640 medium (Sigma-Aldrich) supplemented with 10% FBS (v/v) (Life Technologies) at 37 °C. Wright-Giemsa staining was used to identify individual cell types and confirm neutrophil isolation and purity (routinely > 98%) [12].

Resolvin perturbation experiments

After cell isolation, cultured neutrophils from T2D subjects (N = 3) and healthy controls (N = 3) were treated with one of four doses of exogenous RvE1 (0 nM, 1 nM, 10 nM, or 100 nM) in RPMI and cultured for 1 h at 37 °C. Following incubation, cells were spun down and divided in half for protein and RNA-seq analyses, respectively. Cell culture supernatant was frozen at − 80 °C and extracted for subsequent cytokine analyses. For the perturbation experiments, gene expression was compared between exogenous RvE1 dose-response treatments in 1) T2D neutrophils only, 2) healthy neutrophils only, and 3) T2D versus healthy neutrophils.

Cytokine analyses

Frozen (− 80 °C) cell culture supernatant from the perturbation experiments (N = 3 T2D and N = 3 healthy subject cultured neutrophils perturbed at 0 nM, 1 nM, 10 nM, or 100 nM RvE1) was brought up to room temperature in RPMI and assayed by the Invitrogen human inflammation 20-plex ProcartaPlex cytokine panel (Thermo Fisher Scientific) using a Luminex 200 instrument (Luminex). The 20 assayed cytokines were: MIP-1α, IL-1β, IL-4, IP-10, IL-6, IL-8, IL-10, IL-12p70, IL-13, IL-17A, IFN-γ, GM-CSF, TNF-α, MIP-1β, IFN-α, MCP-1, P-Selectin, IL-1α, sICAM-1, and E-Selectin. Following manufacturer protocols, all samples were run on a plate with 7 standards (diluted 1:4) and a control (RPMI only) [52]. All samples, standards, and controls were run in duplicate.

Quality control (QC) steps were conducted using xPONENT 4.2 software (Affymetrix eBioscience, San Diego, USA) according to manufacturer recommendations. Briefly, any standards with < 70 or > 130% recovery of beads were invalidated. Samples were also checked to ensure they had a bead count of > 30 beads recovered (all samples had > 100 beads recovered; none were excluded at this step). Results were reported as average pg/mL for all measured cytokines following QC. Values at or below the lower limit of quantification (LLOQ) for each cytokine (based on the standard curve after QC) were reported at the LLOQ (the average value of the lowest validated standards). LLOQ for the cytokines (in pg/mL) were as follows: MIP-1α = 7.18, IL-1β = 8.84, IL-4 = 155.02, IP-10 = 4.67, IL-6 = 29.64, IL-8 = 2.47, IL-10 = 8.37, IL-12p70 = 12.23, IL-13 = 5.61, IL-17A = 30.68, IFN-γ = 12.28, GM-CSF = 231.09, TNF-α = 9.74, MIP-1β = 30.76, IFN-α = 2.10, MCP-1 = 15.01, P-Selectin = 1089.78, IL-1α = 2.92, sICAM-1 = 445.56, and E-Selectin = 2894.92.

To compare cytokine levels between groups, we conducted unpaired t-tests with a significance threshold of p < 0.05 and no assumption of consistent standard deviation.

RNA sequencing analysis

RNA-seq was conducted identically for both the serum and cell culture neutrophil experiments. Briefly, RNA was trizol extracted and quantified by Bioanalyzer (Agilent) to confirm RNA-quality (RIN score cutoff of 7.5). Library preparation was conducted using a TruSeq stranded mRNA kit (Illumina), with quality confirmed by Bioanlyzer and sequencing conducted with a Nextseq 500 high output v2 kit on a MiSeq instrument (Illumina). Resulting RNA-seq data was evaluated through a standard bioinformatics pipeline. Briefly, paired-end Illumina sequencing reads were mapped to GRCh38, after quality filtering and trimming, and quantified using STAR v.2.2.0.1 [53] with all default parameters (except “outFilterMismatchNoverLmax” set at 0.05). Reads mapped to each gene in the genome were quantified using htseq-count v.0.9.1 [54] with the following parameters: “-f bam -r pos -s yes”.

Differential gene expression analysis was conducted using edgeR v.3.18.1 [55]. Genes were excluded from subsequent analysis if they did not have a count per million (CPM) > 4 in at least 4 samples. The weighted trimmed mean of M-values (TMM) method [55] was used for normalization. Common dispersion was estimated with parameters method = “deviance”, robust = TRUE. Trended dispersion and tag-wise dispersion was estimated with default parameters. Gene counts were fitted with a negative binomial generalized log-linear model. Likelihood ratio tests were conducted for the specified contrasts of the coefficients, and multiple testing was subjected to a false discovery rate (FDR) correction.

Genes of interest were highlighted for visualization using R library “beanplot” [56]. Heatmaps highlighting significant genes used relative abundances built using R library “vegan” [57] and rendered using R library “heatmap.2”.

Pathway analyses

We also conducted pathway analyses using iPathwayGuide (Advaita Bio). For the main analysis, we subjected the list of all FDR-corrected p-values to association with pathway networks using a threshold of 0.05 for statistical significance (p-value) and a log fold-change (FC) of expression with an absolute value of at least 0.6. Pathways were analyzed from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Release 84.0+/10–26, Oct 17) [58, 59] and Gene Ontology (GO) Consortium database (2017-Nov6) [19, 20], as well as investigating the network of regulatory relations from BioGRID: Biological General Repository for Interaction Datasets v3.4.154. October 25th, 2017 [60], and diseases from the KEGG database [58, 59].

Over-represented (uncorrected p < 0.05) and statistically significant associations (Bonferroni or FDR-adjusted p < 0.05) with pathway classifications for “disease pathways”, “biological processes”, “molecular functions”, and “cellular components” were captured. The underlying differentially expressed genes from the strongest pathway classifications were rendered and expression levels presented in heatmaps using the R library “heatmap.plus”.

For the heatmaps, additional KEGG terms of interest were selected based on previous literature regarding the development of T2D and dental outcomes [3, 8, 13]. These additional terms were highlighted as part of significant networks from the iPathwayGuide analysis, such as the JAK-Stat signaling pathway. Relative abundances for the subset of genes assigned to these terms in the KEGG database [58, 59] were used for rendering to labelled heatmaps.

Additional file 1: Supplementary Table S1. Over-represented KEGG pathways. Supplementary Table S2. Top 5 statistically significant over-represented diseases. Supplementary Table S3. Average cytokine levels (pg/mL) by matched RvE1 treatment of T2D (N = 3) and healthy (N = 3) neutrophils in cell culture. Supplementary Table S4. Average cytokine levels (pg/mL) by differing RvE1 treatment of T2D (N = 3) and healthy (N = 3) neutrophils in cell culture. Supplementary Fig. S1. Principal component analysis (PCA) plot of type 2 diabetic and healthy subjects. Supplementary Fig. S2. Neutrophil gene expression in type 2 diabetic versus healthy subjects for all 50 significant genes. Supplementary Fig. S3. Top 50 significant neutrophil genes differentially expressed in type 2 diabetic and healthy subjects by top biological pathways. Supplementary Fig. S4. Neutrophil gene expression of type 2 diabetic and healthy subjects by biologically relevant KEGG pathways.