extfieldsofstudy: Biology Medicine

s2fieldsofstudy: Biology Medicine

text:
Subgroup-specific immune and stromal microenvironment in medulloblastoma

ABSTRACT Knowledge on immune and stromal cells in medulloblastoma microenvironment is still limited as previous work was frequently restricted by low sample size and the lack of molecular subgroup information. We characterized 10 microenvironment cell populations as well as PD-L1 from gene expression in 1422 brain tumors and 763 medulloblastomas. All in all, medulloblastomas showed low expression of immune markers. Still, there were substantial differences with a clustering of medulloblastoma subgroups according to their microenvironment profile. Specifically, SHH medulloblastomas displayed strong signatures of fibroblasts, T cells and macrophages, while markers of cytotoxic lymphocytes were enriched in Group 4 tumors. PD-L1 gene expression appeared to be relatively high in single SHH and WNT cases but was undetectable by immunohistochemistry. In addition, two diverse immuno-stromal patterns were identified, indicating distinct types of local tumor immunosuppression, which were primarily controlled by either macrophage and regulatory T cell-mediated mechanisms or immunosuppressive cytokines and checkpoints, respectively. None of the immune cell signatures had an independent prognostic value in the present dataset after multiple testing correction. These results suggest a mild, but subgroup-specific infiltration of immune cells in medulloblastoma.

Introduction
Medulloblastoma is the most frequent malignant brain tumor in children. A comprehensive molecular characterization generating large amounts of genomic, epigenomic, and transcriptomic data has permitted a precise profiling of this entity and the definition of four major molecular subgroups. [1][2][3][4][5] These subgroups and additional genetic alterations as well as clinical factors allow for a precise risk stratification of medulloblastoma patients. 6 The tumor microenvironment is of major importance for tumor growth and progression as well as treatment response in many cancers. [7][8][9] It has further gained interest in the last years with the availability of immune checkpoint inhibitors, which allow for a therapeutic targeting of the interaction between tumor and immune system. 10 Tumor-infiltrating immune cells, principally lymphocytes, have been shown to be of prognostic relevance and predictive for response to chemotherapy in various solid tumors. 11,12 In particular, extensive research has been performed on tumor immunology in glioblastoma 13,14 and the potential for immunotherapy has been discussed since many years in medulloblastoma. 15 However, despite a very deep molecular characterization of this entity, immunotherapeutic approaches are still rare, which is likely due to limited knowledge about the tumor microenvironment and the prognostic relevance of the cells contributing to it.
Using multicolor FACS analysis in a cohort of 6 medulloblastomas, it was shown that medulloblastoma exhibited a different immune microenvironment compared to 3 other pediatric brain tumor entities and the findings were validated by gene expression. 16 Natural killer cell (NK cell) therapy was among the first immunotherapeutic approaches that have been proposed for medulloblastoma (see 17 for review). Further, Margol et al. 18 showed that tumor-associated macrophages were more numerous in SHH medulloblastoma than in Group 3 and Group 4 tumors and proposed macrophages as therapeutic targets.
With the avenue of immune checkpoint inhibitors, the PD-1/PD-L1 axis started being studied in medulloblastoma. Pham et al. 19 developed murine SHH and Group 3 models and showed that immune infiltrates were more present in SHH tumors. However, there were more PD-1C CD8C T cells in the Group 3 animals. Interestingly, treatment with PD-L1 and CTLA4 inhibitors provided a survival benefit in the Group 3 animals only. In a cohort including 16 tumors, Murata et al. reported that high expression of PD-L1 (seen in 56% of the cases) was associated with low CD8C T cell infiltration and poor prognosis in medulloblastoma. 20 In contrast, Majzner et al. 21 (40 cases) and Vermeulen et al. 22 (29 cases) did not see PD-L1 expression in any of the studied medulloblastomas. None of the last three studies found subgroup-specific patterns of tumor-infiltrating lymphocytes (TILs).
Most of these results are mainly based on preclinical models or on immunohistological profiling of tumors. This yields to limitations in sample size and reduced statistical power as well as small numbers of markers that can be studied simultaneously. Furthermore, some of the immunohistological results on the expression of PD-L1 in medulloblastoma were contradictory. Associations with survival in larger cohorts have not yet been published.
Recently, large gene expression datasets of medulloblastoma have been released, including one dataset of 763 tumors. 4,5 Also, computational methods have been proposed to characterize the cellular composition of the tumor microenvironment from gene expression data. One of the main methods, CIBER-SORT, 23 uses a deconvolution algorithm that allows for the quantification of 22 immune subpopulations. The algorithm was primarily designed to infer the relative proportion of immune cells and successfully established different immune cell populations as prognostic markers in various cancers. 24,25 Furthermore, the recently developed Microenvironment Cell Populations-counter 26 (MCPcounter) uses characteristic gene signatures that are specifically associated with 8 immune cell populations (T cells, CD8C T cells, cytotoxic lymphocytes, B lineage cells, NK cells, monocytic lineage cells, myeloid dendritic cells, neutrophils) and 2 stromal cell populations (fibroblasts and endothelial cells) and allows for their absolute quantification.
Hence, we set out to analyze the tumor microenvironment of medulloblastoma using gene expression data. This allows to overcome sample size limitations and to study several immune cell populations simultaneously. We used MCPcounter as it also quantifies two stromal cell populations and yields more consistent results if there are only few infiltrating immune cells, which is frequently the case for pediatric brain tumors. We first compared medulloblastoma to seven other major brain cancers. Further, we used a large published gene expression dataset to infer subgroup-specific differences and associations with survival.

Results
Medulloblastoma shows low expression of PD-L1 and microenvironment cell signatures in comparison with seven other brain tumors Although a more precise division into existing molecular tumor subgroups was not possible due the lack of respective information, this should give a rough overview on microenvironment signatures in MB compared to other brain tumors.
The expression scores for all cell populations and PD-L1 were highly significantly different between entities (all p <5.4e-8, see supplement). The overall expression of immune and stromal markers in medulloblastoma was low compared to the other entities (Fig. 1A, B). In the clustering analysis, the expression patterns of medulloblastoma were most similar to the embryonal AT/RT and ependymomas (Fig. 1A). The astrocytic tumors built one cluster and were the only ones to show high expression levels of PD-L1 and markers for cells from the monocytic lineage, e.g. macrophages. There were lower levels of cytotoxic lymphocytes but higher levels of neutrophils in GBM than in the lower grade gliomas. MGM showed a distinct expression pattern to the other tumor entities with large numbers of fibroblasts as well as myeloid dendritic cells and NK cells.
Overall, medulloblastomas had the lowest gene expression level of PD-L1, cytotoxic lymphocytes, and endothelial cells of all studied entities (Fig. 1B). However, there was substantial variation among different medulloblastomas for PD-L1 and the microenvironment cell populations. The expression levels were not uniformly low, indeed, single cases of medulloblastoma showed higher expression than the median of all tumors for all MCPcounter populations and PD-L1.
Hierarchical clustering of all samples revealed a clear-cut aggregation of meningeal tumors, but the remaining entities did not form sharply separated groups (Suppl. Fig. 1). In particular, medulloblastomas did not form a distinct cluster but were mixed with the other entities. Hence, although most medulloblastomas showed lower levels of immune infiltration and stromal cells than the remaining brain tumors some tumors exhibited higher expression of certain cell populations similar to brain tumors from other entities.
Expression of PD-L1 and microenvironment signatures is highly subgroup-specific Next, we studied the microenvironment populations (MCP) and PD-L1 expression in a gene expression dataset of 763 medulloblastomas with known molecular subgroup (70 WNT, 223 SHH, 144 Group 3 and 326 Group 4). 5 The clustering analysis of the 10 MCPs separates most of the SHH and WNT tumors from most of the Group 3 and Group 4 tumors ( Fig. 2A). There is one cluster composed mostly of SHH and WNT tumor that shows relatively high PD-L1 expression.
The expression scores were strongly significantly different between subgroups for PD-L1 and all MCPs with the exception of the B lineage and NK cells (all p <3.85e-10, Fig. 2B). Median PD-L1 expression was highest in WNT tumors (4.06, log 2 -scale) followed by SHH tumors (3.67) and lowest in Group 4 tumors (3.54, fold change WNT/ Group 4: 1.4). The distribution of PD-L1 expression was heavy-tailed and few tumors mostly from the SHH and the WNT subgroup showed comparably high expression ( Fig. 2A). Indeed, 20% of WNT, 6.7% of SHH, 2.7% of Group 3, and 0.3% of Group 4 tumors had expression of PD-L1 higher than 4.5, while 5.7% of WNT, 2.7% of SHH, 1.4% of Group 3 and 0.3% of Group 4 tumors had expression higher than 5.0. Still, in a series of 10 WNT and 36 SHH medulloblastomas we detected small islets of perivascular PD-L1 positive cells, but none of the samples showed membranous PD-L1 expression in more than 1% of tumors cells by using antibodies against PD-L1 (data not shown).
SHH tumors were associated with the most distinctive pattern of microenvironmental cells (Fig. 2B). In particular, they had larger numbers of T cells, fibroblasts, and cells from the monocytic lineage and lower numbers of neutrophils than the tumors from the other subgroups. Group 3 tumors had low numbers of T cells and macrophages, but the highest number of CD8C T cells. Group 4 tumors had the lowest number of fibroblasts and largest number of cytotoxic lymphocytes.

