The chromodomain proteins, Cbx1 and Cbx2 have distinct roles in the regulation of heterochromatin and virulence in the fungal wheat pathogen, Zymoseptoria tritici

Heterochromatin is characterized by specific histone post-translational modifications such as the di- and tri-methylation of histone H3 on lysine 9 (H3K9me2/3), which direct the recruitment of ‘reader’ proteins to chromatin. In the fungal phytopathogen, Zymoseptoria tritici, deletion of the H3K9 methyltransferase gene kmt1, results in a global increase in the expression of transposable elements (TEs), genome instability and loss of virulence. Here we have identified two Z. tritici chromodomain proteins, Cbx1 and Cbx2, that recognise H3K9me modifications. Cbx1 is a Heterochromatin Protein 1 homolog that binds H3K9me2/3 in vitro and associates with heterochromatic loci in vivo. Transcriptomic analysis also indicates that Cbx1 and Kmt1 regulate overlapping sets of protein-encoding genes. However, unlike Δkmt1 mutants, Δcbx1 strains do not exhibit a global increase in TE expression and have only a partial reduction in virulence, suggesting the existence of additional H3K9me reader proteins. Accordingly, we have identified a fungal-specific chromodomain protein, Cbx2, that binds H3K9me3 in vitro. Strikingly, the growth defects of Δcbx1 Δcbx2 double mutants closely resemble those of Δkmt1 consistent with Cbx1 and Cbx2 playing redundant roles in gene silencing. Overall, the data suggest that key functions of H3K9me modifications are mediated by a combination of Cbx1 and Cbx2.


