Loss of Atg2b and Gskip Impairs the Maintenance of the Hematopoietic Stem Cell Pool Size

ABSTRACT A germ line copy number duplication of chromosome 14q32, which contains ATG2B and GSKIP, was identified in families with myeloproliferative neoplasm (MPN). Here, we show that mice lacking both Atg2b and Gskip, but not either alone, exhibited decreased hematopoiesis, resulting in death in utero accompanied by anemia. In marked contrast to MPN patients with duplication of ATG2B and GSKIP, the number of hematopoietic stem cells (HSCs), in particular long-term HSCs, in double-knockout fetal livers was significantly decreased due to increased cell death. Although the remaining HSCs still had the ability to differentiate into hematopoietic progenitor cells, the differentiation efficiency was quite low. Remarkably, mice with knockout of Atg2b or Gskip alone did not show any hematopoietic abnormality. Mechanistically, while loss of both genes had no effect on autophagy, it increased the expression of genes encoding enzymes involved in oxidative phosphorylation. Taken together, our results indicate that Atg2b and Gskip play a synergistic effect in maintaining the pool size of HSCs.

MPNs (11). Genetic analyses of the families revealed that germ line duplication of ATG2B and GSKIP conferred a risk of familial myeloid malignancies, and overexpression of ATG2B and GSKIP in HSCs enhanced hematopoietic progenitor differentiation (11). Both genes were also highly expressed in de novo acute myeloid leukemias (AMLs) (12). Meanwhile, a genetic analysis of the pedigree of a North American family with MPNs recently showed that germ line duplication of ATG2B and GSKIP genes was not required for the development of familiar myeloid malignancy syndromes associated with duplication of chromosome 14q32 (13). The involvement of ATG2B and GSKIP in hematopoiesis and the development of MPNs and/or AMLs is thus controversial.
The protein encoding ATG2B is a core ATG protein that participates in autophagy (14). Autophagy is a system in which an isolation membrane/phagophore that forms in the vicinity of the endoplasmic reticulum (ER) is elongated and a portion of its cytoplasm is sequestered into autophagosomes, which are then transported to lysosomes for degradation (15). The proteins involved in the formation of autophagosomes are called core ATG proteins, each of which consists of six functional units. When the ULK1 protein kinase complex1 is translocated to the ER subdomain, PI3K complex I is recruited and PI(3)P production increases. ATG18 homologues WIPI1 to WIPI4 bind to PI(3)P and accumulate with their binding partners, ATG2A and ATG2B. ATG2A and ATG2B may have redundant autophagy-related functions in mammals (16), and both connect the ER to the isolation membrane/phagophore and transport lipids (16,17). ATG9A, a membrane protein, transiently accumulates in the isolation membrane/phagophore and scrambles phospholipids transported by ATG2 from the ER to the cytoplasmic layer of the isolation membrane/phagophore (18,19). ATG12 and ATG5 are covalently bound to each other via a ubiquitin-like conjugation reaction (20,21). The ATG12-ATG5 conjugate forms a complex with ATG16L and localizes to the isolation membrane/phagophore, which determines the location of amide bond formation between LC3 family proteins and phosphatidylethanolamine (PE). An LC3-PE protein called LC3-II localizes on the isolation membrane/phagophore and promotes autophagosome formation through membrane perturbation (22).
Gskip encodes glycogen synthase kinase 3b (GSK3b) interaction protein (GSKIP). This protein is an A-kinase anchoring protein for GSK3b and protein kinase A (PKA) that regulates or facilitates their kinase activity toward their targets (23)(24)(25). GSKIP directly interacts with GSK3b, a component of the canonical Wnt signaling pathway, and plays a critical role in embryonic development and functions as a negative regulator of GSK3b (26). Gskip-deficient mouse embryos showed incomplete closure of the palatal shelves accompanied by delayed ossification along the fusion area of the secondary palatal bones, probably due to modulation of GSK3b, causing lethality at birth (27).
Copy number variation of chromosome 14q32 upregulates the gene expression of WIPI1, an ATG18 homologue, and is associated with increased conversion of LC3-I to LC3-II, suggesting that autophagic activity is increased in patients with MPNs (11). ATG2B is directly involved in autophagy, and GSK3b also regulates autophagy through AMPK activation and/ or mTORC1 inactivation (28). However, no studies have demonstrated a biological association between ATG2B and GSKIP. Here, we investigated how simultaneous loss of the Atg2b and Gskip genes affected maintenance of the pool size of HSCs in mice, as well as autophagy and gene expression in Atg2B Gskip-deficient human leukemia cell lines.