Immuno-stromal classification of medulloblastoma
To get a more precise overview of the complex relations between molecular subgroups and immuno-stromal context, we computed a two-dimensional embedding of 98 MCP markers using t-SNE 27 yielding an almost complete separation of the molecular subgroups (Fig. 3A, Suppl. Tables, Suppl. Fig. 2). Using unsupervised k-means clustering, we identified 5 major clusters with distinct immuno-stromal The bands indicate the median, the lower and the upper quartile. Medulloblastoma shows rather low expression scores for most cell types and PD-L1 compared to other brain tumors but the variance is large. microenvironment (Fig. 3B). Furthermore, four previously published signatures 28 were used to identify immunosuppressive mechanisms (see methods and supplement). Individual clustering of distinct signature sets revealed two major immuno-stromal patterns (Fig. 3C). Type I pattern represented a primarily monocytic/macrophage supported tumor environment, with a major fibroblastic stromal component and infiltration of regulatory T cells. Type II pattern on the other hand was defined by a stronger infiltration of cytotoxic CD8C T cells and neutrophils with a simultaneous activation of immunosuppressive cytokine and checkpoint pathways, which include CD274 (PD-L1), CTLA4 and TGFB1. Endothelial-, B-, NK-and myeloid dendritic cells did not group into these two patterns (Fig. 3C, lower panel). The type I pattern, found in cluster 2 and 4 ( Fig. 3B), can be assigned to SHH-, WNT-and approximately one quarter of Group 4 medulloblastoma. Type II pattern was representative for Group 3 and the remaining Group 4 tumors. Interestingly, the unsupervised k-means clustering demonstrates that the two distinct immuno-stromal patterns separate Group 4 tumors into two subgroups (cluster 1 and 4), indicating two potentially diverse mechanisms of immunosuppression within the microenvironment of Group 4 medulloblastoma.

Associations of microenvironment signatures with survival
The dataset was split into subgroups prior to survival analysis since molecular subgroups are associated with a significantly different prognosis. The overall association of immune cells with survival was only weak. There were no significant associations with survival in WNT tumors, which might be explained by the good prognosis associated with a low number of events and low statistical power (Fig. 4). The expression score for fibroblasts in SHH tumors in univariate analysis was the only significant finding after multiple testing correction using the Benjamini-Hochberg (BH) method (p D 0.003, BH-p D 0.033, Fig. 4). However, in multivariate analyses adjusted for age, histology, metastatic status, and gender, it narrowly missed significance (p D 0.078). There were no significant associations between MCP-markers or PD-L1 in Group 3 and Group 4 tumors after multiple testing correction. In univariate analysis of Group 3 and Group 4 tumors, myeloid dendritic cells marker expression was associated with better survival (p D 0.028, BH-p D 0.32 and p D 0.05, BH-p D 0.36). In multivariate analysis, this association also lost significance (p D 0.32 and p D 0.41, respectively). There were no significant associations between survival and the expression of the markers for the monocytic lineage (macrophages) or NK cells in any of the subgroups.

