Reduced expression of OXPHOS and DNA damage genes is linked to protection from microvascular complications in long-term type 1 diabetes: the PROLONG study

Type 1 diabetes is a chronic autoimmune disease requiring insulin treatment for survival. Prolonged duration of type 1 diabetes is associated with increased risk of microvascular complications. Although chronic hyperglycemia and diabetes duration have been considered as the major risk factors for vascular complications, this is not universally seen among all patients. Persons with long-term type 1 diabetes who have remained largely free from vascular complications constitute an ideal group for investigation of natural defense mechanisms against prolonged exposure of diabetes. Transcriptomic signatures obtained from RNA sequencing of the peripheral blood cells were analyzed in non-progressors with more than 30 years of diabetes duration and compared to the patients who progressed to microvascular complications within a shorter duration of diabetes. Analyses revealed that non-progressors demonstrated a reduction in expression of the oxidative phosphorylation (OXPHOS) genes, which were positively correlated with the expression of DNA repair enzymes, namely genes involved in base excision repair (BER) machinery. Reduced expression of OXPHOS and BER genes was linked to decrease in expression of inflammation-related genes, higher glucose disposal rate and reduced measures of hepatic fatty liver. Results from the present study indicate that at transcriptomic level reduction in OXPHOS, DNA repair and inflammation-related genes is linked to better insulin sensitivity and protection against microvascular complications in persons with long-term type 1 diabetes.


Type 1 diabetes is a chronic autoimmune disease requiring insulin treatment for survival. Prolonged duration of type 1 diabetes is associated with increased risk of microvascular complications. Although chronic hyperglycemia and diabetes duration have been considered as the major risk factors for vascular complications, this is not universally seen among all patients.
Persons with long-term type 1 diabetes who have remained largely free from vascular complications constitute an ideal group for investigation of natural defense mechanisms against prolonged exposure of diabetes. Transcriptomic signatures obtained from RNA sequencing of the peripheral blood cells were analyzed in non-progressors with more than 30 years of diabetes duration and compared to the patients who progressed to microvascular complications within a shorter duration of diabetes. Analyses revealed that non-progressors demonstrated a reduction in expression of the oxidative phosphorylation (OXPHOS) genes, which were positively correlated with the expression of DNA repair enzymes, namely genes involved in base excision repair (BER) machinery. Reduced expression of OXPHOS and BER genes was linked to decrease in expression of inflammation-related genes, higher glucose disposal rate and reduced measures of hepatic fatty liver. Results from the present study indicate that at transcriptomic level reduction in OXPHOS, DNA repair and inflammation-related genes is linked to better insulin sensitivity and protection against microvascular complications in persons with longterm type 1 diabetes.
Type 1 diabetes is a life-long disease characterized by the destruction of pancreatic beta cells and lack of insulin production. Type 1 diabetes affects patients of any age, but it is the most prevalent form of diabetes for people younger than 20 years old comprising over 85% of all diabetes cases in this age category worldwide 1 . With insulin treatment, the life expectancy of the people with type 1 diabetes has become comparably longer in relation to Figure 1. Pathways involved in the progression to microvascular complications (red color) and protective mechanisms (blue color) in patients with prolonged type 1 diabetes. Altered glucose metabolism is associated with a disbalance between glycolysis and oxidative phosphorylation (OXPHOS). A shift towards enhanced OXPHOS in rapid progressors is associated with increased uncoupling between oxidation and phosphorylation, generation of reactive oxygen radicals (ROS) in mitochondria, induction of DNA damage response pathways and pro-inflammatory activation of the immune cells. These mechanisms are involved in progression to microvascular damage and metabolic violations increasing risk of PDR and DKD. On the contrary, nonprogressors are characterized by a moderate reduction in expression of OXPHOS genes that correlated with decreased DNA damage and immune cell activation, and linked to better insulin sensitivity and lower measures of fatty liver. These mechanisms may protect persons with type 1 diabetes from microvascular injury and hyperglycemia-related organ damage. The increased expression of TDG and PARP4 in non-progressors highlights the contributing role of both epigenetic and post-translational modification mechanisms in preventing development of diabetic complications. OXPHOS oxidative phosphorylation, ROS reactive oxygen species, UPC2 uncoupling protein 2, FIS1 mitochondrial fission 1 protein, NF-kB nuclear factor kappa B, TNFR4 tumor necrosis factor-alpha receptor 4, IL-1 interleukin 1, TDG thymine DNA glycosylase, PARP4 poly(ADPribose) polymerase family member 4, PDR proliferative diabetic retinopathy, DKD diabetic kidney disease, RP rapid progressors to PDR and/or DKD, NP non-progressors (created with BioRender.com).
Differential expression and functional enrichment of blood transcriptome. A large number of genes was found to be differentially expressed in blood cells between the two groups. A gene was labeled as differentially expressed between the groups after multiple testing adjustment of p-values using false discovery rate (FDR) < 0.05. Overall, there were approximately 5000 differentially expressed genes (DEGs) detected in age-adjusted DE analysis at FDR ≤ 0.05 (Fig. 2). The number of DEGs substantially reduced to 50 genes after including HbA 1c as a covariate in the model. The DEGs between NPs and RPs (FDR < 0.05) are shown in Tables S1 and S2.
To provide insights into the underlying biological processes that were associated with the DEGs between NPs and RPs, we performed pathway and ontology term enrichment analyses. The results of over representation analyses for gene ontology (GO) terms in three categories as biological process (BP), molecular function (MF) and cellular component (CC) and KEGG, and gene set enrichment analyses for GSEA gene sets (Hallmark, Reactome and Transcription Factor Targets) are presented in Supplementary Tables S3 to S11 (Fig. 3).
Notably, antioxidant response components such as NFE2, GPX1, GPX4 genes and overall Reactome Detoxification of Reactive Oxygen Species Pathway (p = 0.006), which is a pathway of antioxidant response, were downregulated in NPs. The lower OXPHOS expression in NPs was associated with reduced expression of apoptosis regulators (BAD, BBC3, BCL2L12, BAX) and uncoupling protein (UCP2) dissociating oxidation and phosphorylation processes. No pathways or terms for apoptosis reached statistical significance.
Additionally, the top terms from gene set enrichment analyses remained the same after HbA 1c adjustment (Table S6-S11). The associations between individual genes in canonical OXPHOS pathway (KEGG) with age, duration of diabetes, HbA 1c , complications type as retinopathy and nephropathy, as well as and BMI are shown in Fig. S4-S12. The change of expression was independent of BMI and age for most of the genes (Fig. S3).  Table S7) and specifically base-excision repair (BER) (KEGG, p = 4.1 × 10 -5 , Reactome, p = 2.1 × 10 -4 ). Expression of genes involved in DNA repair processes was downregulated in NPs. Out of 32 BER genes present in KEGG BER pathway, 50% (15) were downregulated in NPs and positively correlated with OXPHOS genes (Fig. 4). Only 2 genes, poly [ADP-ribose] polymerase 4 (PARP4) and thymine DNA glycosylase (TDG) were upregulated, and the latter was an exception demonstrating strong negative correlation with OXPHOS and most of the genes encoding mitochondrial enzymes (Fig. 4).