INTRODUCTION
also regulated by heterochromatin (Chujo & Scott, 2014). These findings have led to signal ( Fig 2C). Next, chromatin immunoprecipitation (ChIP) assays were used to investigate the ability of Cbx1 to associate with H3K9me-enriched regions of the 138 genome ( Fig 2D). A strong enrichment of Cbx1-GFP was observed at a TE 139 (DTH_element 299 5_ZTIPO323) that is known to be associated with H3K9me 140 (Schotanus et al., 2015). Furthermore, a similar enrichment of Cbx1-GFP was also 141 found at a H3K9me-marked subtelomeric region (Chromosome 1: 161011-175642). 142 Importantly, Cbx1-GFP enrichment was not detected at the euchromatic (H3K4me3-143 associated) genes, actin (Mycgr3G105948) and GAPDH (Mycgr3G99044). Notably, 144 GAPDH, is located adjacent to an H3K9me-marked DNA transposon, suggesting that 145 the resolution of the assay was sufficient to distinguish between neighbouring 146 H3K9me-marked and non-marked genes. Taken together the data indicate that Cbx1 147 is an HP1 ortholog and is likely to function in the recognition of H3K9me2/3 in Z. tritici.  151 The role of HP1 in the fitness of Z. tritici was investigated by generating cbx1 152 deletion strains. For comparison, we also constructed strains lacking the H3K9 153 methyltransferase gene, kmt1 and confirmed the loss of H3K9me3 marks in these 154 mutants ( Fig S1A). Initially, the in vitro growth and stress-sensitivity profiles of the 155 Δcbx1 and Δkmt1 mutants relative to the IPO323 reference strain were determined. 156 As previously reported, deletion of kmt1 resulted in a slow growth phenotype (Möller isolates and the biological replicates for all strains, indicative of low variation ( Fig S2).
Next, hierarchical clustering was employed to provide an overview of the global 188 similarities between the transcriptomes of the sequenced strains. This revealed that 189 the transcript profiles of Δcbx1 and Δkmt1 mutants exhibit an overall similarity. Indeed,190 the Δcbx1 mutant profiles were found to be more similar to the Δkmt1 strain than to 191 the reference strain ( Fig 5A). that were differentially expressed in Δcbx1 but not Δkmt1 and vice versa (Table S3 212 and S4). The existence of these non-overlapping sets of DE genes may, at least in 213 part, explain the differences in the phenotypes associated with Δcbx1 and Δkmt1 214 mutants. 215 The similarity of the RNA-seq data from the Δkmt1 strain generated in this study 216 and the Zt09-Δkmt1 strain (Möller et al., 2019) was also analysed. For this comparison 217 the more stringent cut-offs (4-fold change in expression, adjusted p value < 0.001) To determine whether this was also the case for 234 mutants lacking cbx1, normalised read counts were mapped from genes on accessory 235 chromosomes. As expected a significant increase in read counts from accessory chromosomes was observed in the Δkmt1 mutant compared to IPO323 (p < 0.014; Δcbx1 background ( Fig 6A). Therefore, loss of Cbx1 is not sufficient for a global 239 increase in expression from genes on accessory chromosomes.

240
The accessory chromosomes of Z. tritici are highly enriched in TEs and 241 previously it has been shown that expression from these elements is suppressed by 8.25e-09 ANOVA), however in contrast, no significant global increase the Δcbx1 246 mutant was detected ( Fig 6B). Indeed, hierarchical clustering revealed that with 247 respect to the profile of TE expression, the Δcbx1 mutant is more similar to the 248 reference IPO323 strain than to the Δkmt1 mutant ( Fig 6C). Therefore, although global 249 silencing of TEs in Z. tritici requires Kmt1, and by implication H3K9me2/3, it is not 250 dependent upon recognition of these histone modifications by the HP1 protein, Cbx1. TEs. Genetic instability around such loci in filamentous fungi has been proposed to 256 drive rapid evolution and aid niche adaptation (Dong et al., 2015, Faino et al., 2016, 257 Laurent et al., 2018. Therefore, we analysed the expression of all genes within 2 kb 258 of a TE. In total 1505 genes were identified as 'TE-associated' of which 184 were 259 differentially expressed in Δcbx1 and 205 in Δkmt1, a highly significant enrichment in 260 both cases (p < 4.35e-07 and p < 9.66e-08 respectively). Furthermore, 114 TE-associated genes were found to be commonly differentially expressed in Δcbx1 and 262 Δkmt1 (p < 5.55e-65) ( Fig 6D). The expression profile of these genes was also found 263 to be highly similar between the Δcbx1 and Δkmt1 mutants and indeed all but 4 genes 264 exhibited similar expression patterns (Fig 6E). GO-term analysis of the differentially 265 expressed TE-associated genes revealed that 8 of the 114 were annotated as having 266 functions relating to secondary metabolism. However, the majority of these were 267 'orphan' genes with no assigned GO terms. This is not unexpected given that TE-268 associated and heterochromatic loci are known to be enriched with 'orphan' genes in 269 a variety of plant pathogenic fungi (Dong et al., 2015). The removal of HP1 tends to result in the upregulation of genes that are associated 274 with H3K9me2/3-marked chromatin (Chujo & Scott, 2014, Reyes-Dominguez et al., 275 2010, Soyer et al., 2014. To investigate whether this is the case in Z. tritici, genes 276 that were partially (>1 bp) or completely associated with H3K9me were identified 277 through analysis of published ChIP-seq data (Schotanus et al., 2015). The 278 relationship between H3K9me-associated genes and genes that are differentially 279 expressed in Δcbx1 was then determined. A total of 247 genes were found to be 280 partially associated with H3K9me (>1 bp) while only 112 were completely associated 281 with this modification. Of the partially H3K9me-associated genes only 31 were 282 differentially expressed in Δcbx1, an overlap which was just statistically significant (p 283 < 0.023) (Fig 7A). The overlap between the completely H3K9me-associated genes 284 and Δcbx1 DE genes was also modest (17 genes, p < 0.016) (Fig 7B). At first glance 285 this weak correlation is surprising, however it has previously been observed that deletion of kmt1 is not sufficient for the upregulation of the majority of H3K9me-287 associated genes in Z. tritici (Möller et al., 2019) and analysis of the Δkmt1 RNA-seq 288 data generated in this was consistent with these findings. Only a very modest overlap 289 was observed between Δkmt1 DE genes and genes that are fully associated with 290 H3K9me and no significant overlap was observed with partially associated genes (  Cbx2, a fungal-specific CD protein that binds to H3K9me3. 296 Comparison of the phenotypes of Δkmt1 and Δcbx1 mutants suggested that some of 297 the downstream effects of H3K9me2/3 histone modifications are likely to be mediated 298 independently of the HP1 homolog Cbx1. One explanation for this would be that Z. 299 tritci has additional H3K9me2/3 reader proteins. Therefore, we used BLAST analyses 300 to search for further proteins with the potential to bind H3K9me2/3 PTMs and identified 301 five hypothetical proteins with CD domains (as predicted by ExPASy Prosite and or 302 Pfam). None of these proteins contained a recognisable CSD, consistent with Cbx1 303 being the sole HP1 isoform in Z. tritici. Four of the hypothetical CD proteins were 304 eliminated from further analysis for one or more of the following reasons, (i) they 305 exhibited similarity to retroviral/retrotransposon integrases, (ii) the CD domain lacked 306 critical key aromatic methyl-lysine caging residues or (iii) they were encoded on an 307 accessory chromosome. The remaining hypothetical protein (Mycgr3G108849, 308 hereafter called Cbx2) was predicted to be 703 amino acids in length and have two 309 CD domains in the C-terminal region ( Fig 8A and Fig S4). BLAST analyses revealed 310 that organisms that encode proteins with homology to Cbx2 extending beyond the CD domains are limited to species in just a few fungal families (principally the the broadly conserved HP1 family member Cbx1, Cbx2 is a fungal-specific CD protein.

314
Sequence analysis revealed that both Cbx2 CDs possess the conserved 315 'aromatic cage' residues that facilitate methyl-lysine binding ( Fig S4) and furthermore, 316 chromodomain 1 (CD1) was predicted to be acidic, a characteristic of HP1-type  The histone peptide binding assays suggested that Cbx2 has the potential to 325 function as an effector of H3K9me3 PTMs and so Δcbx2 strains were generated.

326
Comparison of the Δcbx2 mutant with the IPO323 reference strain indicated that loss 327 of Cbx2 does not result in any detectable reduction in fitness or stress resistance (Fig   328   S6). Furthermore, wheat infection assays revealed that, unlike the Δkmt1 and Δcbx1 329 strains, Δcbx2 strains exhibited no obvious reduction in virulence (Fig 9A and B).

330
Leaves treated with Δcbx2 mutants developed disease symptoms at a very similar rate 331 to those treated with the reference IPO323 strain and there was no major difference 332 in the numbers of pycnidia at 21 dpi ( Fig 9C). As such, loss of Cbx2 alone does not 333 obviously impact the growth of Z. tritici either in vitro or in planta. This is perhaps not 334 surprising, as when we analysed the expression of cbx2 using our RNA seq data, we found that this gene was expressed at similar levels to the H3K9 methyltransferase kmt1, but at only ~2.9% of the level of cbx1. 337 We hypothesized that Cbx2 may co-operate with Cbx1 but that effects of cbx2 338 deletion may be masked when cbx1 is present. As a test of this, a double deletion H3K9me readers in these organisms.

383
The transcriptomic analysis revealed a highly significant overlap in the 384 differentially expressed genes in the Δkmt1 and Δcbx1 backgrounds. Nonetheless, a set of genes was identified whose expression was dependent upon Cbx1, but independent of Kmt1 (Table S3). This suggests that Cbx1 has functions that are 387 independent of H3K9me2/3. Consistent with this, H3K9me-independent roles for HP1  (Cam et al., 2008, 463 Choudhary et al., 2020, Mamillapalli et al., 2013, Sun et al., 2020.

464
Analysis of the genomes of fungal phytopathogens has revealed that effector 465 genes tend to be associated with rapidly evolving regions of the genome that are   beads. 100 μL of 10% Chelex® 100 (w/v) was also added to 10 μL of the input fraction.

576
All samples were then boiled at 100°C for 12 minutes. 2.5 μL 10 mg/mL proteinase K 577 was then added, and the samples incubated at 55°C for 30 minutes. Samples were 578 then boiled at 100°C for 10 minutes, after which the beads and Chelex® were pelleted 579 and 60 μL of the supernatant transferred to a clean tube. The input and IP fractions 580 were diluted by 1:000 and 1:5, respectively. 2 μL of diluted input and IP template were  Table S5. genes marked by the specified histone modification (completely or < 1 bp association).

627
Read counts were imported into R, normalised and subject to differential expression The sequence data that support the findings of this study are available at the NCBI