Progression to Islet Destruction in a Cyclophosphamide-Induced Transgenic Model: A Microarray Overview

The BDC2.5/NOD transgenic line has been described previously (3). The mice used in the present experiments were from >25 backcrosses to the NOD/Lt genetic background and were all 6–8 weeks of age and nondiabetic. Both male and female mice were used. Diabetes was induced by an intraperitoneal injection of a single dose of CY (200 mg/kg in PBS).

Islets were isolated by the Islet Preparation Core of the Juvenile Diabetes Research Foundation (JDRF) Center for Islet Transplantation at Harvard using a standard technique (15). Briefly, ice-cold collagenase (Liberase; Roche) was infused into the common bile duct of an anesthetized mouse. The pancreas was then removed and incubated at 37°C for ∼20 min with gentle agitation. Further separation was accomplished with a Histopaque density gradient (Histopaque-1077; Sigma). The islets were handpicked under a stereomicroscope to ensure a pure islet preparation. Total RNA was isolated from these islets using Trizol followed by ethanol precipitation.

Before amplification of islet mRNA or to confirm the DNA chip analysis, transcripts of some genes were quantitated with quantitative real-time PCR. Amplification of specific genes was monitored during the PCR utilizing an internal fluorescent probe on an ABI Prism 7700 Sequence Detection System (Applied Biosystems). For each gene of interest, a TaqMan primer/probe set was designed such that both the amplified fragment and the probe spanned an intron to avoid spurious amplification from genomic DNA. For gene expression analysis, cDNA was made from 100 ng of total pancreatic islet RNA in a 20-μl reverse transcriptase reaction. Six microliters of this reverse transcription product were used for triplicate 20-μl TaqMan reactions. Transcripts of the genes of interest (interleukin [IL]-18, IL-12 p40, interferon [IFN]-γ, and Cxcl9) were normalized relative to transcripts of the ubiquitous housekeeping gene glycerol phosphate dehydrogenase (GPDH). For each experiment, a standard curve was generated for every gene based on δ-Ct (the number of cycles required to generate ∼50% of the maximum amount of amplified product, set during the linear phase of amplification). Normalization of the relative amounts of a gene of interest (X) from RNA from a CY-treated (cy) versus a PBS-treated (ctl) mouse were calculated as: X_cy/X_ctl (normalized) = (X_cy/X_ctl) × (GPDH_ctl/GPDH_cy).

Islet RNA was amplified using a standardized protocol (MessageAmp aRNA kit; Ambion). Briefly, cDNA was synthesized from 1 μg total RNA template by reverse transcription primed by a hybrid oligonucleotide containing oligo-T and T7 RNA polymerase promoter sequences. Double-stranded DNA was then synthesized by incubation with RNase A followed by extension with polymerase I. Multiple copies of antisense RNA were then produced with T7 RNA polymerase (∼100-fold amplification). This RNA was then primed with random hexamers, and the entire amplification process was repeated, producing biotinylated RNA with 10,000-fold amplification of the original mRNA. Biotinylated cRNA was then prepared with biotinylated ribonucleotides in another T7 RNA polymerase reaction (Enzo BioArray HighYield RNA Transcript Labeling Kit; Affymetrix).