Moderate reduction of OXPHOS and BER gene expression is associated with better insulin sensitivity and overall metabolic balance.
We have previously demonstrated that NPs were characterized by higher insulin sensitivity as compared to RPs 10 . In the present study, the top results from the GSEA transcriptional factor targets enrichment analyses identified HMGA1-an important regulator of insulin action and function 15 being statistically and significantly upregulated in NPs (Table S11). Additionally, differentially   (Table S9). Moreover, expression of inflammatory pathway genes that have a role in COVID-19 infection (KEGG, Coronavirus disease-COVID-19, p = 2.6 × 10 -8 ) was also decreased (Table S5). This pathway became available after extensive research on COVID-19 during last year and includes many components of immune response, which shows downregulation of general inflammation in NPs. Importantly and consistently, inflammation-related pathways remained significantly overrepresented after HbA 1c adjustment also in the GO terms (Table S4).

Discussion
The pathophysiology of microvascular damage in patients with type 1 diabetes is complex and multifactorial. Among numerous pathogenic factors, chronic hyperglycemia, endothelial dysfunction, and low-grade inflammation are considered to play a critical role in microcirculatory disturbances and development of vascular complications 16,17 . In physiological situations, there is an equilibrium between the microcirculation and a rate of metabolism in tissues, while this relationship is impaired in the organs affected by diabetes environment exhibiting reduced blood supply (hypoxia) and as a consequence a high metabolic rate and energy consumption.
The key findings in the present study revealed synergistic reduction in expression of OXPHOS, DNA repair and genes regulating immune response pathways in persons with long-term type 1 diabetes who remain largely free from vascular complications. Clinically, this was linked to higher insulin sensitivity and overall better metabolic balance in persons who remain free from complications. OXPHOS is one of the core systems of the metabolic control of cells, which impacts cell survival, growth, proliferation, functioning and death. OXPHOS generates a substantial amount of superoxide derived from electron leaks while substrates are metabolized. Complexes I and III of mitochondrial respiratory chain generate O 2− that is quickly transformed by superoxide dismutase (SOD) into H 2 O 2 with its further neutralized by catalase and glutathione peroxidase 18 . Therefore, mitochondrial energy release comes with the cost of generation of ROS, and consequently oxidative damage. There is a direct link between the intensity of OXPHOS and ROS production 19 and subsequent damaging effects on cells 20 . ROS produced through OXPHOS have also been implicated in neurodegenerative diseases 21 . The retinal pigment epithelium was found to be prone to DNA damage due to oxidative stress 22 . Increase of ROS via mitochondria was found to be an important contributor to the age-related macular degeneration 23 and kidney damage and inflammation in persons with diabetes 24 . However, the evidence of antioxidants' benefits on DKD progression to end-stage renal failure was inconclusive 25 . This could partially be attributed to the protective properties of ROS as signalling molecules under normal physiological conditions 26 . Our present findings are in line with these notions indicating that modest downregulation of OXPHOS may be beneficial to preserve physiological functions of ROS in non-progressors and thereby contributing to organ protection from oxidative stress.
Another role of ROS beyond oxidative cell damage is its genotoxic effect causing DNA damage and altering gene expression 12 . Repair of damaged DNA involves various mechanisms, including base-excision repair (BER), homologous recombination (HR), nucleotide excision repair (NER), mismatch repair (MMR) and non-homologous end joining (NHEJ) 27 . BER is considered to be a key pathway for recognition and repair of diabetes-induced persistent ROS-induced DNA damage 28 . Although the role of OXPHOS-born ROS and DNA damage in the risk of development of diabetic vascular complications is widely recognized, results from the present study illuminate the important crosstalk between mitochondrial transcriptome and the DNA repair machinery in persons with type 1 diabetes resistant to vascular complications. About 40% of the genes encoding BER pathway enzymes were found to be upregulated in RPs and highly correlated with OXPHOS genes. The simultaneous increase in www.nature.com/scientificreports/ both may cause more DNA damage than dysfunctional OXPHOS alone since it was reported that single-strand breaks can be generated by ineffective BER or oxidative stress 29 . However, we cannot rule out that increased expression of OXPHOS genes could be a result of DNA damage as a consequence of defective functioning of DNA N-glycosylase proteins 30 . Besides, a BER-associated DNA-glycosylase, TDG, which plays a key role in active DNA demethylation by recognizing and excising 5-formylcytosine (5fC) and 5-carboxylcytosine (5caC) in the CpG sites, and installing an unmethylated cytosine 31 , was upregulated in NPs and negatively correlated with most of the mitochondrial enzymes. On the contrary, expression of another DNA glycosylase gene, NTHL1, which is essential for BER signatures and oxidative damage present in pancreatic cancinomas 32 showed reduced expression in NPs. Hence, the OXPHOS and DNA repair enzymes are highly interconnected. Their lower expression in NPs contributing to cellular defense against genetic alterations could protect cells against violation of various biological processes such as cell proliferation, survival and apoptosis, cell signaling and transcription affecting cellular functions. An additional explanation can reside in the common regulation of both OXPHOS 33 and BER genes 34 by hypoxia, which is commonly present in diabetes. Hypoxia inducible factor (HIF), which is the main cellular adaptor to hypoxia, is supressed in diabetes 35 . An impaired response of HIF to hypoxia has been described in diabetes with direct pathogenic role on different complications of diabetes 35 . In perfect agreement, gain-of-function genetic variants in the HIF gene that are more active in the diabetic environment are protective for nephropathy and retinopathy in patients with type 1 diabetes 36,37 . Notably, HIF is one of the strongest metabolic signals regulating activity of key glycolytic enzyme pyruvate kinase 2 (PKM2). PKM2 is the key rate-limiting enzyme in the last step of glycolysis and thereby determines glycolytic flux. Recent observations from the Joslin Diabetes Center suggested that enhanced metabolism of intracellular glucose via activation of PKM2 can protect from diabetic nephropathy along with preventing or reversing mitochondrial dysfunction dysfunction 38,39 . By interacting with HIF1a-a master regulator of the reaction to hypoxia, PKM2 plays an important role in metabolic shift of glucose metabolism away from oxidative phosphorylation towards a glycolytic program. Importantly, increased glycolysis could be at play serving both bioenergetics and metabolic signaling functions 40 . Thereby, whether NPs have less suppression of HIF activity (functionally through a better glucose control or genetically through a more active HIF 36,37 ) followed by reduction of OXPHOS and BER pathways warrants investigation.
Another important finding was downregulation of KEGG thermogenesis in NPs, which in addition to downregulated OXPHOS, may cause decreased ROS generation due to oxidation and phosphorylation coupling 41 . Uncoupling proteins are responsible for converting mitochondrial membrane potential to heat production 42 . In our study UCP2 expression was significantly lower in NPs comparing with RPs, reflecting the role of uncoupling in the development of diabetic complications. UCP2, an ubiquitous form of uncoupling proteins, plays an important role in diabetes as contributes to the mitigation of ROS, regulate glucose sensing and insulin release 43 . Due to uncoupling, upregulated expression of OXPHOS genes in RPs can be associated with a decrease of mitochondrial ATP generation stimulating alternative sources of energy production. It is well known that diabetes is associated with changes in levels of lactate and pyruvate, and their ratio in the peripheral blood reflects the balance between mitochondrial oxidation and glycolysis 44 . As previously reported, we observed lower levels of pyruvate in NPs as compared to RPs in this cohort 10 . In the present study, the correlation between differentially expressed OXPHOS genes and pyruvate levels further provide a supportive evidence for an important role of cross-talk between glycolysis and mitochondrial oxidative phosphorylation in the development of vascular complications in persons with type 1 diabetes. Although it is possible that environmental factors might modify this relationship in the present cross-sectional study, another plausible genetic determinant might be genetic intrauterine programming of metabolic responses to hypoxia. Thus, a situation when such a shift of glycolysis and mitochondrial biogenesis has been recognized is during embryonic development whilst trophoblast excluded from the oxygenated maternal blood flow 45 . At this stage, exposure to hypoxia induces stabilization of HIF to support metabolic remodeling towards increased glycolysis acquired during transition of embryonic stem cells 46 . This leads to reduced mitochondrial biogenesis and activation of alternatives to glycolysis, the pentosephosphate-and the sorbitol/polyol auxiliary glucose metabolism pathways-supporting our previous findings in PROLONG 10 and similar to what has been observed in the T1D Medalists' kidneys protected from nephropathy 38,39 .
The third affected process seen in the present study was the downregulation of genes in the NAFLD pathway in NPs. These transcriptomics changes are in support of our previous observations of lower measures of fatty liver in NPs compared to RPs 10 . The key mechanisms of NAFLD are related to insulin resistance and the metabolic syndrome. Our earlier observations of higher insulin sensitivity in NPs 10 was further supported by the enrichment of the upregulated HMGA1-target gene pathway, which is an important transcription factor in regulation of the insulin receptor. Additionally, as discussed above pyruvate is positively and lactate/pyruvate ratio negatively correlated with OXPHOS and BER genes, further emphasizing the interdependency between metabolism (glycolysis), insulin sensitivity and mitochondrial functioning for protection of DNA damage. It is possible that metabolic shift towards glycolytic program and previously proposed increased nucleotide biosynthesis via pentose-phosphate pathway of glycolysis act as denominator(s) of enhanced glucose metabolism and thereby improved insulin sensitivity and lower accumulation of liver fat in NPs. Furthermore, such DNA damage/repair mechanisms are highly correlated to immune response and inflammation, however, it is important to note that the relationship between DNA repair machinery and immune response is bi-directional 47 .
Finally, it is important to emphasize a possible contribution of the present findings to the severity and risk of vascular complications in other types of diabetes. Although pathophysiological factors leading to chronic hyperglycemia in different diabetes subgroups (type 1 and type 2 diabetes, MODY) are distinct, metabolic features such as degree of insulin sensitivity/resistance and fatty liver parameters might underpin common mechanisms related to progression towards or free from vascular complications. In support of this, recent observation of novel class of glucose-lowering agents-inhibitors of SGLT2-known to improve insulin sensitivity by inhibiting Scientific Reports | (2021) 11:20735 | https://doi.org/10.1038/s41598-021-00183-z www.nature.com/scientificreports/ glucose re-absorption in the kidney, shown to be effective not only in patients with type 2 diabetes but also in type 1 diabetes 48 and associated with reduced CKD in non-diabetic population 49 .
Limitations. The analysis was done for cross-sectional transcriptomic measurements. Hence, we could observe association between the forementioned ontology terms/pathways with complication status, but we were not able to study temporal effects neither detect causal relationships for the development of vascular complications. Nevertheless, despite differences in age and HbA 1c between non-progressors and rapid progressors, the main associations between OXPHOS, DNA repair and inflammation-related pathways remain consistently and statistically significant after correcting for these factors as covariates in FDR analyses suggesting independent and robust associations. However, cross-sectional study design only depicts one-time measurements of HbA1c, and do not reflect severe impact of chronic hyperglycemia, which is often described as glycemic memory. Legacy of glycemic memory exhibit long-lasting effects and therefore longitudinal analyses of HbA1c effects would be of high value.