RESULTS
Generation of Atg2b Gskip double-knockout mice. As in the human genome, mouse Atg2b is located next to Gskip on chromosome 12, and both exon1 regions are adjacent (Fig. 1A). To investigate the function of Atg2b and Gskip in mice, a targeting vector was constructed by insertion of loxP sequences prior to and after exon1 of the Atg2b gene (prior to and inside exon1 of the Gskip gene) (Fig. 1B). Germline transmission was confirmed by Southern blotting (Fig. 1C). To delete the neomycin resistance gene, mice with mutant alleles were bred with PGK-FLPo transgenic mice (MGI identifier [ID] 4415609). The transgenic mice expressed the mouse codon-optimized FLP recombinase under the direction of the mouse phosphoglycerate kinase 1 promoter. When crossed with a strain containing an FLP recombination target (FRT) site-flanked sequence, FLPo recombinase activity was detected in all cells, with complete recombinase-mediated excision of the target. The resulting progeny containing the Atg2b Gskip flox allele were bred with EIIa-Cre transgenic mice (MGI ID 2137691) to obtain Atg2b Gskip 1/2 heterozygous mice. EIIa-Cre mice carry a cre transgene under the control of the adenovirus EIIa promoter, which targets expression of Cre recombinase to the early mouse embryo. These mice are useful for germ line deletion of loxP-flanked genes. The heterozygotes were crossbred with each other to generate Atg2b Gskip double-knockout mice. We confirmed loss of Atg2b and Gskip transcripts, as well as that of both proteins, in Atg2b Gskip double-deficient mouse embryonic fibroblasts (MEFs) (Fig. 1D and E). As controls for Atg2b Gskip 2/2 mice, we also developed Atg2b and Gskip single-knockout mice by Cas9-CRISPR technology. We targeted exon40 of Atg2b and exon2 of Gskip so as not to affect the expression of either gene ( Fig. 2A). The generated Atg2b-deficient mice had 16-bp (line 1) and 4-bp (line 2) deletions of exon40, and the Gskip-deficient mice had a 1-bp deletion of exon2 (line1, c.33delT; line 2, c.36delC) ( Fig. 2A). We used line 2 as Atg2b knockout mice and line 1 as Gskip knockout mice, and we verified normal levels of the Atg2b and Gskip proteins in the Gskip and Atg2b knockout mice, respectively (Fig. 2B).
Morphological analysis of Atg2b Gskip double-knockout mice. All Atg2b Gskip 2/2 mice from Atg2b Gskip 1/2 intercrosses died in utero, while Atg2b Gskip heterozygous (Atg2b Gskip 1/2 ) mice were born healthy and fertile, with no noticeable pathological phenotype for at least 2 years. Analysis of Atg2b Gskip 2/2 embryos at different developmental stages showed that the majority of embryos died before embryonic day (E) 16.5 ( Table 1). The Atg2b Gskip 2/2 embryos were not significantly smaller than their heterozygous littermates, but exhibited a severe anemic phenotype, with pale skin and subcutaneous edema along the back. Also, the mutant embryos showed exencephaly (  and B). Meanwhile, mice with single knockout of Atg2b were born at the expectedendelian frequency (Table 1) and were viable and fertile. Like other knockout mice for Atg genes that have homologues (29), these knockout mice showed no noticeable pathological phenotypes for 1 year. Gskip knockout mice were also born at the predicted Mendelian frequency (Table 1) and were viable and fertile, which is inconsistent with the study by Deak et al. showing neonatal death (27). This discrepancy may be due to a difference in the targeting region used. Deak's group designed the targeting vector to delete approximate 1,500 bp, including exon2 and introns 1 and 2 (27). The targeting was accompanied by a 1,500-bp deletion about 15 kbp upstream of exon1 of Atg2b ( Fig.  1), and it may have decreased and/or abolished Atg2b gene expression, as with targeting exon1 of Gskip. These results indicated that knocking out both Atg2b and Gskip, rather than just one of them, causes serious developmental defects in mice.
Phenotype of Atg2b Gskip double-knockout embryos. We focused on hematological studies in the mutant mice since hematopoietic progenitor differentiation is promoted in patients with germ line duplication of ATG2B and GSKIP (11). Liver size in Atg2b Gskip 2/2 mice was smaller than that in control mice (Fig. 3B). Histological analysis using Meyer's hematoxylin and eosin (H&E) staining showed nuclear fragmentation of cells in the livers of double-knockout embryos, but not in the livers of embryos with single knockout of either Atg2b or Gskip (Fig. 3C). To examine whether the loss of both Atg2b and Gskip caused cell death, we carried out immunohistochemical analysis using an antibody against cleaved caspase-3, a hallmark of apoptosis. A marked increase in the number of cleaved caspase-3-positive cells was noted in the livers of Atg2b Gskip 2/2 embryos compared to those of control embryos (Fig. 3D). Furthermore, the viable cells in E13.5 fetal livers were significantly decreased in Atg2b Gskip 2/2 embryos compared to those in Atg2b Gskip 1/1 and Atg2b Gskip 1/2 embryos, and consequently the absolute number of living cells in Atg2b Gskip 2/2 embryos was decreased (Fig. 3E). We also found that the frequency of Ter119-positive erythroid cells in the living cells was significantly reduced in Atg2b Gskip 2/2 embryos compared to those in Atg2b Gskip 1/1 and Atg2b Gskip 1/2 embryos (Fig. 3F). Since erythropoiesis is a matter of the highest priority in fetal hematopoiesis, and presumably in fetal development, we consider that the disturbed erythropoiesis might cause fetal liver hypoplasia and life-threatening anemia in Atg2b Gskip 2/2 embryos.
Concomitant reduction of Atg2b and Gskip genes causes a decrease of hematopoietic stem cells. The number of living lineage-negative (Lin 2 ) cells in Atg2b Gskip 2/2 embryonic livers decreased compared to that in wild-type embryos (Fig. 4A). Accordingly, the number of living Lin 2 Sca1 1 c-Kit 1 (LSK) cells was reduced in Atg2b Gskip 2/2 embryos, although no significant differences in the frequency of LSK cells in Lin 2 cells were observed among the three genotype groups (Fig. 4B). On the other hand, Lin 2 CD34 2 Sca1 1 c-Kit 1 (CD34 neg -LSK) cells, which have the potential for long-term lymphohematopoietic reconstitution activity (30), were markedly reduced not only in number but also in frequency in the living Lin 2 cells in Atg2b Gskip 2/2 embryos compared to those in the wild type (Fig. 4C). It is worth noting that the numbers of living Lin 2 cells and living LSK cells differed between Atg2b Gskip 1/2 embryos; some of these embryos had comparable numbers of cells to those of wild-type embryos, whereas the numbers in other embryos were decreased to the level seen in Atg2b Gskip 2/2 embryos ( Fig. 4A and B, right). Furthermore, the absolute number and frequency of CD34 neg -LSK cells were decreased in Atg2b Gskip 1/2 embryos compared to those in wild-type embryos (Fig. 4C).
To further clarify the hematopoietic phenotype caused by concomitant reduction of Atg2b and Gskip genes, we prepared 18 embryos at E13.5 by in vitro fertilization of male and female Atg2b Gskip 1/2 mice and performed flow cytometric analysis. We divided LSK cells into five fractions, specifically long-term HSCs (LT-HSCs), short-term HSCs (ST-HSCs), and three distinct subpopulations of multipotent progenitor cells (MPP2, MPP3, and MPP4), based on the expression profiles of CD150, CD48, and CD135 (31) (Fig. 4D). MPP2, MPP3, and MPP4 cells are defined as MPPs which are developmentally biased toward erythrocyte/megakaryocyte, monocyte/granulocyte, and lymphocyte lineages, respectively (32). As shown in Fig. 4E, the frequency of LT-HSCs within LSK cells significantly decreased in Atg2b Gskip 1/2 embryos compared to that in wildtype embryos and was further reduced in Atg2b Gskip 2/2 mice (Fig. 4E, left), which was in good agreement with the result shown in Fig. 4C. Thus, HSCs in heterozygous mice appeared to exhibit a haploinsufficiency phenotype. We also found that the frequency of MPP2 cells was significantly reduced in Atg2b Gskip 2/2 mice but was maintained in Atg2b Gskip 1/2 embryos. There were minimal differences in the frequencies of ST-HSCs, MPP3, and MPP4 among the three genotypes. We confirmed that there were no substantial differences in the numbers of HSC subpopulations in embryos with single knockout of either Atg2b or Gskip (Fig. 4F).
We next analyzed Lin 2 Sca1 2 c-Kit 1 hematopoietic precursor (LS neg K) cells and found that the relative percentage of these cells among living Lin 2 cells were largely similar among the genotypes (Fig. 5A), suggesting that the remaining HSCs in Atg2b Gskip 2/2 fetal livers have the ability to differentiate into hematopoietic progenitor cells. The frequencies of granulocyte-monocyte progenitor (GMP) cells, megakaryocyte-erythroid progenitor (MEP) cells, and common myeloid progenitor (CMP) cells gated on Lin 2 cells of Atg2b Gskip 2/2 embryos were comparable to those of Atg2b Gskip 1/2 and wild-type embryos, although the absolute numbers of hematopoietic precursor cells were reduced in Atg2b Gskip 2/2 embryos ( Fig. 5B and C). Predictably, embryos of mice deficient in Atg2b or Gskip alone did not show any abnormalities in hematopoiesis (Fig.  5D). In summary, the concomitant reduction of Atg2b and Gskip genes led to a haploinsufficiency phenotype while complete loss of either the Atg2b or Gskip gene had little impact on the hematopoietic system in mice.
Hematopoietic analysis of adult mice heterozygous for the Atg2b Gskip knockout allele. Despite the finding that the absolute numbers of LSK, LS neg K, CMP, GMP, and MEP cells differed between Atg2b Gskip 1/2 heterozygous mice ( Fig. 4B and Fig. 5A and B), these mice were born at the expected Mendelian frequency. Furthermore, the hematopoietic indices of adult Atg2b Gskip 1/2 mice had the same values as those of Atg2b Gskip 1/1 mice (Fig. 6A).
We next performed flow cytometric analyses of the bone marrow (BM) of adult Atg2b Gskip 1/2 mice. In contrast to the findings of embryonic hematopoiesis, the frequencies of HSC/progenitor subpopulations in individual Atg2b Gskip 1/2 mice were virtually equivalent to those in wild-type mice (Fig. 6B to E). While fetal HSCs have the capacity to rapidly self-renew and produce progeny in order to support hematopoiesis in developing embryos, most adult HSCs are in quiescence (33,34). We speculate that the differences between embryonic and adult HSCs are partly due to differences in the nature of these HSCs.
Profiling of ATG2B-and GSKIP-deficient cells. It has been reported that autophagy is critical for proper hematopoiesis (35,36), HSC mobilization (37), and the survival of adult HSCs in the setting of acute metabolic stress (38). Autophagy also preserves the regenerative capacity of HSCs through metabolic suppression (39). Thus, it is plausible that loss of ATG2B and GSKIP is associated with decreased autophagy, resulting in a deficit in the blood development system. We examined whether the loss of ATG2B, GSKIP, or both had an effect on autophagy. To this end, we used CAS9-CRISPR technology to delete one or both genes in the human leukemia cell line K562 (Fig. 7A). The levels of LC3-II and S349-phosphorylated p62, both of which are degraded by autophagy, were comparable between both mutant K562 cell types and wild-type K562 cells (Fig. 7B). An autophagy flux assay with bafilomycin A 1 (Baf A 1 ), which is an inhibitor of lysosomal acidification, revealed that in both mutant K562 cell types, treatment with Baf A 1 caused upregulation of LC3-II and S349-phosphorylated p62 to a similar extent to that in parental K562 cells (Fig. 7B), indicating intact autophagy.
Finally, we conducted transcriptome sequencing (RNA-seq) analysis with wild-type, ATG2B, GSKIP, and ATG2B GSKIP double-knockout K562 cells. We used two pipelines, Kallisto and RNA-Seq by Expectation Maximization (RSEM), to quantify 35,619 transcripts and 26,670 genes, respectively. Multiple comparisons after analysis of variance (ANOVA) revealed that the numbers of genes differentially expressed in the double-knockout K562 cells were 731 with Kallisto and 730 with RSEM. Hierarchical clustering exhibited gene clusters uniquely expressed in each mutant K562 cell type (Fig. 7C). Ingenuity pathway analysis (IPA) with transcripts and genes specified by Kallisto and RSEM identified a significant increase in oxidative phosphorylation activity in ATG2B GSKIP double-knockout K562 cells ( Fig. 7D; see also Table S1 in the supplemental material). Taken together, we concluded that while loss of ATG2B and GSKIP had a minimal effect on autophagy, it was accompanied by increased gene expression of enzymes involved in oxidative phosphorylation.
Excessive oxidative phosphorylation may lead to an increase in reactive oxygen species (ROS) and finally cause cell death by apoptosis (40). We therefore evaluated the levels of apoptosis in hematopoietic stem/progenitor cells. The concurrent loss of Atg2b and Gskip increased the number of annexin V 1 7-aminoactinomycin D (7-AAD) 1 late apoptotic cells in LSK, CD34 neg -LSK, CMP, GMP, and MEP fractions (Fig. 8). In addition, the numbers of annexin V 2 7-AAD 2 living cells in these lineages were markedly decreased in Atg2b Gskip 2/2 mice compared to those in Atg2b Gskip 1/1 mice, except in the GMP fraction (Fig. 8), suggesting that the apoptotic cell death of HSCs and immature progenitors might be at least in part relevant to the decreased number of hematopoietic cells in Atg2b Gskip 2/2 mice. In contrast to these findings, neither an increase in apoptotic cells nor a decrease in live cells was evident in these lineages in the fetal livers of Atg2b Gskip 1/2 mice (Fig. 8), although the number of HSCs was significantly decreased in Atg2b Gskip 1/2 mice (Fig. 4). We consider that concomitant haploinsufficiency of Atg2b and Gskip may contribute to maintaining the pool size of HSCs.

