ORIGINAL RESEARCH article
RNAseq Profiling of Leukocyte Populations in Zebrafish Larvae Reveals a cxcl11 Chemokine Gene as a Marker of Macrophage Polarization During Mycobacterial Infection
- 1Institute of Biology Leiden, Leiden University, Leiden, Netherlands
- 2ZF-screens B. V., Leiden, Netherlands
Macrophages are phagocytic cells from the innate immune system, which forms the first line of host defense against invading pathogens. These highly dynamic immune cells can adopt specific functional phenotypes, with the pro-inflammatory M1 and anti-inflammatory M2 polarization states as the two extremes. Recently, the process of macrophage polarization during inflammation has been visualized by real time imaging in larvae of the zebrafish. This model organism has also become widely used to study macrophage responses to microbial pathogens. To support the increasing use of zebrafish in macrophage biology, we set out to determine the complete transcriptome of zebrafish larval macrophages. We studied the specificity of the macrophage signature compared with other larval immune cells and the macrophage-specific expression changes upon infection. We made use of the well-established mpeg1, mpx, and lck fluorescent reporter lines to sort and sequence the transcriptome of larval macrophages, neutrophils, and lymphoid progenitor cells, respectively. Our results provide a complete dataset of genes expressed in these different immune cell types and highlight their similarities and differences. Major differences between the macrophage and neutrophil signatures were found within the families of proteinases. Furthermore, expression of genes involved in antigen presentation and processing was specifically detected in macrophages, while lymphoid progenitors showed expression of genes involved in macrophage activation. Comparison with datasets of in vitro polarized human macrophages revealed that zebrafish macrophages express a strongly homologous gene set, comprising both M1 and M2 markers. Furthermore, transcriptome analysis of low numbers of macrophages infected by the intracellular pathogen Mycobacterium marinum revealed that infected macrophages change their transcriptomic response by downregulation of M2-associated genes and overexpression of specific M1-associated genes. Among the infection-induced genes, a homolog of the human CXCL11 chemokine gene, cxcl11aa, stood out as the most strongly overexpressed M1 marker. Upregulation of cxcl11aa in Mycobacterium-infected macrophages was found to require the function of Myd88, a critical adaptor molecule in the Toll-like and interleukin 1 receptor pathways that are central to pathogen recognition and activation of the innate immune response. Altogether, our data provide a valuable data mining resource to support infection and inflammation research in the zebrafish model.
Macrophages are phagocytic innate immune cells that, together with neutrophils, form the cellular arm of the first line of host defense against invading pathogens. The activation of macrophages is initiated by the recognition of microbial and danger signals by Pattern Recognition Receptors (PRRs), such as the Toll-like receptors (TLRs), which recruit MYD88 and related adaptor molecules for signal transduction to MAP kinases and Nuclear Factor κB (NFκB) (1). The signaling pathways downstream of TLRs and other PRRs regulate the transcription of a large number of genes that are involved in signaling between immune cells (cytokine and chemokine genes) and in host defense (1, 2). To exert their anti-microbial function, macrophages employ several mechanisms, such as the production of reactive oxygen and nitrogen species, the production of antimicrobial peptides, and the degradation of microbes through the phagosomal-lysosomal and autophagy pathways (3). Following successful elimination of microbial invaders, macrophages mediate the resolution of inflammation by clearing cellular debris and eliminating the surplus immune cells that are undergoing cell death at the infection site (4). In addition to these primary functions in infection and inflammation, macrophages orchestrate a range of developmental processes. For example, macrophages interact with endothelial cells to support angiogenesis (5, 6), help control definitive hematopoiesis (7), and facilitate electrical conduction in the heart (8). Thus, macrophages are highly versatile cells, not only functioning as central players in the immune system, but also contributing to organismal development and maintenance of homeostasis.
Macrophages can adopt different states of activation, which are classically divided into a pro-inflammatory M1 state and an anti-inflammatory M2 state (9). These states are characterized by distinct cytokine and chemokine expression patterns as well as different metabolic profiles. However, it is clear that many intermediate phenotypes exist and that the distinction between M1 and M2 states is a simplification of how macrophage polarization occurs in different organs and tissues and in response to different stress signals (9). Macrophage polarization plays a major role in the context of disease. Tumor-associated macrophages can acquire anti-tumor or tumor-promoting phenotypes (10), macrophage metabolism is critical in development of atherosclerosis and other cardiovascular diseases (11, 12), and certain pathogens are known to manipulate the macrophage phenotype to their advantage (13–15). Therefore, a better understanding of macrophage polarization and function could open novel therapeutic avenues for diseases related to dysfunction or hyperactivation of this cell type.
The majority of studies on differentiation of macrophage subtypes have been performed in vitro, but recently it has been achieved to image the process of macrophage polarization during inflammation in a living organism (16). To this end, the optically transparent early life stages of the zebrafish were used (embryos and larvae), expressing different fluorescent markers for the macrophage cell type and for a classical M1 marker, tumor necrosis factor alpha (tnfa). Expression of the tnfa reporter was observed in macrophages recruited to sites of injury or sites of Escherichia coli infection. Furthermore, the tnfa-positive macrophages were shown to express other M1 markers (il1b, il6), while tnfa-negative macrophages expressed M2 markers (tgfb1, ccr2, cxcr4b). Dynamic tracing of reporter expression showed that tnfa-positive cells at inflammation sites reverted to a tnfa-negative phenotype during the resolution phase (16). In a follow-up study, a tail fin amputation model was used to show that early recruitment of a tnfa-expressing macrophage subpopulation is required for blastema formation and subsequent fin regeneration (17). Zebrafish models for a wide variety of human diseases have been developed in the recent years (18). Therefore, the tnfa-reporter together with other M1 and M2 lines that are still to be developed, will find many applications to elucidate the functions of polarized macrophage subsets during disease processes.
The most frequently used promoter for driving macrophage-specific expression of fluorescent reporters in zebrafish is that of the macrophage expressed 1 gene (mpeg1.1, hereafter called mpeg1) (19, 20). The mpeg1 gene codes for a perforin-like protein with anti-bacterial function (21). Fluorescent mpeg1 reporter lines have been used to study a diverse range of processes. These include for example, macrophage-endothelial interactions (22), long-distance communication between macrophages and pigment cells (23), the function of tumor-associated macrophages (24), and the role of macrophages in host defense against infections (25). Fluorescent mpeg1 reporter expression in zebrafish embryos and larvae marks the monocytic precursors of macrophages in the blood as well as tissue-resident macrophages, including microglia in the developing brain (19, 26). The expression of mpeg1 reporters is non-overlapping with that of a neutrophil-specific BAC reporter line driven by the myeloperoxidase (mpx) promoter (27), or with a reporter for immature lymphocytes controlled by the promoter of the LCK proto-oncogene, Src family tyrosine kinase (lck) gene (28).
The generation of macrophage and neutrophil reporter lines has gone together with the development of zebrafish models for a variety of human infectious diseases, providing new possibilities to visualize host-pathogen interactions in real time (25, 29). It has been shown that zebrafish embryos rely on macrophages for an effective host defense against different pathogens, such as Staphylococcus aureus (30) and Salmonella enterica servovar Typhimurium (31). In contrast, the ablation of macrophages was found to protect zebrafish embryos against infection with Burkholderia cenocepacia, revealing that macrophages are critical for the proliferation of this pathogen and for the development of a fatal inflammatory response (32). Macrophages play a dual role during infection with Mycobacterium marinum, a pathogen widely used to model human tuberculosis in zebrafish, since it provides access to the early stages of tuberculous granuloma formation that is initiated by the aggregation of infected macrophages (33–35). On the one hand, abundant extracellular growth is observed in macrophage-deficient hosts, indicating that proliferation of M. marinum is restricted when phagocytosed by macrophages (36). On the other hand, infected macrophages, driven by bacterial virulence mechanisms, can migrate into tissues and recruit new macrophages, which promotes the cell-to-cell spreading of M. marinum and the expansion of granulomas (37, 38). Consequently, a mutation in the macrophage-specific chemokine receptor cxcr3.2, which restricts macrophage motility, has a positive outcome on the ability of zebrafish embryos to control M. marinum infection (39).
Despite extensive use of the zebrafish mpeg1 reporter lines for studying macrophage biology, the expression signature of these cells has remained uncharacterized. Here, we isolated mpeg1 expressing cells from transgenic zebrafish larvae by fluorescence activated cell sorting (FACS) and performed RNAseq to investigate what distinguishes the mpeg1-driven expression profile from the signatures of neutrophil and lymphocyte populations isolated from mpx and lck reporter lines. In addition, we determined a core expression set of 744 zebrafish macrophage markers, based on enriched expression in mpeg1-positive cells. We compared this gene set with published RNAseq profiles of human macrophages differentiated in vitro toward M1 or M2 phenotype (40), which showed that zebrafish macrophages express a mixed profile of M1 and M2 markers under unchallenged conditions. We then studied how the expression profile is changed upon M. marinum infection and identified a homolog of the human M1 marker CXCL11 as a robust and specific marker of Mycobacterium-infected macrophages that is induced by myd88-dependent signaling.
Materials and Methods
Zebrafish Husbandry and Infection Experiments
Zebrafish lines in this study were handled in compliance with local animal welfare regulations as overseen by the Animal Welfare Body of Leiden University (License number: 10612) and maintained according to standard protocols (zfin.org). All protocols adhered to the international guidelines specified by the EU Animal Protection Directive 2010/63/EU. All experiments with these zebrafish lines were done on larvae before the free-feeding stage. Zebrafish lines included AB/TL, Tg (mpx:eGFP)i114 (27), Tg (mpeg1:Gal4-VP16)gl24/ (UAS-E1b:Kaede)s1999t (19), Tg (mpeg1:mCherry-F)ump2 (20), Tg (lck:eGFP)cz2 (28), cxcr3.2−/−hu6044, and their cxcr3.2+/+ wildtype controls (39), myd88−/−hu3568 and their myd88+/+ wildtype controls (41). Embryos were grown at 28.5°C in egg water (60 μg/ml Instant Ocean sea salts). Mycobacterium marinum M or its RD1-deficient (ΔRD1) isogenic strain (42) containing pSMT3-mCherry (43) was grown and prepared for injections as described (44) and microinjected into the caudal vein of embryos at 28 h post fertilization (hpf) using, where not differently specified, a dose of 100–125 colony-forming units (cfu) per embryo. After injection, embryos were transferred into fresh egg water and incubated at 28°C for 4 or 5 days before collection. Proper infection was controlled by fluorescent imaging before embryo dissociation.
Embryo Dissociation and Fluorescent Activated Cell Sorting (FACS)
Immune cells from 5 to 6 dpf larvae were isolated by FACS as described previously (45). Briefly, live embryos were rinsed in calcium free Ringer solution for 15 min and then digested with 0.25% trypsin for 60–90 min at 28°C. Digestion was stopped with 1 mM CaCl2 and 10% fetal calf serum and the cell suspension was centrifuged and washed in PBS. The cell pellet was resuspended in 1–2 ml in Leibovitz medium L15 without phenol-red with 1% fetal calf serum, 0.8 mM CaCl2, 50 units/ml penicillin and 0.05 mg/ml streptomycin and filtered through a 40 μm cell strainer. FACS was performed at 4°C using a FACS AriaIII (BD Biosciences) with the BD FACSDiva software (version 6.1.3). For collecting mCherry-positive cells a Coherent Sapphire solid-state 561 nm yellow green laser with 36 mW power was used. Laser settings applied were 600 LP, 615/20 BP. For sorting eGFP and Kaede positive cells a Coherent Sapphire solid-state 488 nm laser with 15.4 mW power was used. Laser settings applied were 505 LP, 530/30 BP. mCherry and eGFP/Kaede gates were set-up with non-fluorescent samples and allowed to collect an average of, respectively, 33.8 +/– 16.4 mCherry and 11.6 +/– 4.4 eGFP/Kaede false-positive cells per million of sorted cells. An average of 526 mpeg1:Kaede, 195 mpx:eGFP and 983 lck:eGFP positive cells in 5 dpf embryos and 1,826 mpeg1:Kaede and 5,482 mpeg1:mCherry positive cells at 6 dpf were collected per million of sorted cells. For background expression assessment 500,000 non-fluorescent cells were sorted for each sample. Cell fractions were separately collected in supplemented L15 medium and RNA isolation was performed directly after sample collection.
RNA Isolation, Illumina Sequencing, and Real Time PCRs
RNA extraction was done using the RNAqueous-Micro Kit (Ambion) or RNeasy mini kit (Qiagen). Quality of RNA used for Illumina sequencing was checked on an Agilent Bioanalyzer 2100 using the RNA 6000 Pico kit (Agilent). RNA sample integrity was assessed based on RIN number and visual inspection of the electropherograms. cDNA synthesis and amplification was performed using the SMARTer Ultra Low RNA Kit for Illumina sequencing (Clontech) according to the manufacturer's protocol. Illumina TruSeq DNA Sample Preparation Kit v2 (Illumina) was used on sheared cDNA to prepare libraries. Three changes were made to manufacturer's protocol: the adapters were diluted 20 times in the adapter ligation step, library size selection was achieved by double Ampure XP purification with a 0.7x beads to library ratio and library amplification was made with 15 cycles. The resulting libraries were sequenced using an Illumina HiSeq2000 with 50 bp paired-end reads for all samples and single-end reads for 6 dpf samples. RNA collected for real time PCR experiments was further purified using column DNA digestion (RNase-Free DNase set, Qiagen). cDNA was prepared using iScript cDNA-synthesis kit (Invitrogen, Life Technologies) and was used as a template for qRT-PCR reaction with iQ SYBR Green Supermix according to the manufacturer's instructions (Bio-Rad Laboratories). Gene expressions were analyzed using the ΔΔCt method and were normalized against ppiab for whole-mount analyses and to eif5 for FACS sorted cells. Primers for these genes were: cxcl11aaFw: ACTCAACATGGTGAAGCCAGTGCT; cxcl11aaRv: CTTCAGCGTGGCTATGACTTCCAT; ppiabFw: ACACTGAAACACGGAGGCAAAG; ppiabRv: CATCCACAACCTTCCCGAACAC; eif5Fw: CAAGTTTGTGCTGTGTCCCG; eif5Rv: AGCCTTGCAGGAGTTTCCAA; ccr2Fw: TGGCAACGCAAAGGCTTTCAGTGA; ccr2Rv: TCAGCTAGGGCTAGGTTGAAGAG; ccr5(ccr12b.2)Fw: GGCTTCCAACATCATCTTCACCCTCAC; ccr5Rv: CTATCATCCGAGTGCGCATGATGG; cxcr3.2Fw: CCTCTGTTGGTAATGCTGTATTGC; cxcr3.2Rv: ACACGATGACTAAGGAGATGA; tfebFw: GCATTACATGCAGCATCGCATGCC; tfebRv: CGTGTACACATCCAAATGACTGCTGG; tfecFw: AACAGTACCTCGCTTTGGGC; tfecRv: CAGTGTTCCCAGCTCCTTGA; ctsl.1Fw: CTGGAGGGACAAGGGCTATG; ctsl.1Rv: CTATGGCAACAGATATGGGGCC; spi1Fw: CCGATGAGGAGTGTATGAGAG; spi1Rv: GCTTGGATGAGAACTGGAATG. Primers for mpeg1 and marco were described in Zakrzewska et al. (46), primers for mpeg1.2 in Benard et al. (21), primers for mpx in Kanwal et al. (47), primers for apoa4 in Ordas et al. (48), primers for mmp9 in Stockhammer et al. (49), and primers for grn1/2, grna, and grnb in Solchenberger et al. (50).
RNA-Sequencing on 20 Sorted Cells
Dissociation and cell sorting of infected embryos were performed as mentioned previously. Twenty cells were sorted directly in cDNA synthesis buffer from the SMARTer Ultra Low RNA Kit for Illumina sequencing (Clontech) and used directly for cDNA synthesis. The resulting cDNA were amplified for 20 cycles and used for library preparation and single-end Illumina sequencing as mentioned above.
RNA-Seq Data Analysis
Image analysis, base calling, and adapter removal were done using the Illumina HCS version 1.15.1. Sequencing depth of all samples was between 6.5 and 29.8 million reads with an average of 14.9 million reads. All reads were aligned to the Ensembl zebrafish genome (Zv9) using Bowtie2 with the option –very-sensitive. On average, 72.5% of the total reads aligned on the genome, with a minimum of 32.3% and a maximum of 91.3%. Aligned reads were mapped to zebrafish transcripts using TopHat with default parameters except for minimal intron length set up to 2 and library type set up to fr-unstranded. A modified version of the Ensembl Zv9_79 annotation with additional manually annotated genes was used (Table S1). Differential expression analyses between non-fluorescent cells and fluorescent positive cells were performed using DEseq package in R and DEseq2 for the 20 cell samples. PCA plots and Pearson correlation HeatMaps were generated with DEseq package build in functions. Networks based on GO-enrichment analysis (GOEA) were produced using BiNGO and EnrichmentMap in Cytoscape (51). Gene Set Enrichment Analyses were performed with the GSEA software from the Broad Institute (52, 53) version 3.0.
The Macrophage-Specific Transcriptome of Zebrafish Larvae
To determine the gene expression signature of macrophages, we used mpeg1-driven reporters expressing Kaede or mCherry fluorescent proteins (19, 20, 25). RNA was extracted from positive and negative fluorescent cell fractions obtained by FACS sorting of single cell suspensions obtained by dissociating 5 or 6 dpf transgenic larvae. Illumina sequencing (RNAseq) of RNA samples from Kaede-labeled macrophages at 5 dpf and 6 dpf and mCherry-labeled macrophages at 6 dpf, each in duplicate, resulted in a total of 6 replicates. Reproducibility between these replicates after alignment and mapping of the reads was high, as shown by calculation of the Pearson correlation coefficient (Figure S1). Principal component analysis (PCA) showed that all macrophage samples segregated clearly from the samples of negative fluorescent cell fractions. Furthermore, PCA indicated minor differences between samples from macrophages with Kaede and mCherry markers and between the Kaede-labeled macrophages at 5 and 6 dpf (Figure S2A). Between 8,177 and 12,162 genes were expressed (Transcripts Per Million (TPM) ≥ 3) in the macrophage populations (Tables S2, S3) and a total of 8,892 genes were shared in at least two of the three conditions. Together, these results indicate that our protocol of RNAseq on FACS sorted cells from zebrafish larvae produces high quality results.
We performed differential expression (DE) analysis on the duplicates from the three different conditions by comparing results from fluorescence positive cell fractions with the related fluorescent negative cell fractions. By selecting an adjusted p-value threshold of 0.01, we detected a similar number of genes expressed specifically in Kaede-labeled macrophages at 5 and 6 dpf (400 and 414 genes, respectively), whereas more genes were detected in the 6 dpf mCherry-labeled macrophages (1,318 genes). Comparison between the different conditions showed a high overlap between the enriched gene sets (Figure 1A). To produce a complete and accurate description of the macrophage transcriptome in the larvae, we selected the genes that met the significance threshold in at least two out of the three conditions. This dataset, hereafter named zebrafish macrophage (zfM) core expression set, contains 437 genes (Figure 1A; Table S4).
Figure 1. Core expression of macrophages in zebrafish larvae. (A) Overlap of the genes enriched (P-adj < 0.01) in zebrafish larval macrophages (zfM) at 5 or 6 dpf from mpeg1-driven Kaede and mCherry reporter fish. The zfM core expression is defined as the overlap of two or three mpeg1 positive populations and includes 437 genes. (B) Expression table of immune cell specific genes. First, second and third columns correspond, respectively, to macrophage expression from 5dpf mpeg1:Kaede, 6dpf mpeg1:Kaede, and 6dpf mpeg1:mCherry reporter larvae. Colored cells correspond to genes enriched in the corresponding sequencing data (log2 (fold change) > 1, P-adj < 0.01) whereas gray cells correspond to non-enriched genes. Numbers are expression levels expressed in TPM in fluorescence positive cells. (C) Network visualization of GO enrichment analysis of genes from the core macrophage expression data set using BiNGO and EnrichmentMap. Magenta nodes are terms found exclusively in the zfM core dataset and green nodes are found in both zfM core and neutrophil data sets. Node size corresponds to the number of genes associated to the enriched GO term and edge size to the similarity coefficient between two nodes.
The zfM core dataset includes the main genes that are known to be specific for macrophages and myeloid cells in zebrafish larvae (Figure 1B). For example, in addition to mpeg1 itself, the macrophage-specific genes csf1ra, mhc2dab, the myeloid genes spi1a and b, and the pan-leukocyte markers coro1a, ptprc, and ptpn6 were detected (25, 46). In addition, we validated the macrophage-enriched expression of several chemokine receptors (ccr2, ccr5 also known as ccr12b.2, cxcr3.2) and lysosomal genes (tfeb, tfec, ctsl.1) by RT-qPCR (Figure S3A). Network visualization of Gene Ontology (GO) terms revealed enrichment of biological processes linked to immune response and inflammation, antigen processing and presentation, signal transduction, response to bacterium, peptidoglycan catabolic process, chemotaxis, and proteolysis (Figure 1C). Similarly, GO terms for molecular function were clearly linked to immune cell function and defense mechanisms (Figure S4).
Among the genes from the zfM dataset, 31% (137 out of 437) corresponded to uncharacterized proteins or non-coding RNAs. Manual annotation showed that many of these uncharacterized sequences belong to large immune-related protein families, including the immunoglobulins (36), the C-type lectins (13), and NACHT/LRR proteins (8) (Table S4). Other uncharacterized genes were also related to immunity, such as genes coding for proteins with chemokine/interleukin-like domains, chemokine receptor like domains, interleukin receptor-like domains, complement domains, leukotriene receptor like domains, or MHC class II alpha and beta chains. In addition, 17 genes correspond to non-coding RNAs, long-intronic-non-coding RNAs or processed transcripts, of which the possible role in immunity is of interest for further study.
Comparison of Zebrafish Macrophage and Neutrophil Expression
For comparison with the macrophage transcriptome, we studied the neutrophil transcriptome by sequencing the fluorescent cell population extracted from 5 dpf mpx:gfp larvae (two replicates). A total of 8,627 genes were found expressed in neutrophils (TPM ≥ 3, Tables S2, S3). Selection of differentially expressed genes revealed a data set composed of 227 neutrophil-enriched genes (P-adj < 0.01) (Table S5). Among these genes, 111 (49%) are shared with the zfM core expression dataset (Figure 2A). The neutrophil markers lyz and mpx were detected enriched at a high level in neutrophil population, although a low level of these transcripts could be detected in zebrafish macrophages (Figure 2A) as well as in human macrophages (40). Several macrophage markers could be detected in neutrophils, but none were significantly enriched (Figure 2A). The macrophage and neutrophil specific signatures of several genes were validated by RT-qPCR (Figure S3B).
Figure 2. Comparison of macrophage core dataset with neutrophils and lymphoid cells expression sets. (A) Overlap of the genes from the zfM core data set (magenta) and the genes enriched [log2 (fold change) > 1, P-adj < 0.01] when comparing mpx:gfp positive and negative cells from 5 dpf transgenic fish (green). Below is represented an expression table of a selection of immune cell specific genes. Magenta cells and green cells represent genes enriched in macrophages or neutrophils, respectively, whereas gray cells represent non-enriched genes. Numbers are expression levels expressed in TPM in fluorescence positive cells. (B) Network visualization of GO enrichment analysis of genes from the neutrophil expression data set using BiNGO and EnrichmentMap. Green nodes are terms found exclusively in neutrophil expression dataset and magenta nodes are found in both core macrophage and neutrophil data sets. Network legend is similar to Figure 1C. (C) Overlap of the genes from the myeloid data set corresponding to genes found either in the core macrophage or the neutrophil expression data sets (magenta) and the genes enriched (log2 (fold change) > 1, P-adj < 0.01) when comparing lck:gfp positive and negative cells in 5 dpf lck:gfp transgenic fish (blue). Below is represented an expression table of a selection of immune cell specific genes. Magenta cells and blue cells represent genes enriched in myeloid or lymphoid cell populations, respectively, whereas gray cells represent non-enriched genes. Numbers are expression levels expressed in TPM in fluorescence positive cells (for myeloid cells, expression levels in macrophage and neutrophil populations were averaged). (D) Network visualization of GO enrichment analysis of genes from the lymphoid expression data set using BiNGO and EnrichmentMap. Blue nodes are terms found exclusively in lymphoid expression dataset and magenta nodes are found in both myeloid and lymphoid data sets. Network legend is similar to Figure 1C.
GO analysis identified different biological processes specific for neutrophils grouped into protein metabolic process and lipid metabolic process (Figure 2B). GO terms associated to signal transduction, antigen processing and presentation, peptidoglycan catabolic process, and chemotaxis are enriched in macrophages but not in neutrophils, confirming the functional differences between the two myeloid lineages (Figure 1B). Many GO terms are shared between the two cell populations. However, their contents often are composed of different protein families (Figure 3). For example, proteolysis appears to be a major group in both cellular lineages, but macrophages express cathepsin coding genes (ctsc, ctsh, ctssb2, ctsz) whereas neutrophils express proteinases from the carboxypeptidase (cpa5), the elastase (ela2, ela2l, and ela3l), the chymase families as well as trypsin.
Figure 3. Different gene sets expressed in macrophage, neutrophil, and lymphoid cell population. Expression table of selected genes belonging to different GO terms. First, second and third columns correspond to macrophage, neutrophil, and lymphoid cell populations, respectively. Magenta, green, and blue cells represent genes enriched [log2 (fold change) > 1, P-adj < 0.01] in macrophages, neutrophils, or lymphoid cells compared to respective fluorescent-negative background cells. Numbers are expression levels expressed in TPMin the fluorescence positive cell fractions. Gray cells are non- enriched genes (P-adj ≥ 0.01).
Comparison of Zebrafish Myeloid and Progenitor Lymphoid Expression
While it is clear that zebrafish larvae do not yet have a functional adaptive immune system (54, 55) the precise state of development of both the innate and (immature) adaptive immune system remains unknown. We therefore used the lck:GFP transgenic line to study the transcriptome of immature lymphoid cells that are present in the thymus at 5 dpf. At this stage, T-lymphopoiesis is independent from hematopoietic stem cells (HSC) and arises from aorta endothelium, whereas HSC-dependent lymphopoiesis begins from 8 dpf onwards (56). RNAseq analysis showed that 8,547 genes are expressed in this cell population (TPM ≥ 3, Tables S2, S3) and 1,328 genes were enriched in the lymphoid cell progenitor compared to non-fluorescent cell background (Table S6). Comparison of genes enriched in lymphoid cells with genes specific for myeloid cells (i.e., genes enriched in either macrophages or neutrophils) showed a small overlap between the two populations (Figure 2C). None of the macrophage or neutrophil specific markers was detected in the lymphoid cell transcriptome. On the opposite, several known lymphocyte markers were detected only in this cell population (Figure 2C). GO analysis showed also little overlap with processes detected in myeloid cells (Figure 2D). Surprisingly, very few terms were associated to immune function, except for chemotaxis and response to organic substances, also shared with the myeloid lineage. However, these groups are composed of different genes. For example, the proteases expressed by lymphoid cells mainly belong to the proteasome (Figure 3). A closer look to immunity-related genes showed the presence of chemokines and chemokine receptors as well as MHC class II genes, but often expressed at a low level compared to myeloid cells (Figure 3). The main enriched GO terms in the lymphoid cells were found to be associated with metabolic processes, and cell cycle (Figure 2D), which might reflect the immature status of this cell population in developing zebrafish larvae.
Similarity Between the Zebrafish Macrophage Transcriptome and Human Polarized Macrophage Transcriptomes
By real time imaging of macrophages in a dual fluorescent mpeg1 and tnfa reporter line evidence has been obtained that zebrafish larvae differentiate M1 and M2 like polarized macrophages in response to wounding and infection (16). We found that the known M1 (il1b, tnfa/tnfb) and M2 (cxcr4b, il10, ccr2) markers were expressed in the zfM core expression set. Additionally, the zebrafish homologs of human M1-markers CXCL11 (cxcl11aa), MMP9 (mmp9) and TNFRSF1B (tnfrsf1b) and the M2-markers ALOX5AP (alox5ap), MARCO (marco), and TGFB1 (tgfb1b) were also detected in our zfM core dataset (Table S4).
To investigate further the similarities with human macrophages, we compared our transcriptomic data with RNA sequencing data published by Beyer et al. (40), in which transcriptomes of in vitro polarized M1 and M2 macrophages were analyzed.
By using Biomart (http://www.ensembl.org/) combined with custom annotations, we retrieved 7,963 human homologs of the genes expressed in mpeg1-positive macrophages. Approximately three quarters of these homologs were also found among the genes expressed in human M1 cells (70.9%) and among the genes expressed in M2 cells (72.1%) (TPM ≥ 3) (Figure 4A). Similar proportions of the human homologs of the zfM core expression set (for which 294 human homologs were identified, see Table S7) were found in human M1 (69.0%) and M2 (68.4%) macrophages (Figure 4A). Among those genes, only 4 were expressed exclusively in human M1 macrophages and 7 exclusively in M2 macrophages. These observations suggest that zebrafish macrophages are composed of a mixed population of M1 and M2 type macrophages.
Figure 4. Zebrafish larval macrophages have a gene signature similar to human M1 and M2 polarized macrophages. (A) Overlap of the genes detected in human M1 (blue) or M2 (green) polarized macrophages and of the human homologs of the zebrafish core macrophage data set (magenta). Black number correspond to the comparison between genes expressed in human or zebrafish cells (TPM ≥ 3) and red numbers correspond to the comparison of gene expressed in human cells and specifically enriched in zebrafish macrophages [log2 (fold change) ≥ 1, P-adj < 0.01]. (B) Network visualization of GO enrichment analysis of human homologs of zebrafish macrophage enriched genes not detected in the dataset of human M1 and M2 in vitro polarized macrophages using BiNGO and EnrichmentMap. Red nodes represent GO terms. Network legend is similar to Figure 1C. (C) Volcano plot showing the P-value (-log10-transformed) as a function of the fold-change (log2-transformed) between human M1 and M2 gene expression level of the gene set from Beyer et al. (40). Red dots are genes with a human homolog detected in the zebrafish macrophage core dataset. (D) Gene Set Enrichment Analysis (GSEA) plots of gene expression changes in human M1 in vitro polarized macrophages compared to human M2 in vitro polarized macrophages from Beyer et al. (40). Gene sets used for the analyses are genes expressed in zebrafish macrophages (TPM ≥ 3) (top) and genes from the zebrafish macrophage core dataset (lower).
A total of 87 homologs from the zfM core expression dataset were not present in the M1 or M2 polarized human macrophage datasets (Figure 4A). Among these were the known M1 marker Interleukin 12B and the M2 marker mannose receptor C type 1 (MRC1). Other genes detected exclusively in the zfM expression set were associated to the molecular function cell adhesion and differentiation, catabolic processes, and response to bacterium (Figure 4B). Genes coding for the peptidoglycan recognition protein Pglyrp family, and cytokine receptors (Ccr9a, Il22ra2) are present in these categories. The absence of these genes in the human M1 and M2 sets might be due to low expression levels in in vitro cultured cells or inaccurate orthology detected in zebrafish.
Finally, we computed the differential expression between human M1 and M2 polarized macrophages and searched whether the genes from the zfM core expression set were preferentially associated to either M1 (log2FC > 1) or M2 (log2FC < −1) signal. The results show that 48 genes from the zfM core dataset were associated to M1-enriched genes and 54 were associated to M2-enriched genes (Figure 4C).
We also used Gene Set Enrichment Analysis (GSEA) to compare the set of genes expressed (TPM ≥ 3) in zebrafish macrophages with the differential expression between human M1 and M2 polarized macrophages. The analysis showed a preference for M2-enriched genes, although this enrichment was not significant (FDR > 0.05) (Figure 4D, top). Focusing on the zfM core gene set also showed no clear enrichment for either M1- or M2-enriched genes (Figure 4D, lower).
Altogether, our results indicate no clear polarization of the zebrafish macrophages, suggesting the presence of both M1 and M2-typed macrophages in unchallenged larvae.
Effect of M. marinum Infection on the Zebrafish Macrophage Transcriptome Profile
As zebrafish larval macrophages display mixed M1 and M2 characteristics, we tried to induce a shift in activation phenotype by infecting mpeg1:mCherry embryos at 1 dpf with GFP-labeled Mycobacterium marinum (Mm), an intracellular pathogen of macrophages. Transcriptomes of infected and uninfected macrophages were profiled 5 day post infection (6 dpf).
When retrieving samples from infected larvae, only a small number of double positive cells were collected over a long sorting period, inducing variation between replicates. To minimize the differences, we reduced the sorting time and the number of steps in the protocol by collecting 20 infected and uninfected cells from Mm-infected larvae directly in cDNA synthesis buffer and by proceeding immediately to cDNA synthesis and amplification without RNA extraction. These modifications of the protocol led to reproducible results (Figures S2, S3). Differential expression analysis between infected and uninfected macrophages identified 330 upregulated and 139 downregulated genes (P-adj < 0.05) (Figure 5A; Table S8). GO analysis identified two terms enriched in the downregulated genes only: cell cycle and blood vessel morphogenesis. This group, often related to immune cell migration, included the genes flt1/VEGFR1, known to be expressed in mouse M2 macrophages in vitro (57, 58), and ptprja/CD148, expressed by human macrophages under exposure to LPS and other TLR-ligands but repressed under CSF-1 treatment (59). Performing GO analysis on the human homologs of this set of genes identified the terms transcription coactivator activity, NADP or NADPH binding and serine hydrolase activity associated to upregulated genes (Figure 5B) and protein localization to downregulated genes.
Figure 5. M. marinum infected macrophages exhibit a change in gene expression toward M1-polarization. (A) MA-plot showing the fold change (log2-transformed) between gene expression in Mm-infected and non-infected mpeg1:gfp positive cells from a 6 dpf embryos 5 days after injection of 100–125 cfu of Mycobacterium marinum M strain containing pSMT3-mCherry as a function of the normalized average count between the two conditions (log10-transformed) as calculated with DEseq2. Turquoise: log2FC ≥ 1 and P-adj < 0.05, red: log2FC ≤ −1, and P-adj < 0.05. (B) Network visualization of GO enrichment analysis of human homologs of up-regulated genes in infected macrophages compared with uninfected macrophages from infected larvae using BiNGO and EnrichmentMap. Network legend is similar to Figure 1C. (C) Volcano plot showing the P-value (–log10-transformed) as a function of the fold-change (log2-transformed) between human M1 and M2 gene expression level of the gene set from Beyer et al. (40). Turquoise and Red dots are genes with a zebrafish homolog respectively up- and down-regulated in the infected macrophages compared with the non-infected macrophages from M. marinum infected larvae. (D) Gene Set Enrichment Analysis (GSEA) plots of gene expression changes in human M1 in vitro polarized macrophages compared to human M2 in vitro polarized macrophages from Beyer et al. (40). Gene sets used for the analyses are human homologs of genes found up-regulated (log2FC ≥ 1 and P-adj < 0.05) (top) and down-regulated (log2FC ≤ −1 and P-adj < 0.05) (lower) in macrophages upon Mm infection as described in (A). (E) Table presenting the zebrafish genes expressed in M1 and M2 macrophages studied in Nguyen-Chi et al. (16). First column indicates their presence in our zfM core dataset. Second column indicates their enrichment (log2Fold Change) in Mm-infected macrophages compared to uninfected macrophages. Turquoise: log2FC ≥ 1 and P-adj < 0.05, light turquoise: log2FC ≥ 0, red: log2FC ≤ 0.
Our analysis revealed several Mm-induced genes that could play important roles in host defense. These include for example CIITA, the master transactivator of MHC class II gene expression, which has previously been described to be important for limiting M. tuberculosis infection in mice (60). Another Mm-induced gene is the mpeg1-family gene mpeg1.2, which we have previously shown also to be inducible by Salmonella infection (21). The mpeg1 genes encode proteins of the perforin family with proposed anti-bacterial functions in macrophages that require further mechanistic dissection (21). On the other hand, other overexpressed genes could be more beneficial for the survival of the bacteria. The gene nsfb, the zebrafish homolog of the human N-ethylmaleimide sensitive factor, has been proposed to promote the fusion of phagosomes containing live Salmonella with the early endosome and repress their transport to lysosomes (61), whereas the acap1 gene promotes Salmonella invasion (62).
To explore the possible polarization of the Mm-infected zebrafish macrophages, we compared the differentially expressed gene set with the transcriptomes of M1 and M2 in vitro polarized macrophages reported by Beyer et al. (40). We found that 18 M1-enriched genes (log2FC > 1) were overexpressed in infected macrophages and 2 were downregulated whereas 26 M2-enriched genes (log2FC < −1) were upregulated and 6 were downregulated (Figure 5C). GSEA showed no significant association of either the upregulated or downregulated genes with either M1- or M2-enriched genes (Figure 5D, FDR > 0.05).
One of the most highly induced gene in infected macrophages was cxcl11aa, a zebrafish homolog of the gene for human CXCL11 (Figure 5A), a proinflammatory chemokine that is a typical M1 marker (63). We recently showed that this chemokine is important during Mm infection in zebrafish for the recruitment of macrophages and dissemination of the bacteria (39). Furthermore, expression of tnfa appeared to be highly upregulated in infected macrophages. Tnfa is one of the main markers of M1 activated macrophages in human and has been used as a marker for M1-like activated macrophages in zebrafish larvae (16). Other known zebrafish M1-like activated macrophage markers are non-significantly overexpressed (il1b, tnfb), or barely expressed (il6). On the other hand, the known zebrafish M2-like markers are either expressed at a low level (tgfb1a, il10) or not significantly downregulated (cxcr4b, ccr2) (Figures 5A,E).
We can conclude that the strong induction of two important proinflammatory markers, cxcl11aa and tnfa, and the downregulation of genes associated to M2 polarization as detected by GSEA indicate that Mm-infected macrophages display M1 rather than M2 characteristics.
The cxcl11aa Gene Expression as a Robust Marker of Mycobacterium-Infected Macrophages
Among the chemokine and cytokine genes expressed in Mm-infected macrophages, cxcl11aa emerged as the most reproducible infection marker from the RNAseq analysis, showing significantly higher induction (average log2 (fold change) = 8.6, P-adj < 0.001) than tnfa (average log2 (FC) = 5.6, P-adj = 0.03) in all replicates. In order to confirm the Mm-inducible expression of cxcl11aa in macrophages, we FACS-sorted mpeg1:mCherry positive cells from Mm-infected and mock-injected larvae and quantified the level of cxcl11aa expression by real time PCR. In uninfected conditions, the expression of cxcl11aa was significantly enriched in the mCherry-positive macrophage cell fraction compared with the unlabeled cell fraction (Figure 6A). During infection, the expression levels of cxcl11aa were strongly upregulated in macrophages but not in the unlabeled cell fraction. We found that the level of this infection-induced and macrophage-specific expression of cxcl11aa is high enough to be detectable in total RNA samples from whole larvae and that cxcl11aa induction did not require the bacterial locus RD1 (Region of Difference 1), a pathogenicity locus encompassing the secretion system of ESAT-6 (Early Secreted Antigenic Target 6 kDa), which is associated with mycobacterial virulence and formation of tubercular granulomas (Figure 6B) (65). The induction of cxcl11aa was also independent from the host cxcr3.2 gene, which encodes the receptor for Cxcl11aa (Figure 6B) (39). Next, we asked whether cxcl11aa induction requires the central immune mediator Myd88, which links pathogen recognition by Toll-like receptors and Il1β-mediated inflammation to activation of the transcription factor Nfκb (66). Therefore, we quantified the expression levels of cxcl11aa in myd88 mutant larvae. Since myd88 mutants display an increased infection level when infected with the same initial infection load as wild type controls, we compensated this with a reduced inoculum to obtain a similar infection level at 4 dpi. Both with the reduced and the regular inoculum, myd88 mutants displayed a marked incapability to upregulate cxcl11aa (Figures 6C–E), indicating that Myd88-dependent signaling is key to upregulate macrophage expression of cxcl11aa during Mm infection.
Figure 6. The expression of cxcl11aa is upregulated in macrophages upon infection and requires an active Myd88-immune signaling. (A) Expression of cxcl11aa in FACS-sorted macrophages (MΦ, mpeg1:mCherry positive) and its infection-dependent induction (relative to negative/Mock fraction). (B) Mm- (or Mock-) injected larvae (>100 per replicate per condition) were dissociated at 5 dpi. Induction of cxcl11aa does not require the RD1 pathogenicity locus and mutants of the cognate receptor of cxcl11aa (cxcr3.2−/−) are still able to upregulate cxcl11aa at comparable levels to wt. (C–E) RNA was isolated from pools of >10 whole larvae collected at 4 dpi. Eight hundred CFU of RD1 mutant bacteria vs. 100 CFU of wildtype Mm were injected to reach a comparable infection level at 4 dpi. Dependency of cxcl11aa induction from myd88. Expression levels (C), representative burden analysis (D) and representative burden pictures (E) derive from larvae collected at 4 dpi. RNA was isolated from pools of >10 whole larvae. Each point in (D) represents 1 infected larva from a representative pool. Two hundred CFU of wildtype Mm were injected in myd88+/+ larvae vs. 100 and 200 CFU injected in myd88−/− larvae to reach a comparable infection level at 4 dpi. Quantification of total bacterial pixels was obtained using dedicated bacterial pixel count program (64). Scale bar in (E) 200 μm. Statistical significance was analyzed by one-way ANOVA with Sidak post-hoc correction on ln(n)-transformed relative induction folds (real time PCRs) or untransformed data (infection burden). Significance (P-value) is indicated with: ns, non-significant; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. Error bars: mean ± s.e.m.
Zebrafish larvae provide unique possibilities for real time visualization of macrophage responses during developmental and disease processes. However, it has remained unknown how the expression profile of larval macrophages compares to the profiles of human M1 and M2 in vitro polarized macrophage subsets, which are commonly considered as a reference for pro-inflammatory or anti-inflammatory activation states. Here we used RNAseq analysis of FACS-sorted cell fractions to determine the expression profile of macrophages isolated from mpeg1 reporter lines, which are widely used for imaging studies in zebrafish due to the highly specific labeling of the macrophage lineage. We demonstrate the unique signature of the mpeg1 reporter cells by comparison with the RNAseq profiles of neutrophils, marked by the mpx reporter, and progenitor lymphocytes, marked by the lck reporter. We detected expression of homologs of human M1 as well as M2 markers in the mpeg1 reporter cells, indicating that zebrafish larval macrophages have the potential to differentiate into both directions. Finally, to demonstrate polarization of macrophages under challenged conditions, we achieved an RNAseq analysis of low numbers of mpeg1-positive macrophages infected with a mycobacterial pathogen. The profiling of these infected macrophages revealed downregulation of M2 markers, while M1 markers were upregulated, with strongest induction of a homolog of the human M1 marker CXCL11.
Adult mpeg1 reporter fish have previously been used to determine the transcriptome of microglia, the brain-resident macrophage population (67). Other fluorescent reporter lines for different immune cell types from the myeloid and lymphoid lineages have recently been used to determine single-cell transcriptomes of cells sorted from hematopoietic organs (spleen and kidney marrow) of adult fish (68–70). Our study is the first to report on the transcriptome of larval macrophages. This dataset provides a useful new data mining resource that will facilitate genetic analyses of macrophage-specific genes in zebrafish larval models for development and disease. A dual-fluorescent reporter line with mpeg1-labeled macrophages and tnfa as a marker for M1 phenotype has been used to demonstrate that injury and infection can induce M1 polarization of macrophages in zebrafish larvae (16). While the tnfa reporter does not show detectable fluorescent gene expression in the absence of wounding or infection stimuli, we could detect a basal level of tnfa expression in our RNAseq data of macrophages from unchallenged zebrafish larvae. Furthermore, our RNAseq data set of 437 enriched macrophage markers contains il1b, a M1 marker that was reported to be induced by injury in tnfa-positive macrophages, but also tgfb1b, a M2 marker that is expressed at higher levels in tnfa-negative macrophages. Single cell sequencing would be required to determine if all macrophages express these M1 and M2 markers at low levels or that distinct macrophage subsets exist already under unchallenged conditions. A comparison with RNAseq profiles of in vitro differentiated human M1 and M2 macrophages provided further evidence that the transcriptome of zebrafish larval macrophages displays a mixed M1 and M2 signature (40). Whereas, our results do not allow to conclude if two distinct populations of macrophages similar to human M1 and M2 polarized macrophages exist in zebrafish larvae, we identified several specific genes that suggest the presence of these different populations and that could be used to expand the repertoire of zebrafish transgenic reporter lines for investigating macrophage polarization in vivo during immune challenge in the zebrafish model.
Fluorescent reporters driven by the mpeg1 and mpx promoters distinguish specifically between macrophages and neutrophils (19, 27). In agreement, we did not detect mpeg1 gene expression in neutrophils from mpx reporter fish. However, we detected low levels of expression of mpx and other common neutrophil markers in macrophages from mpeg1 reporter fish. This indicates that the RNAseq procedure is highly sensitive and suggests that post-translational mechanisms contribute to regulating the specificity of innate immune cell types. We found that approximately half of the genes that show enriched expression in neutrophils also show enriched expression in macrophages. However, an obvious difference between the two innate immune cell types is that genes involved in antigen presentation and processing were detected only in macrophages. Other notable differences were found within the families of proteinases. The neutrophil RNAseq data reported here have been data mined to investigate the expression of the major classes of drug transporters in zebrafish larvae, providing useful information for optimizing screening approaches for anti-inflammatory drugs (71).
The enriched gene sets of larval macrophages and neutrophils consist for more than 80% of transcripts that are not detected in progenitor lymphocytes isolated from lck reporter fish. Similarly, the enriched gene set of lck-labeled lymphocytes consists for 92% of transcripts that are not expressed in the myeloid lineage. It is well-known that the adaptive immune system in zebrafish larvae is not yet mature and that full immunocompetence, including antibody production, is achieved only by 3–6 weeks of development (54, 55). However, it is not known at which stage of larval development the first interactions between antigen-presenting cells and T-lymphocytes take place. We found that mpeg1-labeled macrophages from 5 day old larvae express the major histocompatibility class II gene mhc2dab, which is even earlier than the observed expression of a fluorescent mhc2dab reporter that labels putative dendritic cells scattered throughout the skin of larvae from 9 days onwards (72). The presence of mhc2dab and transcripts of other genes involved in antigen presentation and processing in larval macrophages suggest that communication with T-lymphocytes could take place already at stages where zebrafish larvae are generally believed to rely exclusively on innate immunity. In support of this hypothesis, we found that larval lymphocytes express the Cd4 marker for helper T-cells and the co-stimulatory receptor Cd28, which are required for macrophage activation. Furthermore, the expression of a perforin gene (prf1.7) is indicative of the development of cytotoxic T-cells. However, since there was no detectable expression of Cd8, it is unlikely that cytotoxic T-cells are already functional in 5 day old larvae.
To investigate how larval macrophages respond to an intracellular infection with mycobacteria, we determined the expression profile of Mm-infected mpeg1 reporter cells. This RNAseq analysis was challenging due to the low numbers of infected cells that could be obtained by FACS sorting. Infected macrophages have a lifespan of <5 h (73), which likely is an important contributing factor to the difficulty of isolating Mm-infected cells. While different types of macrophage polarization have been reported for in vitro cultured macrophages infected with mycobacteria (15), it is not understood how these pathogens affect macrophage polarization during different stages of tuberculosis disease in vivo. Our RNAseq analysis was performed at a stage where infected macrophages have aggregated into inflammatory infection foci, which are regarded as the earliest developmental stages of tuberculous granulomas (74). We observed that homologs of M2-enriched transcripts of human cells were preferentially down-regulated in M. marinum-infected zebrafish macrophages, whereas several M1-enriched transcripts were highly upregulated. Therefore, although no clear polarization was observed, our analysis suggests that macrophages shift toward M1 phenotype in Mm-infected zebrafish, which are used to model tuberculosis. Our results show an important modification of the macrophage transcriptome upon mycobacterial infection and unravel several targets that can be studied to better understand the molecular mechanisms involved in the host-pathogen interaction.
An important question is whether part of the observed expression changes in Mm-infected macrophages might be triggered by bacterial virulence factors or that all changes represent a general host defense response that is mounted against pathogenic as well as non-pathogenic mycobacteria. Irrespective of the answer to this question, it can be argued that some of the induced genes benefit the pathogen rather than the host. For example, we detected induced expression in Mm-infected macrophages of genes (nsfb, acap1) that promote the survival of bacteria in phagosomes (61, 62). A gene that is induced strongly and reproducibly among all replicates, cxcl11aa, could have both host-beneficial and host-detrimental effects. This gene, which is a homolog of the human M1 marker CXCL11, encodes a chemokine that mediates macrophage recruitment to infection foci through interaction with chemokine receptor Cxcr3.2, the zebrafish counterpart of human CXCR3 (39). While a certain level of macrophage recruitment during Mm infection is necessary to restrict infection (75), Mm bacteria also take advantage of the arrival of new macrophages at infection foci as this promotes spreading of the infection (37). In line with these considerations, we have previously found that deficiency in the receptor for Cxcl11aa, Cxcr3.2, limits the expansion of Mm in granulomas (39). A similar phenotype has been found upon depletion of Mmp9, another host factor required for macrophage recruitment (38). Therefore, high and sustained induction of cxcl11aa is likely to have an adverse effect on the control of Mm infection by the zebrafish host. On the other hand, the robust induction of this M1 polarization marker makes the cxcl11aa gene a prime candidate to expand the collection of zebrafish reporter lines for studying macrophage activation.
In conclusion, the transcriptome analyses reported here present a unique and detailed genetic profile of zebrafish larval immune cells, thereby providing a valuable resource that can be data mined to verify the expression of specific genes in the profiled cell types or to identify novel genes of interest and potential cell-specific markers. In future work single cell RNA sequencing technology will be useful to interrogate the heterogeneity in expression profiles of resting and activated macrophages.
Data Availability Statement
The sequencing data for infected samples have been submitted to the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE68920. The sequencing data for uninfected samples were made previously available under accession number GSE78954. The sequencing data for human macrophages are available under the accession number GSE36952.
This study was carried out in accordance with the local animal welfare regulations as overseen by the Animal Welfare Body of Leiden University (License number: 10612). The protocol was approved by the international guidelines specified by the EU Animal Protection Directive 2010/63/EU.
JR contributed to conception and design of the study, performed and analyzed experiments, and wrote and edited the manuscript. VT contributed to conception and design of the study, performed and analyzed experiments, and wrote and edited a section of the manuscript. AZ and ZK contributed to conception and design of the study and performed experiments. HJ performed and analyzed experiments. FS performed experiments. AM contributed to conception and design of the study, wrote and edited the manuscript, and acquired funding. All authors contributed to manuscript revision, read and approved the submitted version.
Conflict of Interest Statement
HJ was employed by company ZF-screens B.V. (Leiden, Netherlands).
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The authors thank Michiel van der Vaart for critical reading of the manuscript, Graham Lieschke (Monash University), Georges Lutfalla (University of Montpellier), Steve Renshaw (University of Sheffield), and David Langenau (Harvard Stem Cell Institute) for zebrafish reporter lines, Guido de Roo, and Sabrina Veld for support with FACS sorting, and the fish facility team members for zebrafish care. This work was supported by the Marie Curie Initial Training Network FishForPharma (PITNGA-2011-289209) and the project ZF-HEALTH (HEALTHF4-2010-242048) funded by the European Commission 7th Framework Programme, and by the SmartMix programme of the Netherlands Ministry of Economic Affairs and the Ministry of Education, Culture, and Science. Additionally, ZK was supported by the Higher Education Commission of Pakistan, FS was supported by a CONACYT fellowship from the Mexican Government, and AZ was supported by a Horizon grant of the Netherlands Genomics Initiative.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2019.00832/full#supplementary-material
Figure S1. Correlation between RNA-sequencing samples. (A–C) HeatMap of Pearson correlation coefficient of 5 and 6 dpf mpeg1:gal4; UAS-Kaede, and mpeg1:mCherry positive and negative samples (A), 5dpf mpx:gfp and lck:gfp positive and negative samples (B) and 20 cell samples of mpeg1:gfp positive cells infected or not with M. marinum GFP (C).
Figure S2. Principal Component Analysis of RNA-sequencing samples. (A–C) Principal Component Analysis of 5 and 6 dpf mpeg1:gal4; UAS-Kaede, and mpeg1:mCherry positive and negative samples (A), 5dpf mpx:gfp and lck:gfp positive and negative samples (B), and 20 cell samples of mpeg1:gfp positive cells infected or not with M. marinum GFP (C).
Figure S3. RT-qPCR validation of macrophage and neutrophils markers by RT-qPCR. (A) Expression of chemokine receptors and lysosomal genes in FACS-sorted macrophages (MΦ, mpeg1:mCherry positive) and compared with non-fluorescent cells (neg). All genes tested are present in the zf core dataset. RT-qPCR were performed in two independent biological replicates. Y-axis is log10-transformed for ease of visualization. (B) Expression of macrophage and neutrophil marker genes in FACS-sorted macrophages (MΦ, mpeg1:mCherry positive) and neutrophils (NΦ, mpx:GFP positive), and compared with their respective non-fluorescent cells (neg). The genes mpeg1, mpeg1.2, marco, spi1a, grn1/2, and grna were found enriched in the macrophage transcriptome dataset, mpx and mmp9 were found in both macrophage and neutrophil datasets, whereas apoa4 and grnb were found in none of them. RT-qPCR results are similar to RNA-seq results, however neutrophils qPCR results provide more significantly enriched genes. RT-qPCR were performed in three independent biological replicates. Statistical significance was analyzed by one-way ANOVA with Bonferroni post-hoc correction on ΔCt before transformation. Significance (P-value) is indicated with: ns, non-significant; *P < 0.05; **P < 0.01; ***P < 0.001. Error bars: mean ± s.e.m.
Figure S4. Molecular function associated to the zebrafish core macrophage expression dataset. Network visualization of GO analysis enrichment (molecular function category) of genes from the core macrophage expression data set using BiNGO and EnrichmentMap. Node size corresponds to the number of genes associated to the enriched GO term and edge size to the similarity coefficient between two nodes.
Table S1. Additional gene annotation. List of manually annotated genes added to the Ensembl annotation version 79 from the genome version Zv9.
Table S2. Statistics on RNAseq samples. Top: TPM table of macrophage, neutrophil, and lymphoid cell population. TPM for each positive and negative samples were computed by using the longest transcript for each gene. Presented results are average between each replicate.
Table S3. Summary of gene expression and differential expression analysis. Level of gene expression was distributed based on average TPM values among non-expressed (TPM < 3), moderately expressed (3 ≤ TPM < 100) and highly expressed (TPM ≥ 100) for the positive samples and an average of all the negative samples present in this analysis. Fold Change (FC) between positive and negative samples were computed using DESeq as described in the material and methods.
Table S4. Zebrafish core macrophage expression dataset. The 437 genes enriched (log2FC ≥ 1 and P-adj < 0.01) in at least two of the three mpeg1 positive cell populations compared with their respective negative cell background represent the zebrafish core macrophage expression data set. Associated gene names and descriptions in bold are additional manual annotations to the Ensembl annotation. For each gene, one relevant term from the Biological Process GO is presented.
Table S5. Gene expression dataset in neutrophils. List of genes enriched (log2FC ≥ 1 and 0 and P-adj < 0.01) in mpx:gfp positive neutrophils compared with fluorescent negative cells. Associated name and description in bold are additional annotations to the Ensembl annotation. For each gene, one relevant term from the Biological Process GO is presented.
Table S6. Gene expression dataset in lymphoid cells. List of genes enriched (log2FC ≥ 1 and P-adj < 0.01) in lck:gfp positive lymphoid cells compared with fluorescent negative cells. Associated name and description in bold are additional annotations to the Ensembl annotation. For each gene, one relevant term from the Biological Process GO is presented.
Table S7. List of the human homologs from zebrafish genes. List of the human homologs from the genes enriched in the zebrafish core macrophage expression dataset.
Table S8. Differential expression analysis between Mycobacterium marinum infected and uninfected zebrafish macrophages. List of genes upregulated (FC > 0 and P-adj < 0.05) and downregulated (FC < 0 and P-adj < 0.05) in mpeg1:mCherry and Mm-GFP double positive macrophages compared with mpeg1:mCherry only positive macrophages. Associated name and description in bold are additional annotations to the Ensembl annotation.
5. Leibovich SJ, Polverini PJ, Shepard HM, Wiseman DM, Shively V, Nuseir N. Macrophage-induced angiogenesis is mediated by tumour necrosis factor-alpha. Nature. (1987) 329:630–2. doi: 10.1038/329630a0
6. Fantin A, Vieira JM, Gestri G, Denti L, Schwarz Q, Prykhozhij S, et al. Tissue macrophages act as cellular chaperones for vascular anastomosis downstream of VEGF-mediated endothelial tip cell induction. Blood. (2010) 116:829–40. doi: 10.1182/blood-2009-12-257832
7. Travnickova J, Tran Chau V, Julien E, Mateos-Langerak J, Gonzalez C, Lelievre E, et al. Primitive macrophages control HSPC mobilization and definitive haematopoiesis. Nat Commun. (2015) 6:6227. doi: 10.1038/ncomms7227
9. Murray PJ, Allen JE, Biswas SK, Fisher EA, Gilroy DW, Goerdt S, et al. Macrophage activation and polarization: nomenclature and experimental guidelines. Immunity. (2014) 41:14–20. doi: 10.1016/j.immuni.2014.06.008
13. Lugo-Villarino G, Verollet C, Maridonneau-Parini I, Neyrolles O. Macrophage polarization: convergence point targeted by Mycobacterium tuberculosis and HIV. Front Immunol. (2011) 2:43. doi: 10.3389/fimmu.2011.00043
14. Richardson ET, Shukla S, Sweet DR, Wearsch PA, Tsichlis PN, Boom WH, et al. Toll-like receptor 2-dependent extracellular signal-regulated kinase signaling in Mycobacterium tuberculosis-infected macrophages drives anti-inflammatory responses and inhibits Th1 polarization of responding T cells. Infect Immun. (2015) 83:2242–54. doi: 10.1128/IAI.00135-15
15. Trivedi NH, Yu JJ, Hung CY, Doelger RP, Navara CS, Armitige LY, et al. Microbial co-infection alters macrophage polarization, phagosomal escape, and microbial killing. Innate Immun. (2018) 24:152–62. doi: 10.1177/1753425918760180
16. Nguyen-Chi M, Laplace-Builhe B, Travnickova J, Luz-Crawford P, Tejedor G, Phan QT, et al. Identification of polarized macrophage subsets in zebrafish. Elife. (2015) 4:e07288. doi: 10.7554/eLife.07288
17. Nguyen-Chi M, Laplace-Builhe B, Travnickova J, Luz-Crawford P, Tejedor G, Lutfalla G, et al. TNF signaling and macrophages govern fin regeneration in zebrafish larvae. Cell Death Dis. (2017) 8:e2979. doi: 10.1038/cddis.2017.374
19. Ellett F, Pase L, Hayman JW, Andrianopoulos A, Lieschke GJ. mpeg1 promoter transgenes direct macrophage-lineage expression in zebrafish. Blood. (2011) 117:e49–56. doi: 10.1182/blood-2010-10-314120
20. Bernut A, Herrmann JL, Kissa K, Dubremetz JF, Gaillard JL, Lutfalla G, et al. Mycobacterium abscessus cording prevents phagocytosis and promotes abscess formation. Proc Natl Acad Sci USA. (2014) 111:E943–952. doi: 10.1073/pnas.1321390111
21. Benard EL, Racz PI, Rougeot J, Nezhinsky AE, Verbeek FJ, Spaink HP, et al. Macrophage-expressed perforins mpeg1 and mpeg1.2 have an anti-bacterial function in zebrafish. J. Innate Immun. (2015) 7:136–52. doi: 10.1159/000366103
22. Gerri C, Marin-Juez R, Marass M, Marks A, Maischein HM, Stainier DYR. Hif-1alpha regulates macrophage-endothelial interactions during blood vessel development in zebrafish. Nat Commun. (2017) 8:15492. doi: 10.1038/ncomms15492
25. Torraca V, Masud S, Spaink HP, Meijer AH. Macrophage-pathogen interactions in infectious diseases: new therapeutic insights from the zebrafish host model. Dis Model Mech. (2014) 7:785–97. doi: 10.1242/dmm.015594
26. Svahn AJ, Graeber MB, Ellett F, Lieschke GJ, Rinkwitz S, Bennett MR, et al. Development of ramified microglia from early macrophages in the zebrafish optic tectum. Dev Neurobiol. (2013) 73:60–71. doi: 10.1002/dneu.22039
28. Langenau DM, Ferrando AA, Traver D, Kutok JL, Hezel JP, Kanki JP, et al. In vivo tracking of T cell development, ablation, and engraftment in transgenic zebrafish. Proc Natl Acad Sci USA. (2004) 101:7369–74. doi: 10.1073/pnas.0402248101
30. Prajsnar TK, Cunliffe VT, Foster SJ, Renshaw SA. A novel vertebrate model of Staphylococcus aureus infection reveals phagocyte-dependent resistance of zebrafish to non-host specialized pathogens. Cell Microbiol. (2008) 10:2312–25. doi: 10.1111/j.1462-5822.2008.01213.x
31. Samrah Masud TKP, Gerda EML, Marianne B, Vincenzo T, Van Der Vaart M, Annemarie HM. Macrophages target Salmonella by Lc3-associated phagocytosis in a systemic infection model. Autophagy. (2018) 15:796–812. doi: 10.1080/15548627.2019.1569297
32. Mesureur J, Feliciano JR, Wagner N, Gomes MC, Zhang L, Blanco-Gonzalez M, et al. Macrophages, but not neutrophils, are critical for proliferation of Burkholderia cenocepacia and ensuing host-damaging inflammation. PLoS Pathog. (2017) 13:e1006437. doi: 10.1371/journal.ppat.1006437
36. Clay H, Davis JM, Beery D, Huttenlocher A, Lyons SE, Ramakrishnan L. Dichotomous role of the macrophage in early Mycobacterium marinum infection of the zebrafish. Cell Host Microbe. (2007) 2:29–39. doi: 10.1016/j.chom.2007.06.004
38. Volkman HE, Pozos TC, Zheng J, Davis JM, Rawls JF, Ramakrishnan L. Tuberculous granuloma induction via interaction of a bacterial secreted protein with host epithelium. Science. (2010) 327:466–9. doi: 10.1126/science.1179663
39. Torraca V, Cui C, Boland R, Bebelman JP, Van Der Sar AM, Smit MJ, et al. The CXCR3-CXCL11 signaling axis mediates macrophage recruitment and dissemination of mycobacterial infection. Dis Model Mech. (2015) 8:253–69. doi: 10.1242/dmm.017756
42. Volkman HE, Clay H, Beery D, Chang JC, Sherman DR, Ramakrishnan L. Tuberculous granuloma formation is enhanced by a mycobacterium virulence determinant. PLoS Biol. (2004) 2:e367. doi: 10.1371/journal.pbio.0020367
43. Van Der Sar AM, Abdallah AM, Sparrius M, Reinders E, Vandenbroucke-Grauls CM, Bitter W. Mycobacterium marinum strains can be divided into two distinct types based on genetic diversity and virulence. Infect Immun. (2004) 72:6306–12. doi: 10.1128/IAI.72.11.6306-6312.2004
44. Cui C, Benard EL, Kanwal Z, Stockhammer OW, Van Der Vaart M, Zakrzewska A, et al. Infectious disease modeling and innate immune function in zebrafish embryos. Methods Cell Biol. (2011) 105:273–308. doi: 10.1016/B978-0-12-381320-6.00012-6
45. Rougeot J, Zakrzewska A, Kanwal Z, Jansen HJ, Spaink HP, Meijer AH. RNA sequencing of FACS-sorted immune cell populations from zebrafish infection models to identify cell specific responses to intracellular pathogens. Methods Mol Biol. (2014) 1197:261–74. doi: 10.1007/978-1-4939-1261-2_15
46. Zakrzewska A, Cui C, Stockhammer OW, Benard EL, Spaink HP, Meijer AH. Macrophage-specific gene functions in Spi1-directed innate immunity. Blood. (2010) 116:e1–11. doi: 10.1182/blood-2010-01-262873
47. Kanwal Z, Zakrzewska A, Den Hertog J, Spaink HP, Schaaf MJ, Meijer AH. Deficiency in hematopoietic phosphatase ptpn6/Shp1 hyperactivates the innate immune system and impairs control of bacterial infections in zebrafish embryos. J Immunol. (2013) 190:1631–45. doi: 10.4049/jimmunol.1200551
48. Ordas A, Kanwal Z, Lindenberg V, Rougeot J, Mink M, Spaink HP, et al. MicroRNA-146 function in the innate immune transcriptome response of zebrafish embryos to Salmonella typhimurium infection. BMC Genomics. (2013) 14:696. doi: 10.1186/1471-2164-14-696
49. Stockhammer OW, Zakrzewska A, Hegedus Z, Spaink HP, Meijer AH. Transcriptome profiling and functional analyses of the zebrafish embryonic innate immune response to Salmonella infection. J Immunol. (2009) 182:5641–53. doi: 10.4049/jimmunol.0900082
50. Solchenberger B, Russell C, Kremmer E, Haass C, Schmid B. Granulin knock out zebrafish lack frontotemporal lobar degeneration and neuronal ceroid lipofuscinosis pathology. PLoS ONE. (2015) 10:e0118956. doi: 10.1371/journal.pone.0118956
51. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. (2015) 43:D447–452. doi: 10.1093/nar/gku1003
52. Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. (2003) 34:267–73. doi: 10.1038/ng1180
53. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. (2005) 102:15545–50. doi: 10.1073/pnas.0506580102
54. Lam SH, Chua HL, Gong Z, Lam TJ, Sin YM. Development and maturation of the immune system in zebrafish, Danio rerio: a gene expression profiling, in situ hybridization and immunological study. Dev Comp Immunol. (2004) 28:9–28. doi: 10.1016/S0145-305X(03)00103-4
55. Page DM, Wittamer V, Bertrand JY, Lewis KL, Pratt DN, Delgado N, et al. An evolutionarily conserved program of B-cell development and activation in zebrafish. Blood. (2013) 122:e1–11. doi: 10.1182/blood-2012-12-471029
56. Tian Y, Xu J, Feng S, He S, Zhao S, Zhu L, et al. The first wave of T lymphopoiesis in zebrafish arises from aorta endothelium independent of hematopoietic stem cells. J Exp Med. (2017) 214:3347–60. doi: 10.1084/jem.20170488
58. Qian BZ, Zhang H, Li J, He T, Yeo EJ, Soong DY, et al. FLT1 signaling in metastasis-associated macrophages activates an inflammatory signature that promotes breast cancer metastasis. J Exp Med. (2015) 212:1433–48. doi: 10.1084/jem.20141555
59. Dave RK, Dinger ME, Andrew M, Askarian-Amiri M, Hume DA, Kellie S. Regulated expression of PTPRJ/CD148 and an antisense long noncoding RNA in macrophages by proinflammatory stimuli. PLoS ONE. (2013) 8:e68306. doi: 10.1371/journal.pone.0068306
60. Repique CJ, Li A, Brickey WJ, Ting JP, Collins FM, Morris SL. Susceptibility of mice deficient in the MHC class II transactivator to infection with Mycobacterium tuberculosis. Scand J Immunol. (2003) 58:15–22. doi: 10.1046/j.1365-3083.2003.01266.x
61. Mukherjee K, Siddiqi SA, Hashim S, Raje M, Basu SK, Mukhopadhyay A. Live Salmonella recruits N-ethylmaleimide-sensitive fusion protein on phagosomal membrane and promotes fusion with early endosome. J Cell Biol. (2000) 148:741–53. doi: 10.1083/jcb.148.4.741
62. Davidson AC, Humphreys D, Brooks AB, Hume PJ, Koronakis V. The Arf GTPase-activating protein family is exploited by Salmonella enterica serovar Typhimurium to invade nonphagocytic host cells. MBio. (2015) 6:e02253–14. doi: 10.1128/mBio.02253-14
64. Stoop EJ, Schipper T, Rosendahl Huber SK, Nezhinsky AE, Verbeek FJ, Gurcha SS, et al. Zebrafish embryo screen for mycobacterial genes involved in the initiation of granuloma formation reveals a newly identified ESX-1 component. Dis Model Mech. (2011) 4:526–36. doi: 10.1242/dmm.006676
65. Majlessi L, Brodin P, Brosch R, Rojas MJ, Khun H, Huerre M, et al. Influence of ESAT-6 secretion system 1 (RD1) of Mycobacterium tuberculosis on the interaction between mycobacteria and the host immune system. J Immunol. (2005) 174:3570–9. doi: 10.4049/jimmunol.174.6.3570
66. Van Der Vaart M, Van Soest JJ, Spaink HP, Meijer AH. Functional analysis of a zebrafish myd88 mutant identifies key transcriptional components of the innate immune system. Dis Model Mech. (2013) 6:841–54. doi: 10.1242/dmm.010843
67. Oosterhof N, Holtman IR, Kuil LE, Van Der Linde HC, Boddeke EW, Eggen BJ, et al. Identification of a conserved and acute neurodegeneration-specific microglial transcriptome in the zebrafish. Glia. (2017) 65:138–49. doi: 10.1002/glia.23083
68. Athanasiadis EI, Botthof JG, Andres H, Ferreira L, Lio P, Cvejic A. Single-cell RNA-sequencing uncovers transcriptional states and fate decisions in haematopoiesis. Nat Commun. (2017) 8:2045. doi: 10.1038/s41467-017-02305-6
69. Carmona SJ, Teichmann SA, Ferreira L, Macaulay IC, Stubbington MJ, Cvejic A, et al. Single-cell transcriptome analysis of fish immune cells provides insight into the evolution of vertebrate immune cell types. Genome Res. (2017) 27:451–61. doi: 10.1101/gr.207704.116
70. Tang Q, Iyer S, Lobbardi R, Moore JC, Chen H, Lareau C, et al. Dissecting hematopoietic and renal cell heterogeneity in adult zebrafish at single-cell resolution using RNA sequencing. J Exp Med. (2017) 214:2875–87. doi: 10.1084/jem.20170976
71. Foulkes MJ, Henry KM, Rougeot J, Hooper-Greenhill E, Loynes CA, Jeffrey P, et al. Expression and regulation of drug transporters in vertebrate neutrophils. Sci Rep. (2017) 7:4967. doi: 10.1038/s41598-017-04785-4
73. Hosseini R, Lamers GE, Soltani HM, Meijer AH, Spaink HP, Schaaf MJ. Efferocytosis and extrusion of leukocytes determine the progression of early mycobacterial pathogenesis. J Cell Sci. (2016) 129:3385–95. doi: 10.1242/jcs.135194
74. Davis JM, Clay H, Lewis JL, Ghori N, Herbomel P, Ramakrishnan L. Real-time visualization of mycobacterium-macrophage interactions leading to initiation of granuloma formation in zebrafish embryos. Immunity. (2002) 17:693–702. doi: 10.1016/S1074-7613(02)00475-2
75. Pagan AJ, Yang CT, Cameron J, Swaim LE, Ellett F, Liesc0hke GJ, et al. Myeloid growth factors promote resistance to mycobacterial infection by curtailing granuloma necrosis through macrophage replenishment. Cell Host Microbe. (2015) 18:15–26. doi: 10.1016/j.chom.2015.06.008
Keywords: innate immunity, zebrafish, RNAseq, macrophage, mycobacteria, neutrophil, lymphoid progenitor cells
Citation: Rougeot J, Torraca V, Zakrzewska A, Kanwal Z, Jansen HJ, Sommer F, Spaink HP and Meijer AH (2019) RNAseq Profiling of Leukocyte Populations in Zebrafish Larvae Reveals a cxcl11 Chemokine Gene as a Marker of Macrophage Polarization During Mycobacterial Infection. Front. Immunol. 10:832. doi: 10.3389/fimmu.2019.00832
Received: 09 January 2019; Accepted: 29 March 2019;
Published: 17 April 2019.
Edited by:Tiehui Wang, University of Aberdeen, United Kingdom
Reviewed by:Magdalena Chadzińska, Jagiellonian University, Poland
Patricia Pereiro, Spanish National Research Council (CSIC), Spain
Marco Gerdol, University of Trieste, Italy
Copyright © 2019 Rougeot, Torraca, Zakrzewska, Kanwal, Jansen, Sommer, Spaink and Meijer. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Annemarie H. Meijer, firstname.lastname@example.org