1. Characterisation of a tlr2 mutant zebrafish line
The tlr2sa19423 mutant (tlr2-/-) carries a thymine to adenine point mutation that creates a premature stop codon (Fig. 1a), which is located in the sequence coding for the C-terminus of the leucine-rich repeat (LRR) domain. This leads to a truncated protein without the Toll/IL-1 receptor (TIR) domain, which is required for the interaction with Myd88 and Tirap (Mal) [46, 47].
To confirm whether tlr2 mutation blocks its downstream pathway, we analyzed the gene expression profiles of zebrafish treated with the TLR2 agonist, Pam3CSK4 similarly as in our previous work [48], now also including the heterozygote mutant (tlr2+/-). Pam3CSK4 was injected into the blood island of zebrafish embryos at 27 hours post fertilization (hpf). One hour after injection (hpi), we collected samples and performed qPCR to analyze the expression levels of CCAAT/enhancer-binding protein beta (cebpb) and FOS Like Antigen 1a (fosl1a), previously shown to be specific targets of Tlr2 signaling [45]. In wild-type siblings and tlr2+/ sa19423 heterozygotes (tlr2+/-), the expression levels of cebpb and fosl1a, as well as the inflammatory gene il1b were significantly induced upon Pam3CSK4 injection, whereas tlr2-/- showed no significant response (Fig. 1b-d). tlr2+/- animals showed a lower induction of these marker genes than the wild-types, indicating that there is an effect of the tlr2 mutation even in the heterozygote. To confirm that these results are specific to the Tlr2 pathway, we injected flagellin, a Tlr5 agonist, into 27 hpf embryos. Flagellin induced il1b expression, but not cebpb and fosl1a expression in wild-type siblings, tlr2+/- and tlr2-/- larvae (Fig. 1e-g). Overall, these data show that tlr2 mutation specifically blocks the Tlr2 downstream pathways.
To determine if immune cell development was affected by tlr2 mutation, we used the double-transgenic line tlr2+/+ Tg(mpeg1:mCherry-F);TgBAC(mpx:EGFP) and tlr2-/- Tg(mpeg1:mCherry-F);TgBAC(mpx:EGFP) to count the number of macrophages and neutrophils at 2 dpf. The results show that there is no significant difference in the number of macrophages and neutrophils at 2 dpf between the wild type and mutant (Fig. 1j, k).
2. Comparison of gene expression profiles of tlr2 homozygote and heterozygote mutants in the absence of infection
In order to investigate the systemic effects of the tlr2 mutation we compared basal levels of gene expression in the absence of infection between tlr2 homozygote and heterozygote mutants (Fig. 2). We chose to use the heterozygote tlr2+/- line as a control since these are genetically as comparable as possible, considering the fact that random ENU mutations could still contaminate our background even after 3 generations of outcrossing and the use of sibling lines. To further minimize the effect of polymorphisms, we pooled 10 larvae in each of our samples. These results show that there is a large group of genes that are expressed differently even at extremely stringent false discovery rate (FDR) adjusted p-value of 10-10 (Fig. 2). Since the results of the DEseq2 analyses showed such a surprisingly large number of significant differences we also used another statistical method for analyses, called edgeR, that differs in normalization and estimation of the dispersion parameters [49]. EdgeR analyses confirmed the statistical significance of the differences between the homozygote and heterozygote mutants (Fig. 2a and 2b). While no differences in numbers of mpeg-positive macrophages were detected (Fig. 1j, k), the RNAseq analysis showed that the mpeg1 gene is expressed approximately 2 fold higher in the tlr2+/- control than in the tlr2-/- mutant (Supplementary Table 1). This was confirmed by pixel count measurements in the transgenic line tlr2+/+ Tg(mpeg1:mCherry-F);TgBAC(mpx:EGFP) and tlr2-/- Tg(mpeg1:mCherry-F);TgBAC(mpx:EGFP) (Fig. S1 a, b). These results showed lower fluorescence of the mCherry mpeg1 reporter as compared to the wild-type siblings at 2 dpf (Fig. S1 a) whereas no difference was detected for the eGFP mpx reporter (Fig. S1 b). Therefore these in vivo results confirm the RNAseq data. We also examined the expression levels of the tlr2 gene in the homozygote and heterozygote mutants in the RNAseq data (Fig. S2). The results show that, although tlr2 is very lowly expressed, there is no difference in expression levels between the homozygote and heterozygote mutants, indicating that there is no non-sense mediated mRNA decay (NMD).
To investigate whether the basal expression level differences of genes might be relevant to the Tlr2 pathway, we performed GO analysis of a set of genes that were expressed differently with the very low FDR adjust P value of 10-10 (Supplementary Table 2). Gene ontology analyses of the gene sets that differ with an FDR of 10-10 in both analyses show that there is a large enrichment of genes belonging to the GO terms related to neural development (Supplementary Table 2). In addition, we have also analyzed differences with a fold change criterion of two and FDR of 0.05. These analyses indicated that under the GO category “Transcription factor genes” only two categories of genes including the c-Maf transcription factor were present (Supplementary Table 3). Pathway analysis shows that genes that function in glucose metabolism are differentially expressed in the tlr2 homozygote versus heterozygote mutant (Fig. S3).
3. Tlr2-specific gene expression profiles after Mm infection
To study the role of Tlr2 in Mm infection we tested the tlr2 mutant line in comparison with the heterozygote and wild-type sibling controls (Fig. 3). No differences in bacterial burden were observed at 3 dpi (Fig. 3d). However, we found that bacterial burden was significantly higher in tlr2-/- than in tlr2+/- and wild-type larvae at 4 dpi (Fig. 3e). There was no significant difference in bacterial infection burden between the heterozygote mutant strain (tlr2+/-) and the wild-type siblings (tlr2+/+). At 4 dpi, tlr2-/- larvae showed a significantly decreased percentage of survival compared with tlr2+/- and wild-type larvae (Fig. 3f). There was no significant difference in percentage of survival between tlr2+/- and wild-type larvae at 4 dpi. For further analysis of the infection phenotype we performed confocal laser scanning microscopy (CLSM) in the tlr2+/+ Tg(mpeg1:EGFP), tlr2+/- Tg(mpeg1:EGFP) and tlr2-/- Tg(mpeg1:EGFP) transgenic lines in which macrophages are fluorescently labelled (Fig. 4a, b and Fig. S4). The results using fluorescent automated pixel count analyses showed that in the tlr2 mutant the number of bacteria that were not present within macrophages was significantly higher than in the heterozygote and wild-type sibling controls (Fig. 4c). We also observed by manual counting that in the tlr2 mutant a significantly higher number of bacteria in large extracellular clusters (Fig. 4d), indicative for cording structures [50]. Manual counting of infected macrophage clusters showed a significantly lower number of granulomas in the tlr2 mutant (Fig. 4e). However, the percentage of larvae with at least one granuloma was not different from the heterozygote and wild-type sibling controls (Fig. 4f). These results show that the tlr2 mutation results in a defect in the defense response against mycobacteria at 4 dpi.
Next, we set out to assess the general inflammation and specific immune responses in tlr2 mutants. We analyzed these parameters at 4 dpi to correlate transcriptional responses with the first microscopically measured effect of the mutation on the progression of infection process as measured by bacterial burden (Fig. 3e). For this, genes that were shown previously to be specifically or aspecifically responding to Pam3CSK4 [45] were analyzed by qPCR of tlr2-/- and tlr2+/- larvae upon Mm (strain Mma20) infection at 4 dpi: tlr2-aspecific response genes il1b, tnfa, tnfb, irg1l, and tlr2-specific response genes fosl1a and cebpb. Our results show that the induction levels of il1b, tnfb, fosl1a and cebpb in tlr2-/- larvae were significantly reduced when compared to the heterozygotes in the infected condition (Fig. 5). The induction levels of irg1l and tnfa were less clearly affected. tlr2-/- larvae failed to upregulate fosl1a and cebpb expression in response to Mm administration (Fig. 5e, f). In our previous work, we showed that Cxcl11-like chemokines expressed in macrophages play a crucial role in granuloma formation upon Mm infection [51, 52]. We therefore conducted qPCR to assess the expression levels of cxcl11-like genes, previously shown to be induced by infection [51], including cxcl11aa and cxcl11ac (Fig. 5g, h). The expression levels of cxcl11aa and cxcl11ac was significantly higher in tlr2+/- larvae than in tlr2-/- upon Mm infection (Fig. 5g, h). These results indicate that tlr2 mutation results in a defective immune or inflammatory response to Mm infection.
To further study the function of tlr2 in defense against mycobacterial infection, we performed RNAseq of tlr2+/- and tlr2-/- larvae at 4 dpi with M. marinum strain (strain Mma20) and PBS as control. We summarized the number of differential expressed genes (DEGs) according to P-value (Fig. 6a, b) and in volcano plots (Fig. S5). The number of DEGs in tlr2+/- infected with Mm was higher than those in tlr2-/- at any given P-value or false discovery rate less than 0.05, or any given fold-change with a P-value less than 0.05. These data also show that most of the genes downregulated by the Mm infection in the control remained unchanged in tlr2 mutants (Fig. 6a, b). To further analyze these RNAseq data, we chose the genes with a threshold of a P-value less than 0.05 in tlr2+/- with Mm infection (1102 up- and 827 down-regulated genes, Fig. 6a). Then, for these genes, we calculated the fold-change ratio of tlr2+/- versus tlr2-/-, and genes with ratios greater than 2 or less than 0.5 were selected for further analysis (Fig. 6c). As a result, 97 and 92 genes were scored as tlr2 specific up- and down-regulated genes, respectively. Next, we conducted GO analysis (Fig. 6d, e) showing that genes grouped into the immune system category are the most prominently deregulated (36%) in the whole tlr2 up-regulated 97-gene set (Fig. 6d). Within this category we found genes involved in lysosome, chemotaxis, transcription regulation, diverse immunoglobulin domain-containing proteins (dicps) and other immune processes (Fig. 6d and 7a-e). For other categories, many up-regulated genes fell in the categories oxidation-reduction process, DNA repair, transcription regulation and apoptotic process regulation (Fig. S6). In the tlr2 down-regulated 92-gene set, the immune related genes also were the largest portion (15%; Fig. 6e, 7f). Many of these genes are poorly annotated and include genes encoding cysteine proteases, a nitr gene, a mafb transcription factor gene and hsp70. Categories encompassing non-immune related genes are listed in Fig. S7. These results show that, upon infection with Mm, tlr2 mutants show a dampened response of immune genes.
To show the relationship between DEGs, we constructed networks based on common expression targets in the 97 up- and 92 down-regulated genes [53]. The involved networks with the up-regulated genes contain lipc, nr0b2, hmgcr, itln1, fgf23, vdr, irak3 and tlr8 (Fig. S8). In tlr2+/- controls we observed positively regulated hmgcr, fgf23, vdr and tlr8, and negatively regulated nr0b2, whereas tlr2-/- showed an opposite regulation for most of these genes. However, lipc and irak3, were positively regulated in both tlr2+/- and tlr2-/-. For the 92 downregulated genes, a network containing nr0b2 and mafb was constructed (Fig. S9). nr0b2 and mafb were positively regulated in tlr2-/- zebrafish and down regulated in tlr2+/-.
4. Enrichment analysis of Tlr2 specific genes after Mm infection
To link our data to gene sets defined based on prior biological knowledge (including GO), we conducted Gene-Set Enrichment Analysis (GSEA) of the differently expressed genes. This method derives its power by focusing on gene sets, that is, groups of genes that share common biological function, chromosomal location, or regulation [54]. Because the number of predicted gene sets were too large for our analysis method (more than 1,000 with P < 0.05), we focused on the gene sets related to metabolic, immunological and inflammation pathways. GSEA predicted 61 pathways in tlr2+/- and 67 pathways in tlr2-/- zebrafish responsive to Mm infection (Supplementary Table 4A and 4B). Interestingly, most of these pathways are common in tlr2+/- and tlr2-/-, however some of them are detected as being anti-correlated in regulation, including some pathways underlying natural killer cell functions and omega-6-fatty acid metabolism (Supplementary Table 4C). We also performed Sub-Network Enrichment Analysis (SNEA) [55] to identify possible key genes that are responsible for the difference in response of the tlr2+/- and tlr2-/- group to Mm infection (P < 0.05). SNEA predicted 565 and 503 pathways for tlr2+/- and tlr2-/- zebrafish that are linked to the response to infection, respectively (Supplementary Table 4D and 4E), and 264 and 202 of them are specific for the response in tlr2+/- and tlr2-/- fish, respectively (Supplementary Table 4F). Since the RNAseq was conducted with the total RNA from whole body of zebrafish, it is not possible to define which pathways are involved in macrophage functions related to tlr2 expression. Therefore, we analyzed previously published DNA microarray data (GDS4781) of human macrophages transfected with Mtb from Gene Expression Omnibus [56]. SNEA analysis of human macrophages shows that 659 pathways are linked to Mtb infection (Supplementary Table 4G). By comparing the human SNEA result to tlr2+/- specific pathways in zebrafish, 56 pathways were defined as tlr2+/--specific (Supplementary Table 4H). Of these, the pathway of TLR8, which has the lowest P-value in both zebrafish and human enrichment (Supplementary Table 4H), and its network is depicted as example in Figure S10. Overall, these analyses show that tlr2-/- mutants have a strongly altered immune response of which many pathways that are linked to human tuberculosis are differently responding.