Conclusion
Progression to vascular complications in patients with diabetes is a systemic process, which might be monitored by use of transcriptomic alterations in peripheral blood cells. Synchronous elevation of expression of OXPHOS and DNA repair genes may serve as a hallmark of mitochondrial dysfunction associated with increased ROS production, chronic inflammation and elevated lactate-to-pyruvate ratio levels. In contrast, genetically determined moderate reduction in oxidative metabolism due to possible uncoupling towards enhanced glycolytic program in NPs could be linked to lower DNA damage and lower inflammatory activation of blood cells, and overall better insulin sensitivity. Taken together, present findings suggest that reduction in OXPHOS function might be beneficial by protecting against development of vascular complications in patients with long-term type 1 diabetes.

Research design and methods
Study participants. Information on the PROLONG eligibility criteria and complication status was based on initial electronic hospital records 50  Non-progressors (NPs) were defined as those with diabetes duration of more than 30 years and who did not develop any major complications (diabetic nephropathy, proliferative retinopathy, myocardial infarction, stroke or chronic foot ulcer), whereas rapid progressors (RPs) were defined as those, who developed any of the abovementioned complications within 25 years of diabetes duration. NPs who prior to or at the PROLONG study visit were diagnosed with any of the above-mentioned complications were designated "late-progressors" and excluded from the present analyses. PROLONG study was approved by the local ethics committees (PROLONG-Sweden, Regional Ethics Review Board, Department 1, Lund, Sweden, Dnr 777/2009, PROLONG-Denmark, The Capital Region Ethics Committee, Hillerød, Denmark, Dnr H-2-2013-073, for data analysis of the study in UiB, West regional health authority, Norway, REK-2019/1324) and conducted in accordance with local institutional and national regulations.
Definition of diabetic complications. Diabetic nephropathy was defined as the presence of (1) macroalbuminuria ≥ 200 µg/min in a timed overnight urine collection, (2) an albumin/creatinine ratio > 300 mg/g as macroalbuminuria, or based on a documented diagnosis of diabetic kidney disease, which is diagnosed by using eGFR < 60 mL/min. Proliferative diabetic retinopathy was assessed with fundus photography and defined as the presence of proliferative retinopathy in at least one eye and/or laser therapy (panretinal photocoagulation), or non-traumatic blindness.