Biotinylated RNA was fragmented in a proprietary Affymetrix buffer optimized to break down full-length cRNA by metal-induced hydrolysis (https://www.affymetrix.com/download/manuals/expression_print_manual.pdf). The samples were hybridized to Affymetrix MU74v2A microarray chips in an Affymetrix Fluidics Station 400 hybridizer/analyzer, with streptavidin-phycoerythrin to detect biotinylated probe. Fluorescence on the chip was quantified on an Affymetrix Fluidics Station 400 hybridizer/analyzer and digitized using an Affymetrix GeneChip Scanner with Affymetrix GeneChip operating software. The initial reads were processed through the robust multiarray average (RMA) algorithm (implemented on the Array analyzer; Insightful) for probe-level normalization. These primary values were averaged (with outlier elimination) and an approximate measure of significance calculated (Welch’s approximation t test). For analysis of significance, control datasets approximating the values and variances of the real data were generated by random reshuffling of the data (random draws within rows) or as detailed in results. The raw datasets have been deposited on the GEO databank.

Previous experiments have shown that the nondestructive insulitis of adult BDC2.5/NOD mice is extensive, representing two to three times the β-cell volume (3,11,16–19). This infiltrate is heterogeneous and includes large numbers of B- and T-cells, NK (natural killer) cells, macrophages, and dendritic cells. Thus, we chose to evaluate the global response to CY in the infiltrated islet as a whole, rather than in individual cells. To set the stage for the microarray analysis, we investigated the kinetics of diabetes progression after CY treatment in our current BDC2.5/NOD colony and tested RNA preparation methods over that time period. As expected from previous results (11), diabetes was induced in a highly synchronous fashion after administration of a single intraperitoneal injection of CY (200 mg/kg) (data not shown). Islets were prepared at different times after CY treatment (90–120 islets/mouse), using a conventional islet preparation technique: retrograde intraductal injection of collagenase, density gradient fractionation, and picking of isolated islets by micromanipulation. The yield of islets was similar during the first 3 days after treatment, as has been shown previously (11). This result was consistent with parallel histological analyses that demonstrated maintenance of the islet structure up to day 3, but a breakdown after that (11) (data not shown). Because the islets that remained whole and therefore could be isolated by the collagenase procedure at day 4 were likely to represent “survivors” with a delayed evolution, we limited the gene expression analysis to the first 3 days after treatment. Histological examination showed these islets to be morphologically intact and to contain β-cell islands as well as surrounding infiltrate cells. Several RNA preparation techniques were tested; the best yields were obtained with a Trizol method, and analysis by real-time quantitative PCR showed the RNA preparations to be of good quality.

In order to ensure the validity of the microarray analysis, three or four independent experiments were performed at each time point (Table 1). Each experiment also included an untreated group (hereafter referred to as “day 0”). From previous experiments (11), we knew that IL-12 expression should be induced in samples from mice progressing to diabetes. We used RT-PCR to show that IL-12 p40 transcripts were indeed induced in the samples used for microarray analysis (data not shown). From all samples, RNA was prepared, with an average yield of 80 ng, and labeled probes for microarray analysis were synthesized by using biotinylated ribonucleotides in the final RNA polymerase reaction. These probes were then hybridized to Affymetrix Mu74Av2 oligonucleotide microarrays, which represent 8,063 unique genes, of which 1,935 are expressed sequence tags.

Data were processed with a custom analysis package, based on Access and S-Plus software and custom-designed scripts (R.P., unpublished work). In a preliminary step, all raw datasets within a time group were plotted against each other and visually inspected. In most cases, reproducibility was good, except for two datasets that clearly departed from other replicates, likely because of contamination with exocrine tissue; these datasets were removed from further consideration. Probe-level normalization was performed with the RMA algorithm (20) from the Bioconductor site (www.bioconductor.org, implemented on the Insightful Array analyzer), which we found to yield data far less noisy than when processed with Affymetrix MAS-5 software (20) (R.P., J. Luckey, data not shown). For each time point, the set of normalized values were compared with the day 0 set, and a P value of differences was calculated by Welsh’s t test. An aggregate expression value for each feature of the Mu74Av2 chip and each time point was then calculated by simple averaging. Fold changes relative to the untreated day 0 sample were calculated from these mean values. The entire datasets are available as supplemental data in an online appendix (Table S1 [at http://diabetes.diabetesjournals.org]).

Two main strategies were used to assess the significance of the observations. First, since each experiment included a control pool, the study design generated more samples in the control day 0 group than for other time points. In validation analyses, this allowed us to divide the day 0 group into two individual subgroups (four samples each), providing an internal estimation of the intragroup experimental variability. As will be detailed below, a number of changes were observed at day 3 after CY treatment. Figure 1A shows that the differences between the day 3 and day 0 conditions were far more numerous than those between the replicate day 0 groups. This first analysis indicated that most of the variations described below were likely to be true and not merely a consequence of experimental fluctuation. Second, randomized data groups were generated by permutation between the 18 datasets of the expression values for each gene (row-wise permutation) and were then processed following the same methods as those for the true datasets. This generated matrices of fold change and P values for the null hypothesis (no time-dependent expression changes). This randomized dataset was used in unimodal comparisons of fold-change distributions (Table 2, see below). We also took advantage of the kinetic nature of the data. While the limited number of time points precluded the utilization of sophisticated time series techniques, an apparent false-positive rate (FPR) was estimated for the genes of interest by combining the probabilities of coordinated variation in successive time points. For each gene in the real dataset, we counted the number of genes in the random dataset presenting a similar fold-change pattern (day 2 and day 3 only). For the simulation shown in Fig. 1B, close to 100,000 random “gene patterns” were thus queried. The vast majority of the patterns were observed quite frequently (10⁴ times or more), but a few were very rare or absent from the randomly generated data. Not surprisingly, these corresponded to genes showing the greatest fold change at day 3 relative to day 0 (Fig. 1C). From these computations, we derived an estimate of the FPR for each gene that showed a time-dependent variation.

As a first analysis of changes in gene expression during the unfolding of CY-induced diabetes, we counted the overall number of genes that showed a change relative to the starting (day 0) value (Table 2). Overall, the effects of CY treatment mainly became apparent at day 3 after treatment, with a distribution of fold changes that was significantly different from that of the randomized data (Table 2, last column). For instance, an aggregate of 49 genes was induced or repressed by a factor of 3 or more, whereas only five genes showed such variation in the randomized dataset. At earlier times, while some of the changes may have been significant (see below), the overall pattern did not markedly depart from that of the randomized data. Interestingly, a roughly equivalent number of genes was induced or repressed by CY, with a fairly symmetric distribution in Table 2. This indicates that the cellular or genomic changes accompanying progression to diabetes were not limited to induction responses; it was not simply that an inflammatory response developed, but the process also included the downregulation of other functions.

We then analyzed in more detail the kinetics of expression of those genes in which transcripts were clearly altered by day 3 of CY treatment (filtered as described below). Figure 2 displays the distribution of expression at each day after CY treatment, plotted in reference to the day 0 values. Genes that eventually became induced by day 3 are shown in red, and those repressed are in blue. At day 1 (Fig. 2A), the majority of genes in the induced set were not yet affected, still clustering around the diagonal, except for a few repressed genes already showing a clear change. At day 2 (Fig. 2B), most genes in the set began to show clear differences from the day 0 values, albeit still less than the full induction seen at day 3. For the repressed gene set, a few were already underexpressed at day 1, with an amplification of the trend at day 2.

Fig. 3A presents the identity and principal expression values for the 88 genes of the “induced set.” These were selected with a cutoff on fold-change values of 2.1, selected as corresponding to a univariate FPR of 20% or better. (We did not use the t test P value as a filter because the estimate of the variance used in the t test is highly unreliable when dealing with a limited number of replicates, whereas the FPR derived from simulation makes no assumptions.) The t test P value and estimated FPR are also shown for reference in Fig. 3A. Many of these genes exhibited a partial induction at day 2, and the combined analysis in reference to the randomized control shows that their true significance is far greater (FPRs 2.10⁻⁴ to 5.10⁻⁶). This cutoff may eliminate from further examination a number of genes that exhibited significant, but more subtle, variations; however, they fall within a significantly more noisy range of the dataset. The following points can be made from the list in Fig. 3A:

The main cell type expressing each gene is indicated by color coding of the “locus link” column of the Fig., as determined by 1) comparing relative expression in Mu74Av2 chips probed with whole-spleen or islet RNA (L. Poirot, C.B., D.M., unpublished data), from which the islet-to-spleen ratio was a clear indicator of preferential expression in pancreatic or hematopoietic cells; and 2) retrieving information corresponding to the same chip from the GNF Gene Expression Atlas site. The majority of induced genes were typical of those transcribed by hematopoietic cells (myeloid or lymphoid, highlighted in green), and only a minority reflected responses in the target islet cells (highlighted in red). Most of these latter genes are represented by a group of genes that are members of the Reg family (Regα, Regγ, and Pap). These pancreas-specific C-type lectins are associated with conditions of pancreatic stress or active growth (e.g., pancreatitis, partial pancreatectomy), but their exact function and origin remains unknown (21–25). All members of the family were induced in a synchronous fashion: a very slight increase at day 1, a dip at day 2, and robust induction at day 3 (Fig. 3A). This coordinated response likely denotes a stress response by pancreatic cells, possibly involving cell proliferation in response to the initial noxious effects of the CY-induced decompensation, but also perhaps a reaction to the islet isolation procedure. Overall, however, it is clear that most of the action was taking place in the infiltrate because CY-induced diabetes was not initiated in the target cells themselves.

The color-coding of the “gene symbol” column of Fig. 3A, based on text searches through the cMu74Av2 chip annotation fields and PubMed databases, highlights those genes in which expression is controlled to some extent by IFN-γ. Twenty-five of the 88 genes in the induced set belong to this class (more may also be IFN-γ–responsive, but not yet recognized as such). To some extent, the earlier descriptions of inflammatory changes in CY-induced diabetes did foretell the upregulation of a range of inflammatory cytokines and downstream transcriptional programs (11,14,26,27). Yet, the current data highlight the fact that the effects were very much focused on IFN-γ. Much less prevalent were other inflammatory cytokines that one might have expected to find, such as tumor necrosis factor (TNF)-α or other members of the TNF family. There have been conflicting reports on the importance of IFN-γ in the development of diabetes in the NOD mouse, with different gene ablations of the cytokine or its receptor having contrasting effects in backcrossed animals (28–31). IFN-γ blockade by antibody treatment also had a limited effect on disease (11). The overwhelming impact of IFN-γ uncovered here argues for a reevaluation of this question.

These included Ccl2 and -7 and Cxcl9, -5, and -1. These chemokine transcripts showed parallel profiles of induction, i.e., slightly elevated at day 2, more fully induced at day 3 (Fig. 4A), but with a more robust induction for Ccl1 and -5 than for the other chemokine genes. Some of these genes are known to be IFN-γ responsive (but this is not true of all), and their kinetics of induction (contemporary with those of IFN-γ) suggest that other factors might have been responsible for their induction. Interestingly, these induced chemokine genes are clustered in murine genomes, particularly in two regions located on Chr5 and Chr11. Examination of the regions (Fig. 4B) shows that the transcriptional activation did not involve the entire cluster, i.e., only some members were induced, whereas others remained silent or unchanged (Figs. 3 and 4 B and online appendix). Thus, the activation of gene expression did not represent a wholesale activation of chromosomal regions, but rather quite specific effects on particular genes. The disposition of the induced chemokines follows an interesting pattern: in both cases, the most 5′ members of the cluster were those most significantly activated during CY-induced progression to diabetes, whereas there was much less effect on the more distal members of the locus. Is this disposition purely coincidental? It may reflect shared enhancer/response elements, although 130 kb separate Cxcl1 from Cxcl5, which is only 15 kb from the uninduced Cxcl2. It may also reflect the evolutionary history of the loci, inducibility being unequally conferred to duplicated members.

It is important to provide some independent validation of microarray data. In effect, several of the increases had already been validated, since the induction of IL-1α, IL-1β, IFN-γ, and IL-12 have previously been reported (11). In addition, we quantitated transcripts of several of the genes belonging to the induced set (Table 3). Transcripts of IFN-γ and the Cxcl9 chemokine (a.k.a. MIG) genes were measured by real-time quantitative analysis (TaqMan). For both genes, the increase detected on the arrays was mirrored in the real-time quantitative PCR analyses. As is commonly observed, the increase in gene expression was scored as greater in the PCR assay than on the microarrays, which tend to underrepresent the degree of variation; thus, the values in Fig. 3 are likely an underestimate of the true range of variation during progression to diabetes.

More generally, the experiments reported here present a broad perspective on gene inductions that accompany CY-induced diabetogenesis. The comparison with PCR data, and in particular with our prior analysis of cytokine induction after CY treatment (11), also illustrates the limitations of the microarray analysis. Several of the genes in which induction was anticipated were indeed members of the induced set of the present experiments (e.g., IFN-γ, IL-1α, and IL-1β). On the other hand, several others scored as unchanged (IL-12p35, TNF-α, and Fas; see online appendix). Most of these include genes with modest induction scores by PCR analysis, and it has been our experience that microarray analysis very often (but not always) underestimates the magnitude of induction; this is particularly true when the RMA technique is used for probe-level normalization. Others, such as IL-2, were simply not detected in either condition, likely because of the variable efficacy of oligonucleotide chips. Thus, oligonucleotide chips are powerful tools for revealing broad patterns of coordinate changes, but are of doubtful reliability for quantitating a specific gene product.

The 84 genes in which transcripts decreased in the 3 days following CY treatment are listed in Fig. 3B. They are very different in class and composition from the induced gene set.

Twenty-three of the repressed transcripts issued from immunoglobulin loci (heavy or light chains, variable or constant regions) and were among those that showed the most profound reduction, 8 of 10 of the transcripts showing the greatest reduction were from immunoglobulin genes. The reduction was not limited to immunoglobulin genes, a number of other B-cell transcripts were also affected, such as genes encoding molecules involved in antigen presentation (MHC class II, invariant chain), cell surface receptors (CD52, CD79, and CD83), and signal transduction molecules (CD5, Blk, and Dgkα). More than likely, this broad set of decreases reflects a drop in the proportion of B-cells in the insulitic lesion. Insulitis in BDC2.5/NOD mice (and in NOD animals as well) includes a very sizeable proportion of B-cells, which can account for up to 30% of the infiltrate area in immunohistochemical analyses (M.M., R.P., D.M., C.B., unpublished data). CY is known to be particularly toxic to B-cells; treatment with this agent provokes the death of B-cells within 24–48 h after administration (32). Thus, a wave of B-cell death within the infiltrate is likely to be the root of this strong reduction in B-cell–specific transcripts.

None of those corresponded to the primordial function of pancreatic endocrine cells; in particular, insulin transcripts were only reduced from 6,760 to 6,300 and 5,240 to 4,390 units for Ins1 and Ins2, respectively. This finding indicates that a generalized loss of β-cell function had not yet occurred by 3 days after CY treatment. The transcriptional changes that affected islet cells were predominantly centered on secretory or protease control functions (secretogranin, ChromograninA, Kexin2, and Spi family members), perhaps indicative of early changes in secretory pathways.

Given the suggestive genomic clustering of the induced chemokine genes, we analyzed more extensively the genomic distribution of those genes in which expression was affected by CY, clustering them according to their chromosomal position (information from the Affymetrix website [NatAffx.com], complemented by BLAST [basic local alignment search tool] searches on the Ensembl genome browser). Genes were considered to potentially belong to a cluster when they mapped <500 kb apart and were either induced or repressed after CY administration. Several interesting clusters were highlighted by this search (Fig. 5). Some were expected, such as the chemokine cluster described above or the immunoglobulin heavy and light chain clusters. The Reg genes also all map within a 140-kb stretch. Others were somewhat of a surprise and may underscore hitherto unrecognized relationships; for example, the “prostate tumor overexpressed gene 1” and “’kidney-derived aspartic protease-like protein” were coordinately repressed and map close to each other on Chr7. Conversely, the genes encoding Pbef and Gdap10 (pre–B-cell colony-enhancing factor and ganglioside-induced differentiation associated protein 10, respectively) map within 20 kb of each other on Chr12 and were induced in concert. Only in one instance did genes in the same geographical cluster respond discordantly, underscoring the significance of the observations, given that random clustering would be expected to generate equal numbers of concordant and discordant clusters.

These results paint a broad portrait of the molecular processes participating in progression to diabetes subsequent to treatment of insulitic mice with CY. This progression involves the extinction of B-cell transcripts and induction of a response dominated by IFN-γ and the genes it controls. Interestingly, this pattern has significant but incomplete overlap with that detected when comparing the “respectful” and “destructive” islet infiltrates that develop when the BDC2.5TCR transgene is bred onto different genetic backgrounds (L. Poirot, D.M., C.B., unpublished observations). Overabundant IFN-γ gene expression and overexpression of certain IFN-γ–responsive transcripts were also found in the latter context, whereas other genes highlighted in the present study were not differentially expressed, including the chemokine and Schlafen gene families or the B-cell–derived transcripts. Thus, the various aggressive/destructive forms of insulitis observed under different conditions (CY induced or genetic background influenced) may represent a common cellular phenotype at which different immunoregulatory perturbations converge. IFN-γ and the genes it induces would be at the point of convergence, whether resulting from NK-related overexpression or from chemical induction. IFN-γ would then activate the full cytotoxic potential of the various cell types present in the lesion, i.e., T-cells, but also NK cells or macrophages, or even neutrophils, which are secondarily recruited after CY induction (33).

In the case of CY-induced diabetes, one might speculate that the disappearance of B-cells is at the root of the process, either because B-cell depletion results in the release of large amounts of apoptotic cell fragments, thus activating the myeloid cells that abound in the lesion and altering their capacity to interact with/control T-cells, or because B-cells exert a regulatory influence in and of themselves, as was proposed by Turk (34) in the mid-1970s and more recently in a murine model of arthritis (35). Disappearance of B-cells would unleash, then, aggressive insulitis. Although not usually thought as central to the pathogenesis of type 1 diabetes, B-cells do play an important role, as shown by the quasi-absence of autoimmune attacks in B-cell–deficient NOD mice (36–38), although this point is more controversial in the human milieu (39). This is usually interpreted as a perturbation of antigen presentation, which prevents the initial priming of autoreactive T-cells. The present data lead one to speculate that B-cells may also be involved in later stages of the process, when regulatory balances condition the outcome of the autoimmune lesion.

On the other hand, our data do not support the oft-evoked notion that CY exerts its influence by depleting suppressor T-cells. The aggressiveness of the insulitic lesion in BDC2.5 mice is strongly conditioned by regulatory T-cell populations, which are under the control of costimulatory family genes such as ICOS or CTLA-4 (6,40) (A. Herman, C.B., D.M., unpublished observations). In the latter study, diabetes induced by ICOS blockade correlated with clear changes in the “Treg signature,” a set of 150 genes differentially expressed in CD25⁺CD4⁺ regulatory cells (41). In the present data, such changes were conspicuously absent. In addition, we have not observed significant reductions in pancreatic CD25⁺ T regulatory cells after CY treatment (A. Herman, D.M., C.B., unpublished data). Thus, and unless CY affects a very minor subpopulation with regulatory potential, changes in T regulatory populations do not seem prevalent in CY-induced diabetes.

In conclusion, this analysis reveals a landscape dominated, to a striking degree, by IFN-γ and the genes it controls. The model used here certainly represents an extreme form of β-cell destruction, and it will be important to assess whether IFN-γ occupies such a central position in other instances of autoimmune β-cell destruction, and whether the molecular and cellular pathways uncovered in this mouse system play a role in human patients.

This work was supported by grants from the National Institutes of Health and the JDRF to D.M. and C.B. (P01 AI 39671 and 9-1998-1004), by grants from Joslin’s National Institute of Diabetes and Digestive and Kidney Diseases (NIDDK)-funded Diabetes and Endocrinology Research Center (DERC) cores, and by funds from the W.T. Young Chair for Diabetes Research. M.M. was supported by grants from the JDRF and the Endocrine Fellows Foundation.

We thank members of the Benoist-Mathis lab and John Rogus for discussion; Jack O’Neill, Vaja Tchipashvili, and Gordon Weir of the JDRF Center for Islet Transplantation at Harvard’s Islet Preparation Core for provision of islets; Robert Saccone of the DERC Microarray Core for processing and hybridization of the probes; and Judy George and Ella Hyatt for maintaining the BDC2.5 colony.