DISCUSSION
It has been suggested that regions with relatively hypoxic conditions harbor the most quiescent HSCs, rather than other hematopoietic progenitors and lineage-committed precursor cells (41)(42)(43). Glycolysis is the most important metabolic pathway in quiescent HSCs, while increases in reactive oxygen species (ROS) levels due to switching of the energy generation mode to mitochondrial respiration induces stem cell differentiation (44)(45)(46). The metabolic status in HSCs is determined by many signaling pathways (44,47). Therefore, the regulation of metabolic conflict between glycolysis and mitochondrial respiration is a key to balancing quiescence and proliferation/differentiation in HSCs. We showed here that the concomitant reduction of Gskip and Atg2b led to reduced pool size of HSCs in mice, while there were no remarkable changes in the hematopoietic system caused by single loss of either Gskip or Atg2b. Simultaneous knockout of Gskip and Atg2b increased the gene expression of enzymes involved in oxidative phosphorylation. We consider that the energy metabolism pathway of HSCs may be shifted from glycolysis to oxidative phosphorylation in mitochondria by the concomitant deficiency of Gskip or Atg2b functions, thus altering the fate of HSCs.
In this regard, GSKIP is an anchoring protein for GSK3b and PKA, which are involved in the regulation of multiple cellular signaling pathways, including the Wnt/b-catenin and PI3K/AKT/mTOR pathways (23)(24)(25)(26). Accumulating evidence suggests that the Wnt/ b-catenin and mTOR pathways play an important role in determining the fate of nor- Gene clustering. RNA-seq analysis was performed using wild-type (n = 3), ATG2B single-knockout (n = 3), GSKIP single-knockout (n = 3), and ATG2B GSKIP double-knockout K562 cells (n = 3). Hierarchical clustering was performed on differentially expressed genes and is presented as a heatmap. Right and left heatmaps were created using Kallisto and RSEM, respectively. (D) Ingenuity pathway analysis of the RNA-seq data was used to predict signaling pathway activity. The pathways for which the P value was #0.01 and the Z-score could be calculated are shown. The upper axis shows 2log 10 (P value), and the lower axis shows the ratio between the number of differentially expressed genes and the number of genes in each pathway. The Z-score, which predicts activation (positive values) or inhibition (negative values) of canonical pathways, is shown by the color scale. An absolute Z-score of $2 is considered significant. mal HSCs, including quiescence, self-renewal, and differentiation statuses (48,49), and also help determine the nature of MPN stem cells (3)(4)(5)(6). We speculate that although mice with loss of Gskip alone showed no particular phenotype in the hematopoietic system, probably due to functional redundancy in the signaling network in mice, deficiency of Gskip with concomitant ablation of Atg2b may be over the redundancy threshold. Since autophagy is required for erythroid differentiation (35,36), as well as for the maintenance and mobilization of HSCs (37,38), we considered that loss of ATG2B suppresses autophagic activity, although not completely, thus exerting a synergistic effect with Gskip ablation on the pool size of HSCs. However, autophagic activity in ATG2B knockout and ATG2B GSKIP knockout K562 cells was comparable to that in parental K562 cells (Fig. 7B), probably due to redundant function of ATG2A. Actually, we verified the expression of ATG2A protein in ATG2B knockout and ATG2B GSKIP knockout K562 cells and also demonstrated similar copy numbers of both Atg2a and Atg2b transcripts in mouse embryonic LSK cells (data not shown). ATG2A is localized to isolation membranes/phagophores but also to lipid droplets, and deficiency of both ATG2A and ATG2B causes abnormal enlargement of these droplets (14,50). Since lipid droplets have the ability to scavenge ROS (51-53), cells lacking ATG2B and GSKIP might become vulnerable to ROS.
We identified the haploinsufficiency phenotype in Atg2b Gskip 1/2 fetal livers, although we could not observe differences in BM hematopoiesis between Atg2b Gskip 1/2 and Atg2b Gskip 1/1 mice. This might be due to the differences between fetal and adult HSCs-selfrenewal and proliferation are accelerated in fetal HSCs to produce abundant hematopoietic progeny in the developing embryos, while adult HSCs are largely quiescent and only occasionally enter the cell cycle to maintain hematopoietic homeostasis (33,34). Previous work demonstrated that oxidative metabolic pathways were activated in fetal HSCs compared with adult HSCs, and consequently, total ROS levels were significantly higher in fetal livers than in BM (54). We speculate that ROS levels in HSCs in Atg2b Gskip 1/2 embryos are prone to reach the threshold that enables HSCs to change their cell fate, leading to the phenotypic heterogeneity observed in individual Atg2b Gskip 1/2 embryos. In contrast, baseline ROS levels in adult HSCs are predicted to be substantially increased in Atg2b Gskip 1/2 HSCs, with ROS having limited influence. Although we could not identify the haploinsufficiency phenotype in Atg2b Gskip 1/2 adult mice at 10 weeks of age, we speculate that increased baseline ROS levels may facilitate stem cell ageing, accompanied by the characteristic features of less frequent quiescence, myeloid bias, and DNA damage accumulation in HSCs (55).
Our experiments indicated a genetic interaction between Atg2b and Gskip and a synergistic effect of ATG2B and GSKIP on the maintenance of both HSCs and hematopoietic progenitor cells by a mechanism that is currently unknown and does not involve autophagy regulation. Our observations may provide insight into the molecular mechanisms of familial MPN and AML.
Generation of genetically modified mice. The Atg2b-Gskip-targeting vector was constructed by insertion of loxP sequences prior to and after exon1 of the Atg2b gene (prior to and inside exon1 of the Gskip gene). The targeting vectors were electroporated into mouse RENKA ES cells, selected with G418 (250 mg/mL; Invitrogen, San Diego, CA), and then screened for homologous recombinants by Southern blotting. Southern blot analysis in Atg2b Gskip flox mice was performed by digesting genomic DNA with SpeI (TaKaRa Bio, Inc.) or EcoRI (TaKaRa Bio, Inc.), followed by hybridization to detect wild-type 16.5-kb and flox 9-kb bands or wild-type 10.5-kb and flox 8.9-kb bands.
significance of the change in deviance between the four groups was tested by the likelihood ratio test, and post hoc tests between the three knockout groups and the wild-type group were conducted by Wald's test. All P values generated by the above-described tests were adjusted by the Hochberg method (61). Differentially expressed genes present only in double-knockout mice were defined as those that were significant in double-knockout versus wild-type mice and nonsignificant in both of the singleknockout mice versus wild-type mice, with a significance level set to 5%. Transcripts per million (TPM) values of differentially expressed genes occurring only in double-knockout mice were used to explore specially expressed gene clusters for each group by hierarchical clustering, with the Euclidean distance method and Ward clustering method. The clustering results are shown by heatmap. All statistical tests and preparation of graphics were performed with R v4.0.0 software (https://www.r-project.org/). Ingenuity Pathway Analysis (IPA, released September 2021; Qiagen Redwood City, CA) was used to identify significant canonical pathways involving the differentially expressed genes enriched using the above two methods.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, XLSX file, 0.1 MB.