Procedures and measurements.
In the PROLONG study, participants fasted overnight (> 8 h) and were instructed not to take their medication in the morning of the visit. Trained diabetes research nurses and biomedical analysts performed the physical examination following standard operating procedures. On the day of the examination, a signed informed consent was obtained, and blood and urine samples taken. The participants were also asked to fill out a detailed questionnaire. These self-reported records were then validated by the diabetes nurse at the clinical research visit and this information was stored in a database at the Scania University Hospital, Malmö, Sweden. RNA sequencing. TruSeq stranded Total RNA kit from Illumina was used for sequencing library generation. Sequencing was done on NextSeq 500 sequencing machine to create paired-end files. The sequencing was done at Lund University, Malmö, Sweden.
Cell content comparison. The whole blood content of the subjects was compared to in the cell counts from laboratory to check whether the blood cell content differ between two groups. The available cell counts for leukocytes and lymphocytes for the patients (87% of NPs and all RPs) were used to compute difference of cell fractions in the samples. It was found that both lymphocyte and leukocyte percentages in total of the blood cells did not differ significantly between groups (p > 0.5). We used two local tools to compare the cell contents of two groups based on microarray expression signatures for whole blood to confirm this (Fig. S13-S14 www.nature.com/scientificreports/ NASeq was used with the signature coming from ref. 51 and quanTiseq was used with its internal signature 52,53 . Both signatures consistently showed no significant differences between cell contents across two groups, except natural killer (NK) cells, which seems to be slightly higher in NPs. Both signatures were created from microarray expression data and conversion from microarray ids to gene ids for the comparison of our own RNA sequencing results was not complete due to duplicative or non-mapping between two id types.
Processing and analysis of RNAseq data. Total sample size was 154 and the average library size for the included samples were 29.6 ± 15.2 M. The quality of paired-end RNAseq files were checked with MultiQC v1.0 by using fastq files. 22 samples were failed QC and removed from further analysis. Phred score greater than 30 was achieved for each sequencing position for all other samples. There was no adapter contamination shown in QC (adapter contamination < 0.1%), so trimming was not applied in order not to lose information. The alignment of transcripts was done by using Kallisto v0.43.1 54 with GRCh38.p10 reference assembly with Ensembl Homo sapiens v90 annotation as the reference transcriptome 55 . The downstream analysis from this level was done by using R statistical environment 56 and R Studio (https:// www. rstud io. com). The estimated count values of the transcripts given as Kallisto output were converted to the gene level expression values by using Tximport v1.6.04 57 . The genes that had low counts were excluded from the expression matrix with the inclusion criteria of CPM higher than 0.5 for at least the smallest number of the sample subset (n = 34). The differential expression analysis was done with using edgeR v3.20.73 58 . The raw count values were provided to edgeR since the tool handles with the normalizations for sequencing depth, gene length and RNA composition of the libraries (TMM normalization). edgeR handles unbalanced sample sizes robustly and power of the analysis does not depend on the ratio of sample sizes. The condition is used as explanatory variable in a quasi-likelihood negative binomial generalized log-linear model (glmQLFit function) with using age and/or HbA 1c as covariates. Empirical Bayes quasilikelihood F-tests were used to assess the differential expression (glmQLFTest function). Sex adjustment was not done via adding sex as an explicit covariate since Robust algorithm was used to down-weight sex-linked genes. False discovery rate (FDR) correction was done by using Benjamini-Hochberg method within the package. The gene was considered differentially expressed after adjustment for multiple testing FDR ≤ 0.05. 5153 genes were found to be differentially expressed for age-adjusted differential expression analysis (610 of them with at least 50%-fold-change, Fig. 2, Fig. S1) while the number is reduced to 49 if HbA 1c adjustment is applied alongside age (38 genes with at least 50%-fold-change) with false discovery rate of 0.05. There was no missing data for the variables that are included into this analysis.
Average age of NPs was higher than the RPs by study design. Hence, we checked to see whether there is any difference in the gene expression strongly affected by the age difference between two groups. We used an interaction term between age and complication status (NP and RP) to catch this effect, if there was any, in addition to the other covariates. The result showed a significant change in the expression for only two genes (RPL3P4 and CXCL10) were significant for the interaction between age and complication status, which indicates that the effect of the complication status did not vary as a function of age for the expression of the majority of the genes. For the full design formula of age + condition + age:condition, only RPL3P4, HLA-DRB4, MAYDM and H19 genes became significant for the age term (FDR < 0.05). The expression of the last three genes increased with age while the first one decreased.
Functional enrichment analysis. GO and KEGG pathway enrichments were done by using goana and kegga functions 59 . This was done as an overrepresentation analysis that determined over-represented/enriched biological processes and/or pathways in the list of differentially expressed genes. The genes in the relevant GO terms were extracted from gene ontology (https:// www. geneo ntolo gy. org) and KEGG (genome.jp) databases 60,61 . GSEA enrichment tests were done by using hallmark, all transcriptional factor targets and Reactome genesets and camera test 62 (https:// www. gsea-msigdb. org) and employing camera function of edgeR. Gene set enrichment test is the procedure to find any coordinated changes related to a gene set in the current list of genes without taking significance of differential expression into consider. The relevant statistic, which is fold-change in the comparison of gene expression, is aggregated to find the coordinated behavior related with the gene set. Enrichment tests were performed in Jan 2021. Figures 2, 3, 4 and 5 were generated by using R Studio.