Discussion
Medulloblastoma, an embryonal tumor originating from the posterior fossa, is the most common malignant pediatric brain tumor. With multimodal therapeutic concepts including surgery, radiotherapy, and chemotherapy, about two thirds of the patients can be cured. However, the prognosis of very high-risk patients is still poor. 6 Further, many patients suffer from important treatment associated neurocognitive sequelae  28 There are two major patterns (type I C II) potentially associated with immunosuppression in medulloblastoma. Endothelial-, B-, NK-and myeloid dendritic cells did not group into these two patterns.
(reviewed in 29 ). Hence, there is a great need for the development of new therapeutic strategies, and immunotherapy, which has been successfully used in the treatment of many adult solid tumors, is a promising treatment option. A precise understanding of the function of the tumor microenvironment appears key to design medulloblastoma specific immunotherapies. However, there is still not much known about the tumor microenvironment of medulloblastoma. Presumably due to the rarity of medulloblastoma compared to major adult cancers, there are no large immunohistochemical studies characterizing the microenvironment of this entity. Previous reports were limited by low sample sizes and/or focused on only one or two cell populations. 16,[18][19][20][21][22] Hence, to gain a more comprehensive overview, we used a state of the art bioinformatics method to quantify eight immune cell populations and two stromal cell populations as well as PD-L1 in a pancancer dataset of 1422 brain tumors and a medulloblastoma dataset including 763 samples. This allows for the simultaneous characterization of several cell populations and overcomes the described limitations. To our knowledge, this is the largest and most comprehensive study of the microenvironment in medulloblastoma.
The gene signatures proposed by Becht et al. have been carefully selected to be specific for the described immune cell populations and validated in several thousand tumors including a small dataset of medulloblastoma. 26 A study of CD8 and Th17 infiltrates in GBM 30 showed comparable associations with prognosis of immunohistological data and data based on gene expression signatures in a similar approach to those proposed by Becht et al. Further, intratumoral heterogeneity of gene expression in medulloblastoma is very low, 31 which also supports our assumption that transcriptome based analysis is suitable for a global description of the tumor microenvironment. However, it cannot be excluded that tumor cells contribute to the expression score, which might significantly confound the analysis. Indeed, the expression of certain genes might be associated with different cell types in different tissues/organs. Thus, we reviewed all the marker genes manually and excluded genes associated with medulloblastoma tumor biology (see supplement and methods). For example, EOMES is highly expressed in cytotoxic lymphocytes but only weakly in CD4C T cells, 32 making it a suitable marker for this cell population outside of the nervous system. However, it is also associated with embryonal brain development and a specific marker for Group 4 medulloblastoma, 33 which would clearly confound the deconvolution approach. This shows that a careful and tissue-specific selection of the markers is key and should also be performed in further studies.
There have been large numbers of immunohistochemical studies of immune infiltrates particularly in glioma and meningioma including T cells, NK cells and macrophages, with partially contradictory results (reviewed in [34][35][36]. This makes it difficult to establish a clear reference level for immune infiltration in these entities, however, for most of the microenvironment populations used in this study, the expression levels are compatible with previously published results based on immunohistochemistry. In the pancancer analysis, PD-L1 expression was highest in GBM in line with results based on immunohistology. 37 The very distinctive profile of meningiomas with a larger number of fibroblasts and a larger number of NK cells might be linked to the different histological origin of these tumors. Medulloblastoma had the lowest amount of PD-L1 and cytotoxic lymphocytes of all studied entities and overall a very small amount of infiltrating immune cells. This suggests that the tumor either actively modulates the immune response or simply has little immunogenicity as is suggested by the relatively low mutational burden. 38 However, a recent study also demonstrated some form of systemic immune suppression present in medulloblastoma with lower lymphocyte counts and a higher neutrophil-to-lymphocyte count ratio already present at the moment of diagnosis prior to the application of steroids. 39 Our study demonstrated highly significant differences in the pattern of immune and stromal cell signatures as well as PD-L1 between the molecular subgroups of medulloblastoma. Margol et al. 18 were the first to highlight differences in the microenvironment of medulloblastoma subgroups and showed that there were more macrophages in SHH tumors. They also defined a gene expression signature with 25 tumor genes as well as 6 immune genes able to accurately classify medulloblastoma into subgroups. Here, we showed that these differences apply to all immune and stromal cell populations except NK cells and B cells. Further, our analysis shows that an immune signature alone without genes associated with tumor biology yields a good separation of medulloblastoma subgroups in a clustering analysis (Fig. 3, Suppl. Fig. 2). Recently, it could be shown that WNT medulloblastoma but not the other subgroups lack a functional blood-brain barrier. 40 This might explain differences between WNT tumors and the remaining subgroups, but cannot explain the relatively distinct immune profile of SHH tumors.
Since SHH and WNT tumors had the highest expression levels of PD-L1, we also stained WNT and SHH tumors for PD-L1. Although this is the largest series of medulloblastomas stained for PD-L1 currently reported, we did not find membranous PD-L1 expression in more than 1% of tumor cells in any of the studied samples. This confirms previous studies, 21,22 which did not see PD-L1 expression in medulloblastoma and contradicts a small study by Murata et al. 20 reporting PD-L1 expression in 9/16 samples. Hence, immunohistological staining of PD-L1 cannot be used as a predictive biomarker for checkpoint inhibition. However, results of PD-L1 immunohistochemistry are variable and patients exist without detectable PD-L1 expression in tumor cells showing a substantial response to PD-L1 inhibitory therapy. 39 In a murine model of medulloblastoma, disruption of CDK5 was associated with lower PD-L1 expression and tumor rejection mediated by CD4 T cells, indicating biological relevance of the PD-L1 pathway in medulloblastoma. 41 High levels of interferon gamma (IFNg) were present in CDK5 deficient tumors. 41 As we previously reported, TH1 predominance with high expression of IFNg in the peripheral blood was associated with favorable outcome in children with medulloblastoma treated by high-dose chemotherapy followed by autologous stem cell transplantation. 42 This further corroborates the relevance of the IFNg/PD-L1 axis in medulloblastoma although its precise functioning in vivo remains to be elucidated.
SHH tumors exhibited larger numbers of fibroblasts, macrophages and T cells than the other subgroups. The fibroblast score is mostly based on the expression of collagens and other elements from the extracellular matrix, and it cannot be excluded that these molecules are produced by tumor cells. Higher numbers of fibroblasts were associated with significantly better prognosis and might be linked to a better differentiation of these tumors. The larger amount of fibroblasts in SHH tumors may be explained by desmoplastic histology, which is known to be associated with better prognosis 43 and the addition of which as a marker reduces the significance of fibroblast expression in multivariate analysis. Nonetheless, the numerical amount of fibroblasts seems to be an important factor.
We identified two distinct immuno-stromal patterns, defined by an either monocytic/macrophage supported pro-tumorigenic microenvironment versus a pattern which suppresses cytotoxic T cell infiltration by immunosuppressive cytokines and checkpoint pathways. Interestingly, as the diverse patterns, with the exception of the divided Group 4 tumors, are accurately distributed to the known medulloblastoma subtypes, a presumable relation to driver pathways has to be considered. Recent studies support such a pathway-specific immunomodulation. For example, an in vivo colitis study showed that a specific activation of the hedgehog pathway resulted in local immunosuppression mediated by stromal and regulatory T cells. 44 A similar immunomodulation of the tumor microenvironment was described in WNT/b-catenin pathway activation. Here, b-catenin activation supported tumor-associated macrophage and regulatory T cell survival (summarized in 45 ). Furthermore, segregation of TCGA datasets according to an inflamed T cell signature set showed that WNT/b-catenin pathway activation inversely correlated with T cell infiltration. 46 These observations are in line with our finding of a predominantly monocytic/macrophage and Treg supported (type I) immunosuppressive environment within SHH and WNT medulloblastoma. A pathway-specific inhibition in these tumors could therefore result in a pleiotropic effect, simultaneously promoting tumor-specific immune activation. In contrast, according to our data, pattern type II tumors, including most group 3 and 4 medulloblastoma, will presumably benefit to a greater extent from a targeted inhibition of immunosuppressive signaling pathways, which is also supported by preclinical results from murine models. 19 Our analysis showed only weak associations between intratumoral immune infiltrates and survival, which all missed significance after multiple testing correction despite a large sample size. This might be explained by the overall very low immune infiltration, with small differences in the absolute levels of immune infiltrating cells. It corroborates findings by Vermeulen et al., who did not see associations between TILs and prognosis. 22 In contrast, peripheral levels of TH1 cells were associated with better prognosis already in a small cohort of medulloblastoma. 42 Hence, a simultaneous investigation of peripheral and intratumoral immune cells seems to be of key relevance in further studies of medulloblastoma tumor immunology.
In summary, we report that medulloblastoma has the lowest level of PD-L1 and cytotoxic lymphocytes of eight studied brain tumors. There are substantial differences in the microenvironment cell populations between molecular subgroups, which are likely to be associated with different interactions between the tumor and the immune system. In particular, two immunostromal patterns (macrophage and regulatory T cell-mediated mechanisms vs. immunosuppressive cytokines and checkpoints) were identified. Both patterns were closely linked to the molecular subgroups. Upcoming immunotherapeutic treatment strategies will have to consider these distinct microenvironmental phenotypes in order to efficiently promote a tumorspecific immune response in the context of different mechanisms of immune escape in medulloblastoma subgroups.

Material and methods
All data analyses were performed using the statistical programming language R 47 including the packages survival, gplots, bee-swarm and Rtsne. [48][49][50][51] Datasets and preprocessing Raw gene expression data (.CEL files) were downloaded from Gene Expression Omnibus (GEO). 52 For the pancancer analysis, 18 datasets covering 8 different tumors entities measured on the Affymetrix Human Genome U133 Plus 2.0 Array were used (see supplement). To allow for comparisons between different series, fRMA normalization, which includes a batch effect correction, was used as implemented in the R-package fRMA 53 using the input vectors from the R-package hgu133plus2frmavecs. 54 Duplicated samples (defined by a Pearson correlation > 0.995) were removed. This yielded a final dataset of 1422 unique tumor samples (421 glioblastomas (GBM), 282 medulloblastomas (MB), 272 ependymomas (EPN), 99 meningiomas (MGM), 159 pilocytic astrocytomas (PA), 61 anaplastic oligodendrogliomas (AO), 69 anaplastic astrocytomas (AA) and 59 atypical teratoid/rhabdoid tumors (AT/RT). To exclude major remaining batch effects after preprocessing, we compared samples from entities present in multiple gene expression series. Clustering analysis showed that MCPcounter profiles for the same entity were similar with substantial variance between entities (data not shown).
The data published by Cavalli et al. was used for medulloblastoma specific analyses. 5 Raw data were normalized using the RMA method as implemented in the R-package affy 55 with the custom chip definition file hugene11sthsensgcdf (v22.0.0). 56

Computation of immune and stromal signatures
Published signatures of mRNA markers from 10 cell populations (MCPcounter) were used to estimate the content of T cells, CD8C T cells, cytotoxic lymphocytes, B lineage cells, NK cells, monocytic lineage cells, myeloid dendritic cells, neutrophils, fibroblasts and endothelial cells in tumor tissue. 26 All markers were manually reviewed and those associated with medulloblastoma tumor biology were removed to avoid bias (see supplement). MCPcounter scores were defined as the average of the log 2 -scaled values of the corresponding probes or genes (see supplement) for the U133 Plus 2.0 Array and the 1.1 ST array respectively. Expression of PD-L1 was defined as the mean of the two corresponding probes (223834_at and 227458_at) for the U133 Plus 2.0 Array. The resulting scores are expressed in arbitrary units. Immunosuppressive signatures were obtained from Doucette et al. 28 and computed analogously.
Heatmap, differential expression, k-means clustering, t-SNE and survival analysis Medulloblastoma data was trimmed at Z-scores § 4 prior to heatmap analysis. Heatmaps were generated using the average linkage method and Pearson correlations as similarity measure. Significance of differential expression between subgroups of medulloblastoma was assessed using the Kruskal-Wallis test (two-sided).
Immuno-stromal clusters were defined using k-means clustering of the MCPcounter marker genes. The number of clusters (5) was selected using the elbow method. For visualization, t-SNE 27,51 was used with perplexity 100.
Survival analysis was performed using the proportional hazards model as implemented by the coxph function from the R package survival. For multivariate analysis, age, metastatic status, histology and gender were included. Samples with missing clinical data were excluded. Expression values for the different signatures were transformed to 3 quantiles (tertiles) and included as numerical values. Multiple testing correction for the cell populations and PD-L1 was performed using the Benjamini-Hochberg (BH) method where appropriated.

Immunohistochemistry
Immunohistochemical stains were performed on formalinfixed, paraffin-embedded tumor material according to standard protocols on an automated Ventana HX IHC system (Ventana-Roche Medical systems, Tucson, AZ, USA) following the manufacturer's instructions. Diaminobenzidine (Sigma) was used as Chromogen and sections were counterstained with hematoxylin. Antibodies against PD-L1 (E1L3N, Cell Signaling) were used in a dilution of 1:200.
Informed consent was obtained from all patients prior to the analysis.

Disclosure of potential conflicts of interest
No potential conflicts of interest were disclosed.

Funding
This research was supported by the F€ ordergemeinschaft Kinderkrebs-Zentrum Hamburg.


text:
Extraordinary Genome Instability and Widespread Chromosome Rearrangements During Vegetative Growth

The ability to rapidly adapt to changing environments is crucial for the success of pathogens infecting plants and animals. In some eukaryotic pathogens, rapid evolution can be mediated by genome rearrangements, which...

and animals have been shown to rapidly adapt to changing environments, including industrial agriculture and medical treatments, usually in response to adaptive pressure caused by the widespread use of fungicides and pharmaceuticals (Selmecki et al. 2010;Bennett et al. 2014). The genomes of many prominent fungal plant pathogens exhibit high levels of structural variation, including isolate-or lineage-specific regions often characterized by high repeat contents, and accessory chromosomes varying in frequency among individual isolates (de Jonge et al. 2013;Persoons et al. 2014;Plissonneau et al. 2016). Structural variation is often associated with meiotic recombination, but highly variable regions are also found in asexually and mostly asexually reproducing species (Wang et al. 2006;Chuma et al. 2011;Faino et al. 2016). In several plant pathogenic fungi, determinants of virulence, so-called effector proteins, are located in dynamic regions of the genome (Miao et al. 1991;Coleman et al. 2009;Ma et al. 2010). Accessory chromosomes often carry genes that encode virulence determinants and are linked to pathogenicity. Therefore, the absence of a particular chromosome in some species results in avirulent phenotypes (Ma et al. 2010;Tsuge et al. 2016;Vlaardingerbroek et al. 2016;van Dam et al. 2018). Little is known about the mechanisms that drive the dynamics of accessory chromosomes and highly variable regions in genomes of eukaryotic pathogens.
The plant pathogenic fungus Zymoseptoria tritici reproduces both sexually and asexually and causes disease on wheat, especially in Northern Europe and North America, resulting in extensive annual yield loss (Fones and Gurr 2015;Torriani et al. 2015). The Z. tritici reference isolate IPO323 contains eight accessory chromosomes with sizes ranging from 0.4 to 1 Mb, accounting for 12% of the entire genome . Chromosome rearrangements, including complete loss of accessory chromosomes, are a frequent phenomenon in this fungus and have been demonstrated to occur during meiosis (Wittenberg et al. 2009;Croll et al. 2015;Fouché et al. 2018). Infection experiments with Z. tritici strains in which single or multiple accessory chromosomes had been deleted showed that at least some accessory chromosomes encode virulence factors that determine host specificity (Habig et al. 2017). In contrast to the core chromosomes, accessory chromosomes are enriched with transposable elements and have low gene density . Chromatin of accessory chromosomes shows hallmarks of constitutive and facultative heterochromatin, consistent with the observed transcriptional silencing of genes present on these chromosomes (Kellner et al. 2014;Rudd et al. 2015;Schotanus et al. 2015). The centromeres, subtelomeric regions, and telomeric repeats of accessory chromosomes are indistinguishable from those of core chromosomes (Schotanus et al. 2015), suggesting that accessory chromosomes contain all the required regions for proper chromosome segregation. The unusual high number of accessory chromosomes, in combination with the extreme variability in chromosome content among field isolates or progeny from controlled crosses, make Z. tritici an excellent eukaryotic model to study accessory chromosomes and their dynamics (Mehrabi et al. 2007;Wittenberg et al. 2009;Stukenbrock et al. 2010;Goodwin et al. 2011). Here, we provide evidence for unexpectedly high rates of chromosome loss and changes during asexual propagation, both in vitro and in planta. We describe the types of structural rearrangements and overall variation in chromosome stability in this important crop pathogen.

Short-term in vitro growth experiment in liquid culture
Zymoseptoria strains were diluted from glycerol stocks (280°), plated on yeast-malt-sucrose (YMS) agar (4 g yeast extract, 4 g malt extract, 4 g sucrose, and 20 g agar per 1 liter) plates and grown for 7 days at 18°to obtain single colonies. One single colony was picked and suspended in 100 ml of YMS. Three replicate cultures were inoculated with 20 ml of cells from the single colony (50,000 cells). Cells were grown in 25 ml YMS medium at 18 or 28°shaking at 200 rpm and 900 ml of the cultures were transferred to fresh medium after 3-4 days of growth. In total eight transfers were conducted, for a total time course of 4 weeks.

Short-term in vitro growth experiment on agar plates
A single colony derived directly from a plated dilution of frozen stock for Zt09 (IPO323D18) was resuspended in 1000 ml YMS including 25% glycerol by 2 min vortexing on a VXR basic Vibrax at 2000 rpm, and 10-50 ml were replated onto a YMS agar plate. Forty replicates were produced. Cells were grown for 7 days at 18°whereby a random colony (based on vicinity to a prefixed position on the plate) derived from a single cell was picked and transferred to a new plate as described above. The transfer was conducted for a total of four times before a randomly chosen colony of each replicate was PCR screened and their complement on accessory chromosomes characterized as described below.

Screening by PCR
Cells from transfer eight of liquid cultures were diluted and plated on YMS-agar plates to obtain single colonies. To extract DNA, single colonies were suspended in 50 ml of 25 mM NaOH and boiled at 98°for 10 min; 50 ml of 40 mM Tris-HCl (pH 5.5) were added and 4 ml were used as template for the PCR. Primers for the right and left subtelomeric regions and close to the centromere were used for the chromosome-loss screening in our Z. tritici isolates, for Zymoseptoria ardabiliae primers in the center of candidate accessory chromosome unitigs were used. Primers and expected fragment lengths are listed in Supplemental Material, Table S1. Primers were designed with Clonemanager (Sci-Ed Software, Denver, CO) and MacVector and ordered from eurofins Genomics (Ebersberg, Germany).

Plant experiments and pycnidia isolation
Growth conditions for plants (Triticum aestivum) were 16 hr at light intensity of 200 mmol/m 22 s 21 and 8 hr darkness in growth chambers at 20°with 90% humidity. Seeds of the Z. tritici susceptible wheat cultivar Obelisk (Wiersum Plantbreeding BV, Winschoten, The Netherlands) were germinated on wet sterile Whatman paper for 4 days at growth conditions before potting. After potting, plants were further grown for 7 days before infection.
For the phenotypic analyses, three independent experiments were conducted and 21 plants were used per strain per experiment. An 5 cm long section on the second leaf of each plant was infected by brushing a cell suspension with 10 7 cell/ml in H 2 O + 0.1% Tween 20 on the abaxial and adaxial side of the second leaf. Plants were placed into sealed bags containing 1 liter of water for 48 hr to facilitate infection at maximum air humidity. Infected leaf areas were harvested 21 days postinfection for further phenotypic analysis or prepared for pycnidia isolation by surface sterilization with 1.2% NaClO for 2 min, followed by 70% ethanol for a few seconds, and washed twice with H 2 O. Leaves were placed into a sterile environment with maximum air humidity and incubated for 7-14 days at plant growth conditions. Spores that had been pressed out from pycnidia were isolated using a sterile syringe needle and resuspended into 50 ml YMS medium with 50% glycerol. Cells were suspended by 30 min vortexing on a VXR basic Vibrax (IKA, Staufen, Germany) at 2000 rpm, plated on YMS agar, and incubated for 7 days at 18°. Colonies arising from single cells were PCR screened and their complement on accessory chromosomes characterized as described below.

Phenotypic characterization of infected leaves
Harvested leaves were taped to a sheet of paper and scanned using a flatbed scanner at a resolution of 2400 dpi. The scanned leaves were analyzed using ImageJ (Schneider et al. 2012) and a plug-in described previously (Stewart et al. 2016). The measurement of pycnidia/cm 2 was used for further statistical analyses in R (Ihaka and Gentleman 1996). The Wilcoxon rank-sum test and the Holm's correction for multiple testing were applied to assess statistical differences between the samples.

In vitro phenotype assay
A spore solution containing 10 7 cells/ml and 10-fold dilution series to 1000 cells/ml was prepared. To test for responses to different stress conditions in vitro, YMS plates containing NaCl (0.5 and 1 M), sorbitol (1 and 1.5 M), Congo Red (300 mg/ml and 500 mg/ml), H 2 O 2 (1.5 and 2 mM), methyl methanesulfonate (0.01%); an H 2 O-agar plate; and two plates containing only YMS were prepared. Three microliter of the spore suspension dilutions were pipetted on the plates and incubated at 18°for 7 days. One of the YMS plates was incubated at 28°to test for thermal stress responses.
In vitro growth assay YMS liquid cultures containing 10 5 cells/ml of either Zt09 or the chromosome-loss strains Zt09D14 and Zt09D21 were prepared. Three replicates per strain were used in each experiment. The spores were grown for 4 days at 18°at 200 rpm in 25 ml YMS and the OD 600 was measured at different time points throughout the experiment. The data were fitted using the R package growthcurver (Sprouffske and Wagner 2016) and the r values of each replicate used for statistical comparison. The Wilcoxon rank-sum test was used to compare the r values of the different strains.

Pulsed-field gel electrophoresis and Southern blotting
Fungal strains were grown in YMS medium for 5 days. Cells were harvested by centrifugation for 10 min at 3500 rpm. We used 5 3 10 8 cells for plug preparation. The cells were resuspended in 1 ml H 2 O and mixed with 1 ml of 2.2% low-range ultra-agarose (Bio-Rad, Munich, Germany). The mixture was pipetted into plug casting molds and cooled for 1 hr at 4°. The plugs were transferred to 50 ml screw cap Falcon tubes containing 5 ml of lysis buffer (1% SDS, 0.45 M EDTA, 1.5 mg/ml proteinase K; Roth, Karlsruhe, Germany), and incubated in lysis buffer for 48 hr at 55°, replacing the buffer once after 24 hr. Chromosomal plugs were washed three times for 20 min with 1 3 TE buffer before storage in 5 ml of 0.5 M EDTA at 4°. Pulsed-field gel electrophoresis (PFGE) was performed with a CHEF-DR III pulsed field electrophoresis system (Bio-Rad). To separate the small accessory chromosomes, the following settings were applied: switching time 50-150 s, 5 V/cm, 120°angle, 1% pulsed-field agarose (Bio-Rad) in 0.5 3 TBE (Tris/borate/EDTA) for 48 hr. Separation of midsize chromosomes was conducted with the settings: switching time 250-1000 s, 3 V/cm, 106°angle, 1% pulsed-field agarose in 0.5 3 TBE for 72 hr. Saccharomyces cerevisiae chromosomal DNA (Bio-Rad) was used as size marker for the accessory chromosomes, and Hansenula wingei chromosomal DNA (Bio-Rad) for mid-size chromosomes. Gels were stained in ethidium bromide staining solution (1 mg/ml ethidium bromide in H 2 O) for 30 min. Detection of chromosomal bands was performed with the GelDocTM XR+ system (Bio-Rad).
Southern blotting was performed as described previously (Southern 1975) using digoxigenin-labeled probes generated with the PCR digoxigenin labeling Mix (Roche, Mannheim, Germany), following the manufacturer's instructions.
Sequencing and genome comparison DNA for whole genome sequencing was prepared as described in Allen et al. (2006). Library preparation with a reduced number of PCR cycles (four cycles) and sequencing were performed by the Max Planck Genome Center, Cologne, Germany (http://mpgc.mpipz.mpg.de/home/) using an Illumina HiSeq3000 machine, obtaining 303 coverage with 150 nt paired-end reads. Raw reads were quality filtered using Trimmomatic and mapped to the reference genome of Z. tritici or Z. ardabiliae using bowtie2. SNP calling was conducted with samtools and bcftools and further quality filtered. The filtered SNPs were validated by manual inspection. The data were visualized in Integrative Genomics Viewer (http://software. broadinstitute.org/software/igv/) (Thorvaldsdóttir et al. 2013). Detailed commands for each step of the data processing pipeline can be found in the supplemental material.

Data availability
The authors affirm that all data necessary for confirming the conclusions of the article are present within the article, figures, and tables. Strains are available upon request. Table S1 is a list of primers used for screening accessory chromosomes. Table S2 is a list of identified chromosome-loss strains. Table  S3 shows chromosome variation of the sequenced strains grown at 28°. Table S4 lists location and annotation of SNPs and INDELs found in the sequenced strains. Table S5 displays sequence annotation of chromosomal breakpoints in the sequenced strains derived from the temperature stress experiment.
Supplemental figures are available in supporting figures and legends file. File S1 contains additional methods section including code used to analyze the genomic data. Genomic sequencing raw reads (FASTQ files) are available online at Sequence Read Archive under BioProject ID PRJNA428438. Supplemental material available at Figshare: https://doi.org/ 10.25386/genetics.6854819.

Accessory chromosomes are lost at a very high rate in vitro
To assess the stability of accessory chromosomes in Z. tritici we conducted an in vitro long-term growth experiment. We used the Z. tritici isolate IPO323Dchr18, derived from the reference strain IPO323 for which there is a completely assembled genome ). IPO323Dchr18 lost chromosome 18 during previous in vitro propagation, and is here referred to as Zt09 (Kellner et al. 2014).
We propagated fungal cells in liquid culture at 18°, including eight transfers to fresh medium. After 4 weeks (representing 80 cell divisions per cell), 576 single strains originating from three replicate cultures were tested for the presence of the seven accessory chromosomes from the Zt09 progenitor by a PCR assay. The presence of marker regions, located close to the centromere, and the right and left telomere repeats of each accessory chromosome, were tested (Table S1). When the screening by PCR suggested absence of an accessory chromosome, we further validated the result by electrophoretic separation of the fungal chromosomes by PFGE (Figure 1, A and B). Thirty eight of the 576 (7%) tested strains lacked one accessory chromosome. We did not find strains that lacked .1 chromosome, but observe a clear trend where accessory chromosomes 14, 15, and 16 were lost more frequently than smaller accessory chromosomes (Table 1). For example, chromosome 14 was absent in 18 strains, chromosome 15 in eight, and 16 in nine strains. The small chromosomes 20 and 21 were lost only in two and one strain, respectively, whereas chromosomes 17 and 19 were never absent following the in vitro propagation (Table S2 and Table 1).
We considered the possibility that natural selection acted on a large population of cells in the liquid Z. tritici cultures and thereby contributed to the observed nonrandom chromosome losses. To test this, we designed a second evolution experiment to propagate individual cell lineages in the absence of selection. For this experiment, 40 strains originating from one common progenitor (Zt09) were cultivated on plates. Every week, randomly selected colonies derived from single cells were transferred to a fresh plate. These repeated strong bottlenecks allowed us to propagate single lineages of Z. tritici without an effect of natural selection (Lynch et al. 2008). After 4 weeks of growth (including four transfers to new plates), the 40 evolved strains were tested for the presence of accessory chromosomes as described above. Five out of 40 (13%) strains, twice as many as in the previous Figure 1 Screening by PCR and pulsed-field gel electrophoresis (PFGE) to identify accessory chromosome-loss strains in Z. tritici. (A) Screening by PCR of the progenitor Zt09 strain and evolved strains to assess chromosome losses using primer pairs located close to the centromeric region of the respective accessory chromosome. Multiplex PCR was used to simultaneously screen for the presence or absence of all accessory chromosomes in one strain. Absence of a PCR product indicates the absence of the chromosome; chromosome 18 is not present in Zt09. (B) PFGE was conducted to validate the absence of accessory chromosomes in chromosome-loss candidates identified by the initial PCR. Here, the separation of small accessory chromosomes of the progenitor strain Zt09 and four chromosome-loss strains is shown. The corresponding chromosome to each band is indicated on the right. The absence of chromosomal bands confirms the loss of the respective chromosome; chromosome 18 is similar in length to 16 and 17 but is absent in Zt09. PFGE of the accessory (C) and midsize (D) chromosomes of strains originating from the in vitro temperature stress experiment confirms multiple accessory chromosome losses, size alterations, and chromosome fusions (see also Figure  3, Figure S1, and Table S3). The bands for chromosomes 3, 12, 13 and 21 are absent in strain 28-2; however, genome sequencing confirms that the chromosomes are still present in the genome (Figure 2B). The absence of bands on the PFGE can be explained by chromosome size changes whereby chromosomes 3 and 13, and 12 and 21 have experienced chromosome fusions (Figure 3 and Figure  S1). Note a new chromosome band with a size of 1.8 Mb resulting from the fusion of chromosomes 12 and 21. In the strain 28-3, the accessory chromosomes 17 and 19 are shorter than the reference chromosomes due to chromosome breakage (Figure 2B and Table S3). Strain 28-4 has an additional band of 1.2 Mb, representing the fusion of a duplicated chromosome 17 ( Figure S1). All images of stained gels are color-inverted to make differences more obvious. experiment, lacked an accessory chromosome suggesting that selection in a larger population of cells indeed had removed some of the spontaneously occurring chromosome losses. Interestingly, in this experiment we observe the loss of only two different chromosomes: chromosome 16 was lost in four strains, and chromosome 21 in one strain. Our two in vitro experiments suggest that chromosome loss is not entirely random, and that specific chromosomes are lost at a higher rate (Table 1), depending on the in vitro growth conditions. The frequency of chromosome losses during asexual growth in vitro is extremely high in Z. tritici and exceeds previously reported spontaneous chromosome losses in Fusarium oxysporum f. sp. lycopersici by far (Vlaardingerbroek et al. 2016).

Accessory chromosomes are unstable in planta
The growth of Z. tritici in rich medium in vitro is highly distinct from growth of the fungus in its natural environment. Z. tritici is a hemi-biotrophic pathogen infecting the mesophyll of wheat leaves Ponomarenko et al. 2011), where it forms asexual fruiting bodies, pycnidia, containing presumably clonal conidiospores (pycnidiospores). Pycnidia develop in the substomatal cavities, where the environment and nutrient availability differ from our tested in vitro conditions (Rudd et al. 2015;Haueisen et al. 2017). To test whether chromosome loss occurs only in vitro or also during the natural lifecycle of Z. tritici, we assessed the loss of accessory chromosomes in planta. We collected pycnidia from wheat leaves infected with the Z. tritici strain IPO323 and tested single pycnidiospores for the presence of accessory chromosomes. A total number of 968 pycnidiospores originating from 42 separate pycnidia were screened by PCR. In total, 17 strains missing an accessory chromosome were identified (1.7%). Chromosomes 15 and 18 were the most frequently lost chromosomes, while chromosomes 17, 20, and 21 were not lost in any of the strains (Table S2 and Table 1). As for the in vitro experiments, no strains lacking .1 chromosome were found. Interestingly, chromosome 18, which was found to be lost at high rates in our plant experiments, was shown to be frequently absent in field isolates of Z. tritici (Croll et al. 2013;McDonald et al. 2016). In summary, our findings show that chromosome loss in Z. tritici not only results from the nondisjunction of homologous chromosomes during meiosis II, as proposed previously (Wittenberg et al. 2009), but also from frequent chromosome losses during asexual growth and mitotic spore formation in planta. The observed frequency of chromosome losses in planta is not as high as during in vitro growth; however, the number of mitotic divisions to develop pycnidia is presumably lower than during 4 weeks of in vitro growth. Therefore, we propose that the reduced number of chromosome losses in planta rather reflects the reduced number of cell divisions than possible fitness effects of chromosome losses.

Chromosome instability is greatly increased during exposure to heat stress
In their natural environments, pathogens are exposed to various kinds of biotic and abiotic stresses. For example, the local environment on the leaf surface can fluctuate severely in temperature and humidity conditions (Zhan and McDonald 2011). We assessed the effect of an increase in temperature on chromosome stability of Z. tritici. To this end, we cultivated Zt09 in vitro at elevated temperatures, namely an increase of 10°from 18 to 28°during 4 weeks of incubation. Subsequent screening by PCR for accessory chromosomes revealed severe chromosome losses in the tested strains. Out of 188 evolved and tested strains, 148 (80%) were missing at least one accessory chromosome. Chromosome 14 was the most frequently lost chromosome and 34 evolved strains were lacking .1 chromosome. The maximum number of missing chromosomes was six in one strain (Table S2 and Table 1). PFGE revealed that karyotype alterations were not restricted to simple chromosome loss as observed in the in vitro experiments at 18°and the in planta experiments. In contrast, we observed frequent size variation of both core and accessory chromosomes as a consequence of chromosome breakage, fusions, and duplications (see below) ( Figure 1, C and D, Figure S1, and Table S3).
Genome sequencing of chromosome-loss strains reveals chromosome breakage and fusion, but a low number of SNPs In total, we sequenced 19 genomes of strains that had lost different chromosomes or displayed rearrangements of Listed are the number of strains that lost an accessory chromosome and which accessory chromosomes were lost. Sizes of the accessory chromosomes are listed next to the chromosome number. Strains with .1 chromosome lost are listed in brackets.
chromosomes as revealed by PFGE. The strains represented derived strains from the different in vitro and the in planta experiment and the respective progenitor strains; all sequencing results have been deposited in the Sequence Read Archive under BioProject ID PRJNA428438. Overall, we confirmed the absence of complete accessory chromosomes (Figure 2A), and we found very few additional changes. After filtering to exclude reads of poor quality and coverage (see Materials and Methods and Supplemental Material) we identified a total number of nine SNPs in eight sequenced strains originating from the in vitro experiment in liquid culture at 18°and the in planta experiment. Of these, three SNPs were located in coding regions (Table S4). Besides the complete loss of accessory chromosomes and the identified SNPs, we found no other mutations when comparing the progenitor and evolved strains. However, we identified SNPs distinguishing Zt09 from our IPO323 reference strain, and the published genome sequence; most of these SNPs were in noncoding sequences (Table S4). As we show chromosome losses in different strains, we conclude that the few point mutations do not have a measurable effect on chromosome loss. Whole genome sequencing of five strains evolved at 28°v erified the losses of several accessory chromosomes. Furthermore, comparison of the resequenced genomes of the progenitor and the evolved strains revealed substantial size variation resulting from chromosome breakage of core and accessory chromosomes, and duplications of accessory chromosomes ( Figure 2B). Chromosome breaks located close to the ends of chromosomes resulted in shortened chromosomes caused by subtelomeric deletions of 0.2-60 kb (Table S3). Detailed analyses of the distribution of discordant paired-end reads mapping to different chromosomes revealed a high number of reads at the chromosome breakpoints with the respective read mates mapping to telomeric repeats ( Figure 3). This suggested that most chromosome breaks resulted in shorter chromosomes to which telomeres were added de novo. Besides de novo telomere formation, we found evidence for fusion of core chromosomes in strain Zt09 28-2, where discordant read mapping indicated fusion of the right arms of chromosomes 3 and 13 ( Figure 3 and Table S3). PFGE and Southern blots further showed fusion of chromosomes 12 and 21 in Zt09 28-2 ( Figure S1). In a previous study, we reported evidence for the likely fusion of an ancestral accessory chromosome and core chromosome 7 in the reference isolate IPO323 and Zt09 (Kellner et al. 2014;Schotanus et al. 2015). Here, we deduced a fusion of the duplicated chromosome 17 in strain Zt09 28-4 (Figure 1, C and D and Figure S1). Similar fusions of the same accessory chromosome 17 following meiosis had been suggested previously, and "breakage-fusion bridge" cycles (McClintock 1938(McClintock , 1941 were invoked as a mechanism to form a new accessory chromosome (Croll et al. 2013). Based on the position of centromeres in Z. tritici (Schotanus et al. 2015), all three chromosome fusion events that we detected, resulted in the formation of dicentric chromosomes. On two chromosomes involved in the fusion events, we identified the "internal" breakpoints (Table S3) likely reflecting the consequences of chromosome breakage of a dicentric chromosome during mitosis. This indicates that breakage-fusion bridge cycles also occur during mitosis in Z. tritici and can be a mechanism involved in the structural variation detected in our experiments. In total, in the five analyzed genomes derived from strains grown at 28°, we identified the following spontaneous changes: 15 chromosome breakages (12 with de novo telomere formation and three without evidence for new telomeres), three chromosome fusions, three chromosome duplications, and 16 chromosome losses (Table S3). The breakpoints of the structural rearrangements often colocalize with annotated transposons, predominantly class I elements (Table S5). Transposable elements are often associated with structural rearrangements (Faino et al. 2016;Plissonneau et al. 2016) and their activation under stress conditions might contribute to the observed genome instability (Chadha and Sharma 2014). A key observation from our genome data analyses is the frequent involvement of de novo telomere formation and chromosomal fusion as a mechanism to heal broken chromosome ends. While this has been reported for cancer cells (Murnane 2012), it is generally considered a rare event in normal, nontransformed cells. Here we show that de novo telomere formation and chromosome fusion readily occur in a filamentous fungus during growth at temperature stress.

Accessory chromosome instability also occurs in Zymoseptoria sister species
To assess whether the high frequency of mitotic chromosome loss is specific to Z. tritici, we conducted a short-term in vitro evolution experiment on the fungus Z. ardabiliae, a closely related sister species of Z. tritici that infects wild grasses. We used the previously characterized isolate STIR04 1.1.1 (Stukenbrock et al. 2011), here called Za17, in our experiments.
We first generated a high-quality reference genome based on long-read single molecule real time (SMRT) sequencing (Supplemental Material). Analyses of a population genomic dataset revealed at least four accessory chromosomes varying in their frequency among 17 Z. ardabiliae isolates (Stukenbrock et al. 2011;Stukenbrock and Dutheil 2018). Based on PFGE and homology of Z. ardabiliae unitigs to fragments of Z. tritici IPO323 accessory chromosomes (Soderlund et al. 2011), we designed primers to amplify six loci located on putative accessory chromosomes of Za17 (Table S1).
Next, we conducted an in vitro evolution experiment at 18°f or 4 weeks in liquid culture with Za17 as the progenitor strain. Screening of 288 single clones by PCR identified five strains lacking one of the accessory chromosomes (1.7%). We further verified the absence of accessory chromosomes by PFGE for all five strains and whole genome sequencing for three of the five Z. ardabiliae strains ( Figure 2C). By combining PFGE and genome sequencing of chromosome-loss strains we could confirm unitigs 24, 34, 49 and 51 to be accessory sequences pertaining to three accessory chromosomes ( Figure 2C).
These findings resemble the rapid loss of chromosomes in Z. tritici and suggest a common phenomenon, most likely related to mitotic cell division in the two Zymoseptoria species.

Accessory chromosome losses affect in planta phenotypes
We next addressed the effect of spontaneous chromosome loss by comparing the fitness of the progenitor and evolved Z. tritici strains under different in vitro growth conditions, and during infection of a susceptible wheat variety by comparing growth and pycnidia formation. In vitro, the growth rate of strains lacking an additional accessory chromosome (Zt09D14 and Zt09D21) was comparable to the growth rate of the progenitor, Zt09, but with a slight tendency to slower growth of the derived chromosome-loss strains ( Figure S2).
For a more detailed phenotypic characterization, we used the five chromosome-loss strains isolated from pycnidia from the in planta experiment (Table S2). Several stress conditions including osmotic, oxidative, temperature, and cell wall stresses were tested in vitro. We compared the fungal phenotypes of the chromosome-loss strains and the progenitor strain IPO323 and observed no difference in growth rate or colony morphology between any chromosome-loss strain and IPO323 ( Figure 4A). In planta, however, we observed phenotypic differences between the evolved chromosome-loss strains and their progenitor. We measured fitness by counting the number of asexual fruiting bodies formed in stomata on infected wheat leaves and found a slightly higher fitness of the chromosome-loss strains compared to the progenitor strain ( Figure 4B). The difference was most pronounced for Figure 2 Whole genome sequencing confirms loss of accessory chromosomes in Z. tritici and Z. ardabiliae. The genomes of chromosome-loss strains derived from the Z. tritici in vitro evolution experiment, the in planta pycnidiospore isolations, the Z. ardabiliae in vitro evolution experiment, and the respective progenitor strains were sequenced by paired-end Illumina sequencing and mapped to the respective reference isolates, (A and B) IPO323 (Z. tritici) or (C) Za17 (Z. ardabiliae). The losses of entire chromosomes (white boxes) were verified for each strain, but few mutations in form of SNPs, indels, or copy number variation were detected (Table S4). The Z. tritici reference genome consists of whole chromosomes, while the Z. ardabiliae reference genome is composed of unitigs obtained from the SMRT sequencing assembly (ordered by name of the unitig). The read depth is here normalized to 1 X per strain. (B) Genome sequencing of strains derived from the temperature stress experiment revealed, besides verification of accessory chromosome losses, chromosome breaks at the ends of core and accessory chromosomes, and accessory chromosome duplications. Zooming in on the coverage of accessory chromosomes highlights the duplications of chromosomes 17 (Zt09-28-2, Zt09-28-4) and 20 (Zt09-28-5), and the chromosome breakage of chromosomes 17 and 19 (Zt09-28-3). Darker blue shading indicates higher coverage. Most repetitive sequences have a higher coverage than single copy regions, resulting in different shades of blue in the coverage graph. The intense dark blue lines indicate high coverage regions on ribosomal DNA (rDNA) clusters or repetitive DNA due to underestimation of repeats in the reference assembly (for example, the thin black line in chromosome 7 indicates the location of several rDNA cluster repeats, the only such repeats in the annotated genome sequence, but the expected true number of rDNA repeats is 50). Accessory chromosomes or unitigs are highlighted in grey.
the strains lacking chromosome 14 and chromosome 19 where the density of pycnidia was found to be significantly higher (P values of 0.0003 and 0.015, Wilcoxon rank-sum test and Holm's correction for multiple testing) compared to the IPO323 progenitor ( Figure 4B).

Consequences for evolution of accessory chromosomes
Accessory chromosomes were proposed to serve as an "evolutionary cradle" for creating novel virulence genes without risking the disruption of essential genes on the core chromosomes (Croll and McDonald 2012). Indeed, signatures for accelerated evolution were shown by overall higher dN/dS ratios of genes on the accessory chromosomes compared to the core chromosomes (Stukenbrock et al. 2010). However, natural selection acts at the level of individuals and cannot maintain "a playground" for future beneficial effects. Rather, we hypothesize that the dynamic of accessory chromosome loss reflects a dynamic in the selective environment of Z. tritici, where under certain conditions these small chromosomes confer a fitness advantage, while under different Figure 3 Discordant read analysis reveals chromosome fusion and chromosome breakage. (A) Schematic illustration of discordant read mapping in the case of chromosome breakage and fusion. If chromosome breakage is followed by fusion to a different chromosome, the discordantly mapping reads at the breakpoint will have their respective read mate on the chromosome that is fused to the breakpoint. If breakage is followed by de novo telomere formation, the discordant reads have their respective read mates on telomeric repeats on random chromosomes. (B and C) Read mapping to the genome of the reference isolate IPO323. (B) In the experimentally evolved strains of Z. tritici at 28°, the fusion of chromosomes 3 and 13 in the strain Zt09 28-2 is indicated by the increased occurrence of discordant reads at the chromosome breakpoint of the right arm of chromosome 3 and the right arm of chromosome 13. The color of the reads represents the chromosome their respective mate is mapping to. (C) New telomere formation at the breakpoint, here, as example, shown for chromosome 9 in the strain Zt09 28-2. Telomere formation is indicated by a high number of discordant reads, where the read mates are mapping to the telomeric repeats of different random chromosomes (dashed line). A total of 12 breakage events that were followed by de novo telomere formation were detected in the five analyzed strains derived from the experiment at elevated temperature.
conditions they cause decrease in fitness. A fitness advantage of chromosome-loss strains in planta may indicate the presence of avirulence factors on at least some accessory chromosomes. Avirulence factors can be recognized by the host to induce a resistance response resulting in reduction in virulence or even abortion of infection (Petit-Houdenot and Fudal 2017;Zhong et al. 2017). Our previous findings support this notion, as we found small but significant negative effects of accessory chromosomes during host infection (Habig et al. 2017). In contrast, Fouché et al. (2018) detected small, quantitative increases of virulence traits if specific accessory chromosomes were present. However, so far no avirulence factors and very few potential effectors were identified on accessory chromosomes in Z. tritici and the vast majority of annotated genes on accessory chromosomes does not have a predicted function Figure 4 In vitro and in planta phenotype assays of the Z. tritici reference isolate IPO323 and the in planta chromosome-loss strains. (A) Several in vitro stress conditions were tested to assess the effect of accessory chromosome losses on fitness. Five in planta chromosome-loss strains were tested, IPO323 was used as the reference strain. We observed no noteworthy differences in growth rate or colony morphology between reference and chromosome-loss strains. (B) To investigate the effect of accessory chromosome losses on Z. tritici infection of wheat, three independent experiments with the reference strain IPO323 and five chromosomeloss strains were conducted. The fitness of each strain was measured by counting the number of pycnidia per square centimeter on the leaf surface. Statistical analysis (Wilcoxon rank-sum test and Holm's correction, *P # 0.05, ***P # 0.001) identified two chromosome-loss strains (chromosomes 14 and 19 P-values 0.0003 and 0.015) that show a significantly higher number of pycnidia compared to the reference strain IPO323 under the conditions used in this experiment. Grandaubert et al. 2015). In other pathogenic fungi, however, the opposite effect has been demonstrated in several cases. A strong positive effect on virulence of accessory chromosomes has been observed in F. oxysporum f. sp. and Nectria haematococca (Coleman et al. 2009;Ma et al. 2010;van Dam et al. 2018). In these pathogens chromosome losses resulted in complete loss of pathogenicity. While these cases are currently considered the norm, our observations show that responses to accessory chromosome loss can be more varied. Even though this fungus is an important pathogen in all wheat-growing countries, the biology of Z. tritici is still not completely understood. Mating and overwintering in the soil or litter layer remain largely unknown lifecycle stages, and we postulate that they impose different selection pressures on the fungus. During these stages, accessory chromosomes may confer significant fitness advantages that have not been demonstrated yet. Nevertheless, any effects on fitness found to date are rather small. While beneficial effects, even if they are comparably small, might explain the maintenance of accessory chromosomes in Z. tritici populations, negative effects as observed in this study and reported by Habig et al. (2017) raise the question why these chromosomes have not been lost over time, especially as they can easily be lost during meiosis and mitosis. Accessory chromosome can be lost during meiosis (Wittenberg et al. 2009), but interestingly accessory chromosomes have also been shown to follow non-Mendelian segregation and to be transmitted at a significantly higher rate than expected, if one of the two mating partners lacks an accessory chromosome (Fouché et al. 2018;Habig et al., 2018). While there is no comprehensive mechanistic explanation yet, chromosome conservation by re-replication during meiosis may counteract the chromosome loss during vegetative growth and therefore maintain accessory chromosomes in pathogen populations. In Alternaria alternata, the spontaneous loss of a conditionally dispensable chromosome during subculturing has been described for one isolate (Johnson et al. 2001). Furthermore, loss of a conditionally dispensable chromosome after meiosis has been observed in N. haematococca (aka Fusarium solani MPVI) (Miao et al. 1991). A recent study of the plant pathogenic fungus F. oxysporum f. sp. lycopersici describes the spontaneous loss of dispensable chromosomes in vitro in one cell out of 35,000 (Vlaardingerbroek et al. 2016). This rate is much lower than the chromosome-loss rate that we observe in Z. tritici and Z. ardabiliae, where, depending on the growth conditions, between 2 and 80% of the tested strains, lack an accessory chromosome. In combination, these results indicate that chromosome instability is a common phenomenon in fungi, but the stability of chromosomes can vary substantially between different species. Studies in Saccharomyces and Candida species have shown that chromosome loss resulting in aneuploidy is a relatively frequent phenomenon in diploid cells and is advantageous under certain conditions (Selmecki et al. 2010;Kumaran et al. 2013;Bennett et al. 2014;Zhu et al. 2014). The loss of chromosomes in haploid organisms, such as Z. tritici, is expected to have more severe consequences, as genetic information is lost from the cell. We considered, however, that the effect of rapid chromosome loss or instability may be advantageous if these traits are selected for in pathogens like Zymoseptoria.
The temperature stress experiment in vitro showed that genome instability in Z. tritici, at this temperature, is not restricted to the accessory chromosomes, but also involves the core genome. New chromosome formation induced by telomere to telomere fusion followed by breakage of the dicentric chromosome was reported in Cryptococcus neoformans (Fraser et al. 2005) but likely occurred during meiosis rather than mitosis. In the asexual Candida glabrata, clinical isolates frequently inhabit chromosomal rearrangements indicating that structural variation acts as a virulence mechanism (Poláková et al. 2009) and might be a response to stresses such as antifungal treatment or elevated temperature as observed here. The instable regions close to the chromosome ends on the core chromosomes, however, show structural characteristics similar to the accessory chromosomes, such as higher repeat content, lower gene density and specific histone modification patterns (Dhillon et al. 2014;Grandaubert et al. 2015;Schotanus et al. 2015).
The accessory chromosomes of Z. tritici do not differ in terms of centromere or telomere organization from the core chromosomes, but the chromatin is largely transcriptionally inactive and potentially more condensed due to an enrichment with the histone H3 that is trimethylated at lysine 27 (H3K27me3), which covers almost all of the accessory chromosomes (Schotanus et al. 2015). Enrichment of the facultative heterochromatin mark H3K27me3 seems to be a common feature of accessory chromosomes in several fungal pathogens but the functional relevance for stability and transcriptional regulation of these chromosomes has not been addressed so far (Galazka and Freitag 2014;Seidl et al. 2016). H3K27me3 is also enriched in subtelomeric regions of core chromosomes, that we show here as prone to instability. High-resolution microscopic analyses indicate an unusual localization of centromeres in the Z. tritici nuclei, suggesting a distinct spatial organization of different chromosomes (Schotanus et al. 2015). Heterochromatin has been found close to the nuclear periphery, and H3K27me3 has been shown to be involved in facilitating lamina-proximal positioning (Harr et al. 2015). We hypothesize that the heterochromatic structure reflects a distinct physical organization of core and accessory chromosomes and likely the subtelomeric regions in the nucleus of Z. tritici, and that this is correlated to the instability of heterochromatic regions and chromosomes. An Hi-C study on Neurospora crassa showed that absence of H3K27me3 resulted in movement of subtelomeric regions and centromeres from the nuclear periphery into the nuclear matrix . Further experiments focusing on localization of specific Zymoseptoria chromosomes in the nucleus by cytology or Hi-C  should be conducted to address this hypothesis.
In conclusion, using experimental evolution, electrophoretic karyotyping and genome sequencing we showed that the accessory chromosomes of Z. tritici are highly unstable during asexual growth across numerous mitotic cell divisions ("mitotic growth") in vitro as well as in planta. Surprisingly, increasing the temperature from 18 to 28°d ramatically increased overall genome instability in Z. tritici. Besides chromosome losses we observed structural variation in form of chromosome breakage, duplication, and fusion involving both core and accessory chromosomes, and all events were increased near telomeric sequences. In Z. tritici, all of these chromosome abnormalities have thus far been associated with meiosis (Wittenberg et al. 2009;Croll et al. 2013), however, our study highlights an important role of mitotic growth in generating genetic diversity. The variability in chromosome number and structure obtained over a relatively short period of time (during one infection event, or 4 weeks of vegetative growth) correlates with the genomic diversity observed in field populations of this important plant pathogen (McDonald and Martinez 1991;Linde et al. 2002;Zhan et al. 2003;Mehrabi et al. 2007). Our findings demonstrate that even sexual fungal pathogens can accelerate the generation of new genetic diversity by mitosis-associated structural variation. The mitotic events responsible for these rearrangements are still to be characterized.
Lastly, chromosome instability has been observed in various organisms and is a frequent phenomenon in cancer cells (McGranahan et al. 2012). Epigenetic factors and heterochromatin in particular have been shown to affect chromosome stability in human cancer cells (Slee et al. 2012). The accessory chromosome dynamics correlating with distinct chromatin patterns suggest that Z. tritici may provide a new model for mechanistic studies of chromosome instability in cancer cells.