Characterization of immune response against Mycobacterium marinum infection in the main hematopoietic organ of adult zebrafish (Danio rerio)

Contents lists available at ScienceDirect Developmental and Comparative Immunology journal homepage: www.elsevier.com/locate/devcompimm Characterization of immune response against Mycobacterium marinum infection in the main hematopoietic organ of adult zebrafish (Danio rerio) Sanna-Kaisa E. Harjula a , Anni K. Saralahti a , Markus J.T. Ojanen a,b , Tommi Rantapero c , Meri I.E. Uusi-Mäkelä a , Matti Nykter c , Olli Lohi d , Mataleena Parikka e,f , Mika Rämet a,g,h,i,∗ a Laboratory of Experimental Immunology, BioMediTech, Faculty of Medicine and Health Technology, FI-33014, Tampere University, Finland b Laboratory of Immunoregulation, BioMediTech, Faculty of Medicine and Health Technology, FI-33014, Tampere University, Finland c Laboratory of Computational Biology, BioMediTech, Faculty of Medicine and Health Technology, FI-33014, Tampere University, Finland d Tampere Center for Child Health Research, Tampere University and Tays Cancer Center, Tampere University Hospital, FI-33014, Tampere University, Finland e Laboratory of Infection Biology, BioMediTech, Faculty of Medicine and Health Technology, FI-33014, Tampere University, Finland f Oral and Maxillofacial Unit, Tampere University Hospital, P.O. Box 2000, FI-33521, Tampere, Finland g Department of Pediatrics, Tampere University Hospital, P.O. Box 2000, FI-33521, Tampere, Finland h PEDEGO Research Unit, Medical Research Center Oulu, P.O. Box 8000, FI-90014, University of Oulu, Finland i Department of Children and Adolescents, Oulu University Hospital, P.O. Box 10, FI-90029, OYS, Finland ARTICLE INFO Keywords: Zebrafish Mycobacterium marinum Transcriptome analysis Forward genetic screen ABSTRACT Tuberculosis remains a major global health challenge. To gain information about genes important for defense against tuberculosis, we used a well-established tuberculosis model; Mycobacterium marinum infection in adult zebrafish. To characterize the immunological response to mycobacterial infection at 14 days post infection, we performed a whole-genome level transcriptome analysis using cells from kidney, the main hematopoietic organ of adult zebrafish. Among the upregulated genes, those associated with immune signaling and regulation formed the largest category, whereas the largest group of downregulated genes had a metabolic role. We also performed a forward genetic screen in adult zebrafish and identified a fish line with severely impaired survival during chronic mycobacterial infection. Based on transcriptome analysis, these fish have decreased expression of several immunological genes. Taken together, these results give new information about the genes involved in the defense against mycobacterial infection in zebrafish. 1. Introduction In 2017, approximately 10 million people worldwide developed tuberculosis, a pulmonary or disseminated infection caused by Mycobacterium tuberculosis (World Health Organization, 2018). Although M. tuberculosis infection can be cleared by innate immunity (Verrall et al., 2014), it more typically leads to a latent phase. According to current estimations, 23% of the world's population are asymptomatic carriers of the bacteria (World Health Organization, 2018). The currently available Bacillus Calmette-Guérin (BCG) tuberculosis vaccine protects infants from disseminated tuberculosis but is less effective against pulmonary disease in adults. This vaccine is also unable to prevent the reactivation of latent disease (Mangtani et al., 2014;Tang et al., 2016). Thus, BCG vaccination cannot prevent the spread of tuberculosis. Based on genome-wide association studies (GWAS) and candidate gene studies in various human populations, as well as on animal studies, host genetics contributes to the susceptibility of developing active tuberculosis (Berg et al., 2016;Cooper et al., 1993,1997;Filipe-Santos et al., 2006;Harjula et al., 2018a;Sullivan et al., 2005;Tobin et al., 2010;Yim and Selvaraj, 2010). Human polymorphisms in genes encoding human leukocyte antigens (HLA) (Dallmann-Sauer et al., 2018; Yim and Selvaraj, 2010), Interferon (IFN) gamma, Interleukin10 and various chemokines and their receptors (Yim and Selvaraj, 2010), to mention but a few, have been associated with tuberculosis susceptibility. For some of these genes, there is also experimental evidence https://doi.org/10.1016/j.dci.2019.103523 Received 13 September 2019; Received in revised form 9 October 2019; Accepted 14 October 2019 Abbreviations: BCG, Bacille Calmette-Guérin; CFU, colony forming unit; dpf, days post fertilization; dpi, days post infection; GWAS, genome-wide association study; hpi, hours post infection; RIN, RNA integrity number; TL, Tüpfel long fin; wpi, weeks post infection; WT, wild type ∗ Corresponding author. Faculty of Medicine and Health Technology, FI-33014, Tampere University, Finland. E-mail addresses: [email protected] (S.-K.E. Harjula), [email protected] (A.K. Saralahti), [email protected] (M.J.T. Ojanen), [email protected] (T. Rantapero), [email protected] (M.I.E. Uusi-Mäkelä), [email protected] (M. Nykter), [email protected] (O. Lohi), [email protected] (M. Parikka), [email protected] (M. Rämet). Developmental and Comparative Immunology 103 (2020) 103523 Available online 15 October 2019 0145-305X/ © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/BY-NC-ND/4.0/). T from animal models on their importance in disease resistance (Beamer et al., 2008;Cooper et al., 1993;Cyktor et al., 2013;Flynn et al., 1993;Harjula et al., 2018a;Higgins et al., 2009;Peters et al., 2001;Torraca et al., 2015). During the last couple of decades, Mycobacterium marinum infection in both zebrafish embryos/larvae and adult zebrafish has become a widely used model of tuberculosis (Myllymäki et al., 2015,2016). M. marinum is a natural zebrafish pathogen and, similarly to its close relative M. tuberculosis (Stinear et al., 2008), it infects macrophages (Barker et al., 1997;El-Etr et al., 2001) and can cause a latent or slowly progressive infection, which can be reactivated by immunosuppression (Harjula et al., 2018a;Myllymäki et al., 2018;Parikka et al., 2012). Also comparable to tuberculosis, M. marinum infection both in zebrafish larvae and adult zebrafish results in granuloma formation (Davis et al., 2002;Myllymäki et al., 2018;Parikka et al., 2012;Swaim et al., 2006). In fact, the notion that granuloma not only benefits the host but also gives the bacteria a good environment in which to survive and spread the disease, originates from zebrafish studies (Davis and Ramakrishnan, 2009;Volkman et al., 2004). The host response to the M. marinum infection in zebrafish has been studied in several transcriptome analyses (Benard et al., 2016;Hegedüs et al., 2009;Kenyon et al., 2017;Meijer et al., 2005;Ojanen et al., 2019;Rotman et al., 2011;Rougeot et al., 2014,2019;van der Sar et al., 2009;van der Vaart et al., 2012;Veneman et al., 2015). These studies show that the gene expression profile of M. marinum infected zebrafish changes as the infection progresses. Consequently, in zebrafish larvae infected at 1 day post fertilization (dpf), three stages of the transcriptional response have been recognized: early, mid and late phase (Benard et al., 2016). Characteristic of the early phase is a transcriptional response, including upregulation of genes involved in response to bacterium, proteolysis and cellular macromolecular complex assembly and a downregulation of protein folding-associated genes at 2 h post infection (hpi), concurrent with phagocytosis of the bacteria (Benard et al., 2014). During mid-phase, starting at 6 hpi, the number of the differentially regulated genes is low. This is followed by an increase in differentially expressed genes during late phase (4 and 5 days post infection, dpi) (Benard et al., 2016). According to the aforementioned studies, the exact response not only depends on the stage of the infection but also on the virulence of the M. marinum strain (Rotman et al., 2011;van der Sar et al., 2009). In general for adult zebrafish, genes encoding for zinc-finger proteins, immune-related transcription factors and proteins related to apoptosis, among others, are upregulated at the early stages of the M. marinum infection (van der Sar et al., 2009). At later stages, genes encoding for immune-related transcription factors, cytokines, chemokine receptors, complement components, matrix metalloproteases and lysosomal proton transporters, among others, are induced independently of the virulence of the bacteria in the adult zebrafish (Meijer et al., 2005; van der Sar et al., 2009). In the present study, we performed a genome-wide transcriptome analysis using cells from kidney, the main hematopoietic organ in zebrafish, at 14 days post M. marinum infection. In this manner, we focused the analysis predominantly on immune cells. Furthermore, the immune response of wild type (WT) zebrafish was compared to the immune response of a susceptible mutant fish line, identified from an ongoing forward genetic screen. 2. Materials and methods 2.1. Zebrafish lines and maintenance 3 to 15 month-old WT AB and TL (Tüpfel long fin,gja5b t1/t1 ,lof dt2/ dt2 ) zebrafish (Danio rerio) lines from the Tampere Zebrafish Core Facility were used for the experiments. The rag hu1999/hu1999 mutant fish obtained from the Zebrafish International Resource Center (ZIRC, University of Oregon, Eugene, Oregon, USA) at the age of 6–8 months were used as a positive control in the infection experiments. The zebrafish were maintained according to standard protocols (Nüsslein-Volhard and Dahm, 2002). Briefly, the unchallenged zebrafish were kept in a standard flow-through system (Aquatic Habitats, Apopka, Florida, USA) with a light/dark cycle of 14/10 h and fed once a day with SDS 400 (Special Diet Services, Witham, Essex, UK) and once a day with in-house cultured Artemia nauplia. Alternatively, fish were fed once a day with GEMMA Micro 500 (Skretting, Stavanger, Norway). The M. marinum infected fish were kept in a standard flow-through unit (Aqua Schwarz GmbH, Göttingen, Germany) with the same light/dark cycle as the unchallenged fish and fed twice a day with SDS 400 or once a day with GEMMA Micro 500. Zebrafish embryos/larvae were maintained in embryonic medium (5 mM NaCL, 0.17 mM KCl, 0.33 mM CaCl2, 0.33 mM MgSO4, 10% methylene blue at 28.5 °C, fed with SDS100 (Special Diet Services) or GEMMA Micro 75 (Skretting) starting at 5 dpf and transferred to the flow-through system at 6 dpf. The well-being of the fish was monitored daily and the humane endpoint criteria determined in the animal experiment permits were used. The Animal Experiment Board has approved the housing and care of the zebrafish and all the conducted experiments (permits ESAVI/ 4234/04.10.03/2012, ESAVI/6403/04.10.03/2012, ESAVI/10079/ 04.10.06/2015, ESAVI/2776/2019, ESAVI/10539/2019, LSLH-20077254/Ym-23 and ESAVI/10366/04.10.07/2016). Furthermore, the Finnish Act on the Protection of Animals Used for Scientific or Educational Purposes (497/2013) as well as the EU Directive on the Protection of Animals Used for Scientific Purposes (2010/63/EU) were applied during this study. 2.2. The gene-breaking Tol2 transposon based mutagenesis in zebrafish embryos The gene-breaking transposon pGBT-RP2-1 (RP2) vector (Addgene plasmid # 31828; http://n2t.net/addgene:31828; RRID:Addgene_31828) and pT3TS-Tol2 vector (Addgene plasmid # 31831; http://n2t.net/addgene:31831; RRID:Addgene_31831) for the production of tol2 mRNA were received as a generous gift from Professor Stephen C. Ekker's laboratory (Mayo Clinic, Rochester, Minnesota, USA) (Balciunas et al., 2006;Clark et al., 2011). Both plasmids were transformed into E. coli One Shot TOP10 cells (Invitrogen™, Thermo Fisher Scientific, Waltham, Massachusetts, USA). Plasmid DNA was extracted with QIAGEN Plasmid Plus Maxi Kit (Qiagen, Hilden, Germany) and sequenced to confirm successful transformation. In order to produce tol2 mRNA, pT3TS-Tol2 plasmid was first linearized by digestion with FastDigest BamHI (Thermo Scientific, Thermo Fisher Scientific). Then, the linearized plasmid was used as a template for in vitro transcription performed with mMESSAGE mMACHINE Transcription Kit, (Invitrogen™, Thermo Fisher Scientific) according to the manufacturer's instructions. The zebrafish mutagenesis was conducted as previously described (Clark et al., 2011). Briefly, 12.5 pg of RP2 plasmid and tol2 mRNA in 1x PBS with 0.6% phenol red (Sigma-Aldrich, Saint Louis, Missouri, USA) were injected into the cell of 1-cell stage WT AB zebrafish embryos with a borosilicate capillary needle (Sutter Instrument Co., Novato, California, USA) using a PV830 Pneumatic PicoPump (World Precision Instruments, Sarasota, Florida, USA) and a micromanipulator (Narishige International, London, UK). If the mutagenesis was successful, RP2 was randomly inserted into the zebrafish genome causing one or several mutations per fish (Clark et al., 2011). Mutation carrying embryos (the founder fish in the genetic screen, designated as F0) were selected based on the expression of Green fluorescent protein (GFP) under a Lumar V.12 fluorescence stereomicroscope (Carl Zeiss MicroImaging GmbH, Göttingen, Germany) or Nikon AZ100 Fluorescence Macroscope (Nikon, Minato, Tokyo, Japan) (Clark et al., 2011). S.-K.E. Harjula, et al. Developmental and Comparative Immunology 103 (2020) 103523 2 2.3. Generation of the mutant zebrafish lines In order to produce zebrafish lines carrying unknown mutations, F0 founder fish were first crossed to TL. The fish in the resulting F1 generation carried various sets of mutations. These fish were named with a running number (also the final name of the subsequent resultant mutant zebrafish line) and crossed to TL. This resulted in an F2 generation consisting of fish heterozygous for the mutations. The F2 generation was then incrossed to produce an F3 generation with a theoretical Mendelian ratio of WT, heterozygous and homozygous fish for each mutation. Further fish generations were produced by incrossing the previous generation (until F5). In some cases, in order to maintain the line, the fish were again crossed to TL. For the M. marinum infection experiments on embryos, the F7 generation was utilized. Mutation carrying progeny in each generation were selected for based on GFP expression. 2.4. Experimental M. marinum infections The culture of M. marinum strain ATCC 927 and the inoculation into adult zebrafish were done as described previously (Harjula et al., 2018a;Parikka et al., 2012). For these infections, the zebrafish were anesthetized with 0.02% 3-amino benzoic acid ethyl ester (Sigma-Aldrich). Following this, 5 μl of M. marinum suspended in 1x PBS with 0.3 mg/ml phenol red (Sigma-Aldrich) was injected into the abdominal cavity of the fish using a 30 gauge Omnican 100 insulin needle (Braun, Melsungen, Germany). The survival of the fish was followed for 12–14 weeks. During the screen, 2–10 mutant lines were included in one experiment together with the WT line and occasionally rag hu1999/hu1999 mutants. For the lines showing increased or decreased susceptibility, the survival assay was performed at least three times in total. For the zebrafish embryo infections, M. marinum was suspended in 1x PBS with 0.3 mg/ml phenol red (Sigma-Aldrich). At 1 dpf, zebrafish embryos were anesthetized with 0.02% 3-amino benzoic acid ethyl ester (Sigma-Aldrich) and 2 nl of bacterial solution was injected into the caudal vein with a borosilicate capillary needle using a PV830 Pneumatic PicoPump (World Precision Instruments) and a micromanipulator (Narishige International). The survival of the embryos was monitored once a day for seven days. The infection dose (colony forming units, CFU) in both adult and embryo survival assays was verified by plating injected bacterial suspension on 7H10 agar (Becton Dickinson and Company, Franklin Lakes, New Jersey, USA) plates. 2.5. Experimental Streptococcus pneumoniae infections The culture of Streptococcus pneumoniae WT strain TIGR4 (T4), of serotype 4 and sequence type 205, and the inoculation into zebrafish embryos were done as described previously (Aaberge et al., 1995; Rounioja et al., 2012). Briefly, 5% lamb blood agar plates (TammerTutkan maljat Oy, Tampere, Finland) were used to grow T4 in 37 °C and 5% CO 2 overnight. From the plate, T4 was suspended in 5 ml of Todd Hewitt broth (Becton, Dickinson and Company) and grown from OD 620 0.1 to OD 620 of 0.4. After this, bacteria were suspended in 0.2 M KCl with 1% of 70 kDa Rhodamine Dextran (Invitrogen™, Thermo Fisher Scientific). At 2 dpf, zebrafish embryos were anesthetized with 0.02% 3-amino benzoic acid ethyl ester (Sigma-Aldrich) and 2 nl of the bacterial solution was microinjected into the blood circulation valley with a borosilicate capillary needle (Sutter Instrument Co.) using a PV830 Pneumatic PicoPump (World Precision Instruments) and a micromanipulator (Narishige International). The survival was followed once a day for five days. The infection dose (CFU) was verified by plating injected bacterial suspension on 5% lamb blood agar plates (TammerTutkan maljat Oy). 2.6. RNA extraction and quality control For the transcriptome analysis, the RNA was extracted from zebrafish kidney with the Qiagen RNeasy Mini Kit (Qiagen) according to the manufacturer's protocol. The removal of genomic DNA from the samples was done with the RapidOut DNA Removal Kit (Thermo Scientific, Thermo Fisher Scientific). The purity of the samples was checked with a NanoDrop™ 2000 Spectrophotometer (Thermo Scientific, Thermo Fisher Scientific) and the concentration was measured with a Qubit™ RNA BR Assay Kit (Invitrogen™, Thermo Fisher Scientific). The integrity of RNA was checked with a Fragment Analyzer (Advanced Analytical Technologies, Iowa, USA) using the Standard Sensitivity RNA Analysis Kit (Advanced Analytical Technologies) and the PROSize®2.0 Data Analysis Software (Advanced Analytical Technologies). The samples with an RNA integrity number (RIN) ≥ 7.8 were chosen for the RNA sequencing. 2.7. Whole genome transcription analysis by RNA sequencing The preparation of the cDNA library and the RNA sequencing were performed at Novogene, Hong Kong. The cDNA library of 250–300 bp was prepared with the 150 bp paired-end sequencing on the Illumina platform yielding > 20 million reads/sample. 2.8. RNA sequencing data analysis The quality of the reads was inspected using FastQC (Andrews, 2010). The reads were aligned against the GRCz10 reference genome with STAR using default parameters (Dobin et al., 2013). The expressions of genes were quantified with FeatureCounts (Liao et al., 2014) using Ensembl GRCz10.91 as the reference gene set (Hubbard et al., 2002). The normalization of the raw expressions and the differential gene expression analysis were both conducted using R-package DESeq2 (Love et al., 2014). To analyze the transcriptome data further, we included genes which were expressed ≥ |3|-fold between the groups and also had medians of ≥ |3|-fold. Gene ontology enrichment analysis was performed with The Gene Ontology Enrichment Analysis (GO Ontology Database released 1st of January 2019 and Panther Overrepresentation Test released 13th of November 2018) (Ashburner et al., 2000;Mi et al., 2017;The Gene Ontology Consortium, 2019) using the unranked list of upregulated or downregulated genes as the target list and the list of all the protein coding genes present in the RNA sequencing data as a background list. Classification of upregulated genes was based on the data available from Ensembl genome browser (Zerbino et al., 2018) versions 91 (version used for the data analysis) and 95, Reactome (Fabregat et al., 2018), The Zebrafish Information Network ZFIN (Howe et al., 2013), InterPro (Mitchell et al., 2019), The PROSITE database (Sigrist et al., 2013), Pfam 32.0 (El-Gebali et al., 2019), the NCBI Gene database (Gene, 2004), the NCBI BioSystems database (Geer et al., 2010) and the literature. Downregulated genes were classified according to the gene ontology data available in Ensembl, ZFIN and NCBI gene databases. If the gene did not have an official name/description, the preferred name of the gene from the NCBI database was used if available. 2.9. qPCR Reverse transcription of the RNA samples was done using the SensiFAST™ cDNA synthesis kit (BioLine, London, UK). cDNA was then used for determining the relative gene expression levels of the target genes with quantitative PCR (qPCR) using PowerUp™ SYBR®master mix (Applied Biosystems™, Thermo Fisher Scientific). Supplementary Table 1 gives the sequences of the qPCR primers used and the Ensembl gene identification codes for the analyzed genes. 2 −ΔCt method was used for calculating the expression levels of target genes relative to the expression of eef1a1l1 (Tang et al., 2007). M. marinum burden (CFU) S.-K.E. Harjula, et al. Developmental and Comparative Immunology 103 (2020) 103523 3 was quantified from the total DNA of the infected fish with qPCR using SensiFAST™ SYBR®No-ROX (Bioline) as previously described (Parikka et al., 2012), considering 100 CFU as the detection limit. DNA for the analysis was extracted from the abdominal organ blocks (excluding kidney) with TRI Reagent®(Molecular Research Center, Cincinnati, Ohio, USA) according to manufacturer's protocol. A CFX96 qPCR machine (Bio-Rad, California, USA) was used to perform qPCR and the BioRad CFX Manager software v3.1 (Bio-Rad) to analyze the data. Random RNA samples with no reverse transcription and non-template controls were used in qPCR to control for genomic DNA and other contamination. A melt curve analysis and 1.5% TAE agarose gel electrophoresis was used to validate the specificity of the qPCR products. PCR products below the detection limit or the ones with incorrect melt curves were given a Ct value of 40 when the gene expression analysis was performed. 2.10. Statistical analysis and power calculations Sample size calculations were conducted as described earlier (Harjula et al., 2018a). As for the M. marinum quantification, based on the high mortality of mutant463 zebrafish and our previous results (Myllymäki et al., 2018;Parikka et al., 2012), the difference between the groups was this time estimated to be 0.75 unit on the log10 scale resulting in the group size minimum of 7. In the gene ontology enrichment analyses Pvalues were calculated with Fisher's exact test and Bonferroni correction was used for multiple testing. Statistical analyses for other results were conducted with Prism, version 5.02 (GraphPad Software, Inc, California, USA). For the survival analyses a log-rank (Mantel-Cox) test was used. To test whether the bacterial count is increased in the mutant463 line compared to WT and whether the differences in gene expression detected by RNA sequencing are repeated when analyzed with qPCR, a nonparametric one-tailed Mann-Whitney test was used. Pvalues of < 0.05 were considered significant. 2.11. Data management The RNA-sequencing data discussed in this publication (Harjula et al., 2018b,2018c) have been deposited in NCBI's Gene Expression Omnibus (Edgar et al., 2002) and are accessible through GEO Series accession numbers GSE118288 (https://www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?acc=GSE118288) for the WT data and GSE118350 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE118350) for the mutant463 data. 3. Results 3.1. A low-dose infection with M. marinum leads to induction of immune response and changes in expression of genes regulating metabolism in kidney derived cells in adult zebrafish We have previously shown that a low-dose M. marinum infection into the abdominal cavity of adult zebrafish leads to a latent or slowly progressive infection, which can be used to model different phases of tuberculosis (Harjula et al., 2018a;Myllymäki et al., 2018;Parikka et al., 2012). In order to characterize the immune response to low-dose M. marinum infection in WT zebrafish, we infected 5–6 month-old zebrafish with 5–9 CFU of M. marinum and performed a whole-genome level transcriptome analysis from their kidneys. Since the kidney is the main hematopoietic organ in adult zebrafish (Davidson and Zon, 2004), it was selected to characterize the immunological response to the infection, particularly in immune cells. At 7, 14 and 28 days post infection (dpi), we collected kidneys for RNA extraction. From the fish collected at 14 dpi, we selected four males for each group for RNA sequencing. As a negative control, we used four 7 month-old unchallenged WT male fish. The yield from RNA sequencing was 23–28 million reads per sample and 90% of the reads were successfully mapped to the transcript database Ensembl GRCz10.91. In transcriptome analysis, we found 201 protein coding genes to be differentially expressed at 14 days post a low-dose M. marinum infection compared to the unchallenged WT fish (Tables 1–3,Supplementary Tables 2–3). More specifically, 96 of the differentially expressed protein coding genes were upregulated and 105 downregulated. In addition, we identified 21 upregulated and 30 downregulated noncoding RNAs (Supplementary Tables 4 and 5). We performed Gene Ontology Consortium enrichment analysis for both the upregulated and downregulated protein coding genes. According to this analysis, all the significantly enriched biological processes (P< 0.05), among the analyzed 86 upregulated genes recognized by the tool, were related to immune response (Table 1). This emphasizes the feasibility of using kidney tissue to characterize the immune response to infection. Among the most enriched processes were inflammatory response, response to bacterium and chemotaxis (Table 1). Next, we divided all the upregulated genes into subgroups based on the data available the databases listed in Chapter 2.8., as well as on the existing literature (Fig. 1A, Tables 2 and 3). The categories included immune response (43 genes), lipid binding and metabolism (9 genes), other metabolic role (13 genes), other role (27 genes) and unknown role (4 genes) (Fig. 1A, Tables 2 and 3). Downregulated protein coding genes are shown in Fig. 1B and Supplementary Tables 2–3. Among the analyzed 92 downregulated genes, many of the significantly (P< 0.05) enriched processes were Table 1 Enriched processes from gene ontology analysis of the upregulated protein coding genes. The Gene Ontology Consortium Enrichment analysis was performed using upregulated genes as a target list and all the protein coding genes in the RNA sequencing data as a background list. Pvalues were calculated with Fisher's exact test and Bonferroni correction was used for multiple testing. P< 0.05 was considered significant. GO term Description Pvalue GO:0006952 defense response 1.05E-06 GO:0006955 immune response 1.48E-06 GO:1990266 neutrophil migration 1.53E-06 GO:0097530 granulocyte migration 2.57E-06 GO:0097529 myeloid leukocyte migration 4.58E-06 GO:0070098 chemokine-mediated signaling pathway 7.44E-06 GO:0050900 leukocyte migration 7.85E-06 GO:1990869 cellular response to chemokine 9.50E-06 GO:1990868 response to chemokine 9.50E-06 GO:0030593 neutrophil chemotaxis 2.86E-05 GO:0006954 inflammatory response 4.02E-05 GO:0071621 granulocyte chemotaxis 4.67E-05 GO:0060326 cell chemotaxis 8.22E-05 GO:0030595 leukocyte chemotaxis 8.75E-05 GO:0002376 immune system process 1.44E-04 GO:0009617 response to bacterium 1.72E-04 GO:0006959 humoral immune response 4.57E-04 GO:0019221 cytokine-mediated signaling pathway 1.15E-03 GO:0071345 cellular response to cytokine stimulus 1.37E-03 GO:0043207 response to external biotic stimulus 2.06E-03 GO:0051707 response to other organism 2.06E-03 GO:0009607 response to biotic stimulus 2.13E-03 GO:0061844 antimicrobial humoral immune response mediated by antimicrobial peptide 2.42E-03 GO:0034097 response to cytokine 5.04E-03 GO:0006950 response to stress 6.42E-03 GO:0019730 antimicrobial humoral response 6.53E-03 GO:0002440 production of molecular mediator of immune response 1.07E-02 GO:0002377 immunoglobulin production 1.07E-02 GO:0051704 multi-organism process 1.10E-02 GO:0006935 chemotaxis 1.54E-02 GO:0042330 taxis 2.26E-02 GO:0009605 response to external stimulus 2.33E-02 GO:0071222 cellular response to lipopolysaccharide 3.91E-02 GO:0071219 cellular response to molecule of bacterial origin 3.91E-02 GO:0071216 cellular response to biotic stimulus 4.36E-02 S.-K.E. Harjula, et al. Developmental and Comparative Immunology 103 (2020) 103523 4 associated with metabolism, in particular with lipid metabolism (Supplementary Table 2). Since the gene ontology tool did not recognize all the downregulated protein coding genes, they were also manually classified into different biological processes, based on the gene ontology data available. The resultant downregulated subgroups were: immune system (4 genes), lipid metabolism and transport (10 genes), other metabolic process (the largest category with 36 genes), other localization process (8 genes), signaling and regulation (9 genes), other process (10 genes) and an unknown process (28 genes) (Fig. 1B, Supplementary Table 3). Of note, there was a clear expression of several blood cell and hemoglobin related genes in addition to known kidney expressed genes (Elmonem et al., 2018;Song et al., 2004;Walters et al., 2010) such as chemokine (CXC motif), receptor 4b (cxcr4b), aminoleculinate, delta-, Table 2 Immune response -related protein coding genes induced in mycobacterial infection in adult zebrafish kidney at 14 dpi (days post infection). The table shows the fold change in expression in M. marinum infected zebrafish compared to the unchallenged fish (n= 4 in both groups). Fold change represents the fold change between the group medians. The table includes only the genes with at least two samples with ≥ 20 normalized reads after infection and the genes whose expression was induced at least 3.0-fold (calculated both with the DEseq-tool and from the medians). Gene symbol Gene name/description Ensembl gene ID Fold change (median) Reference Acute phase response and antimicrobial activity saa serum amyloid A ENSDARG00000045999 21.8 hamp hepcidin antimicrobial peptide ENSDARG00000102175 20.5 lygl2 lysozyme g-like 2 ENSDARG00000099562 7.4 Mohapatra et al., 2019 leap2 liver-expressed antimicrobial peptide 2 ENSDARG00000104654 3.3 Neutrophil degranulation SERPINB8 (1 of many) zgc:173729 ENSDARG00000057263 8.7 plaub plasminogen activator, urokinase b ENSDARG00000039145 4.2 cst14b.1 cystatin 14b, tandem duplicate 1 ENSDARG00000045980 3.2 serpinb1l4 serpin peptidase inhibitor, clade B (ovalbumin), member 1, like 4 ENSDARG00000096888 3.1 krt5 keratin 5 ENSDARG00000058371 3.0 Inflammasome si:ch211-236p5.3 NACHT, LRR and PYD domains-containing protein 3-like ENSDARG00000086418 13.3 Hu et al., 2017 si:ch211-233m11.1 NACHT, LRR and PYD domains-containing protein 12 ENSDARG00000074653 6.8 fosab v-fos FBJ murine osteosarcoma viral oncogene homolog Ab ENSDARG00000031683 3.4 Chinenov et al., 2001;van Dam et al., 2001;Malik et al., 2017 Immune signaling and regulation chad chondroadherin ENSDARG00000045071 7.5 ccl39.1 chemokine (C–C motif) ligand 39, duplicate 1 ENSDARG00000101041 7.3 cxcl11.1 chemokine (C-X-C motif) ligand 11, duplicate 1 ENSDARG00000100662 6.0 cxcl8a chemokine (C-X-C motif) ligand 8a ENSDARG00000104795 5.3 si:dkey-117a8.4 c3a anaphylatoxin chemotactic receptor-like ENSDARG00000097698 5.0 ccl34b.8 chemokine (C–C motif) ligand 34b, duplicate 8 ENSDARG00000093098 4.0 rxfp1 relaxin family peptide receptor 1 ENSDARG00000090071 3.9 Figueiredo et al., 2006 CABZ01001434.1 ENSDARG00000098602 3.8 BX323596.1 C-X-C motif chemokine 11-6-like ENSDARG00000101138 3.8 cxcl11.8 chemokine (C-X-C motif) ligand 11, duplicate 8 ENSDARG00000095747 3.8 bmp5 bone morphogenetic protein 5 ENSDARG00000101701 3.7 Shih et al., 2017;Rosendahl et al., 2002 selenou1b selenoprotein U1b ENSDARG00000087059 3.6 Guo et al., 2015;Avery et al., 2018 sox2 SRY (sex determining region Y)-box 2 ENSDARG00000070913 3.4 egr3 early growth response 3 ENSDARG00000089156 3.3 Kenyon et al., 2017 Other role in innate immune response si:dkey-21e2.15 mast cell protease 1A ENSDARG00000092788 5.9 foxq1a forkhead box Q1a ENSDARG00000030896 5.5 Earley et al., 2018 p2ry11 purinergic receptor P2Y, G-protein coupled, 11 ENSDARG00000014929 3.8 Berchtold et al., 1999;Adrian et al., 2000 cyp21a2 cytochrome P450, family 21, subfamily A, polypeptide 2 ENSDARG00000037550 3.5 Poliani et al., 2010 MFAP4 (1 of many) microfibril-associated glycoprotein 4-like ENSDARG00000089667 3.1 Walton et al., 2015 Adaptive immune response RAC1 ras-related C3 botulinum toxin substrate 1-like ENSDARG00000099506 64.5 BX649608.2 ENSDARG00000102940 16.5 Mandel et al., 2012 ighv5-5 immunoglobulin heavy variable 5-5 ENSDARG00000096342 10.1 ighv4-1 immunoglobulin heavy variable 4-1 ENSDARG00000096259 9.8 CU896602.3 ENSDARG00000074999 7.0 tmem176l.3b transmembrane protein 176l.3b ENSDARG00000096874 6.0 Louvet et al., 2005;Zuccolo et al., 2010 egr1 early growth response 1 ENSDARG00000037421 5.7 Gomez-Martin et al., 2010 si:dkey-234i14.12 ENSDARG00000097228 5.4 unc119b unc-119 homolog b (C. elegans) ENSDARG00000044362 5.2 Gorska et al., 2004,2009;Gorska and Alam 2012 Stephen et al., 2018 ighv13-2 immunoglobulin heavy variable 13-2 ENSDARG00000096372 5.1 zgc:153659 ENSDARG00000039801 5.0 igl3v1 immunoglobulin light 3 variable 1 ENSDARG00000093258 4.0 S.-K.E. Harjula, et al. Developmental and Comparative Immunology 103 (2020) 103523 5 synthase 2 (alas2), choride channel K (clcnk) and nephrosis 1, congenital, Finnish type (nephrin) (nphs1) (Supplementary Fig. 1) as well as genes typically used as markers for different types of immune cells (Supplementary Fig. 2), indicating that the samples represent the kidney and different types of blood cells. 3.2. Innate immune response is dominant at 14 days post M. marinum infection The largest group of the protein coding genes induced upon infection consisted of genes with a documented or predicted immunological role. To further examine the induced immune response, these genes were divided into subcategories comprising acute phase response and Table 3 Metabolic and other non-immunological protein coding genes induced in mycobacterial infection in adult zebrafish kidney at 14 dpi (days post infection). The table shows the fold change in expression in M. marinum infected zebrafish compared to unchallenged fish (n= 4 in both groups). Fold change represents the fold change between the group medians. The table includes only the genes with at least two samples with ≥ 20 normalized reads after infection and the genes whose expression was induced at least 3.0-fold (calculated both with the DEseq-tool and from the medians). The fold changes are presented as the quotient of medians when the divisor is zero. Gene symbol Gene name/description Ensembl gene ID Fold change (median) Lipid binding and metabolism fabp6 fatty acid binding protein 6, ileal (gastrotropin) ENSDARG00000044566 131.5 dpep1 dipeptidase 1 ENSDARG00000068181 16.8 scd stearoyl-CoA desaturase (delta-9-desaturase) ENSDARG00000033662 6.3 abhd5b abhydrolase domain containing 5b ENSDARG00000100388 5.3 alox5b.2 arachidonate 5-lipoxygenase b, tandem duplicate 2 ENSDARG00000043089 4.4 nceh1a neutral cholesterol ester hydrolase 1a ENSDARG00000020427 3.6 crabp1b cellular retinoic acid binding protein 1b ENSDARG00000035904 3.4 g0s2 G0/G1 switch 2 ENSDARG00000078859 3.3 rbp7b retinol binding protein 7b, cellular ENSDARG00000070486 3.0 Other metabolic role CR559945.2 probable N-acetyltransferase CML5 ENSDARG00000106491 13:0 TPO thyroid peroxidase ENSDARG00000033280 43.0 ca4b carbonic anhydrase IV b ENSDARG00000042293 11.5 tgm5l transglutaminase 5, like ENSDARG00000098837 10.5 arid6 AT-rich interaction domain 6 ENSDARG00000069988 6.3 tgm8 transglutaminase 8 ENSDARG00000097651 6.0 ca9 carbonic anhydrase IX ENSDARG00000102300 5.8 lum lumican ENSDARG00000045580 5.6 hoxb5b homeobox B5b ENSDARG00000054030 4.1 si:dkey-203a12.9 ENSDARG00000104721 3.9 paplna papilin a, proteoglycan-like sulfated glycoprotein ENSDARG00000027867 3.7 si:ch211-243a20.3 si:ch211-243a20.3 ENSDARG00000092240 3.1 prrx1b paired related homeobox 1b ENSDARG00000042027 3.0 Other role si:dkey-7i4.24 myosin heavy chain, clone 203 ENSDARG00000096906 27.8 tnni2a.1 troponin I type 2a (skeletal, fast), tandem duplicate 1 ENSDARG00000045592 16.1 zgc:136930 thread keratin gamma ENSDARG00000055192 13.6 myhb myosin, heavy chain b ENSDARG00000001993 8.8 tcnbb transcobalamin beta b ENSDARG00000091996 7.4 hspb11 heat shock protein, alpha-crystallin-related, b11 ENSDARG00000002204 7.4 mrap melanocortin 2 receptor accessory protein ENSDARG00000091992 6.7 scpp1 secretory calcium-binding phosphoprotein 1 ENSDARG00000090416 6.7 slco2a1 solute carrier organic anion transporter family, member 2A1 ENSDARG00000061896 5.6 CABZ01072043.1 ENSDARG00000102907 5.6 tmem196a transmembrane protein 196a ENSDARG00000013935 5.4 hpcal1 hippocalcin-like 1 ENSDARG00000022763 5.1 asic2 acid-sensing (proton-gated) ion channel 2 ENSDARG00000006849 5.0 krt17 keratin 17 ENSDARG00000094041 4.3 cthrc1a collagen triple helix repeat containing 1a ENSDARG00000087198 4.1 col12a1b collagen, type XII, alpha 1b ENSDARG00000019601 4.1 sh3glb2b SH3-domain GRB2-like endophilin B2b ENSDARG00000035470 4.0 mamdc2b MAM domain containing 2b ENSDARG00000073695 3.9 si:ch211-105c13.3 ENSDARG00000089441 3.8 krt15 keratin 15 ENSDARG00000036840 3.8 tspan34 tetraspanin 34 ENSDARG00000103951 3.7 apnl actinoporin-like protein ENSDARG00000090900 3.4 myoz3a myozenin 3a ENSDARG00000067701 3.3 FO834829.3 ENSDARG00000102718 3.3 im:7150988 ENSDARG00000098058 3.3 BX664721.3 endonuclease domain-containing 1 protein-like ENSDARG00000073995 3.2 si:ch211-150d5.3 von Willebrand factor A domain-containing protein 7-like ENSDARG00000093384 3.0 Unknown role FO704758.1 ENSDARG00000098478 6.3 CABZ01056516.1 ENSDARG00000102467 5.0 si:dkey-248g15.3 ENSDARG00000097959 4.2 si:ch211-113d11.8 ENSDARG00000105494 3.4 S.-K.E. Harjula, et al. Developmental and Comparative Immunology 103 (2020) 103523 6 antimicrobial activity (4 genes), neutrophil degranulation (5 genes), inflammasome (3 genes), immune signaling and regulation (14 genes), other role in innate immune response (5 genes) and adaptive immune response (12 genes) (Fig. 1C, Table 2). Most of these subcategories were related to innate immunity, highlighting that these responses still have an important role in the immune defense at 14 days post M. marinum infection and in the initiation of the adaptive immune response. The innate immunity processes associated with defense response to bacteria were represented in our data. Acute phase protein coding gene serum amyloid A (saa) was induced together with three genes encoding proteins with antimicrobial activity: hepcidin antimicrobial peptide (hamp), liver-expressed antimicrobial peptide 2 (leap2) and predicted lysozyme g-like 2 (lygl2). Out of the other innate defense mechanisms, neutrophil degranulation was evidenced by the induction of five protein-coding genes. An inflammatory response against pathogenic M. tuberculosis includes inflammasome activation (Wawrocki and Druszczynska, 2017). Three of the induced genes were associated with inflammasomes. These were: si:ch211-236p5.3 encoding predicted NACHT, LRR and PYD domains-containing protein 3-like protein; si:ch211-233m11.1 coding for a putative NOD-like receptor, NACHT, LRR and PYD domains-containing protein 12 (Hu et al., 2017;Vladimer et al., 2012); and v-fos FBJ murine osteosarcoma viral oncogene homolog Ab (fosab). The latter codes for an ortholog of a transcription factor that is a part of activator protein 1 (AP-1) protein complex (Chinenov and Kerppola, 2001;van Dam and Castellazzi, 2001) which participates in inflammasome activation (Malik and Kanneganti, 2017). Also related to inflammasome activation, another upregulated gene, egr3, codes for Early growth response 3 transcription factor that positively regulates interleukin 1, beta, (il1b) expression (Kenyon et al., 2017). Some genes upregulated in our data, or their orthologs, are expressed in various innate immune cell types. These included genes expressed in macrophages (forkhead box Q1a,foxq1a (Earley et al., 2018) and one putative microfibril-associated glycoprotein 4-like, MFAP4 variant (Walton et al., 2015)); dendritic cells (purinergic receptor P2Y G-protein coupled, 11,p2ry11 (Berchtold et al., 1999) and putative cytochrome P450, family 21, subfamily A, polypeptide 2,cyp21a2 (Poliani et al., 2010)); and mast cells (putative mast cell protease 1A, si:dkey-21e2.15). Related to cell migration, genes and putative genes encoding chemokines and chemokine-like proteins formed the largest subset of the genes with a postulated or documented role in immune signaling and regulation. On top of this, the human ortholog of upregulated signaling protein relaxin family peptide receptor 1 (rfxp1), LGR7, has a reported role in leukocyte migration (Figueiredo et al., 2006). In addition to the genes involved in innate immune response, twelve genes coding for proteins with a suggested or documented role in adaptive immunity were upregulated upon M. marinum infection. Four of these: predicted CU896602.3,si:dkey-234i14.12,zgc:153659 and immunoglobulin light 3 variable 1 (igl3v1), code for proteins having a postulated role in immunoglobulin production. Five potential genes (or their human orthologs) have suggested roles in T and B lymphocyte signaling and regulation: early growth response 1 (egr1) (Gomez-Martin et al., 2010), immunoglobulin heavy variable 4-1 (ighv4-1), immunoglobulin heavy variable 5-5 (ighv5-5), immunoglobulin heavy variable 13-2 (ighv13-2) and RAC1 coding for a putative ras-related C3 botulinum toxin substrate 1-like protein. 3.3. M. marinum infection changes expression of metabolism associated genes In addition to the immune system-related genes presented in Table 2, nine genes involved in lipid binding or lipid metabolism and 13 genes classed as ‘other metabolic role’ were induced (Table 3). Significantly, the most induced gene, fatty acid binding protein 6 (fabp6) has a postulated role in lipid binding and metabolism. The upregulation of lipid metabolism associated genes due to M. marinum infection is supported by the literature, as the ability of M. tuberculosis to manipulate host metabolism and to use the host lipids for its pathogenic activities Fig. 1. Categorization of mycobacterium-responsive genes. Differentially expressed genes at 14 days post a low-dose (5–9 CFU, colony forming units) M. marinum infection were divided into categories based on their roles and/or the biological processes they participate in. (A, C) The upregulated genes were classified based on their roles as detailed in Materials and methods. (B) Downregulated genes were classified according to the gene ontology process data as detailed in Materials and methods. (C) Upregulated immunological genes were classified as in panel (A). S.-K.E. Harjula, et al. Developmental and Comparative Immunology 103 (2020) 103523 7 and survival has been widely studied in a number of different research models (Barisch and Soldati, 2017;Korb et al., 2016;Stutz et al., 2018). Involvement of lipid metabolism in the response to M. marinum is also apparent when analyzing the set of downregulated genes. This group included nine genes under the gene ontology term ‘lipid metabolism’, including 5 genes encoding apolipoproteins or predicted apolipoproteins and one, apolipoprotein Bb, tandem duplicate 2 (apobb.1), especially associated with lipid transport. Among the most enriched processes of the downregulated genes were lipid localization, cellular lipid catabolic process and lipid transport (Supplementary Table 2). The downregulation of lipid catabolism and transport could be associated with the tendency of mycobacteria to manipulate the host to store the lipids inside the cell for the bacteria to use for its own purposes (Peyron et al., 2008;Russell et al., 2009;Stutz et al., 2018). 3.4. A forward genetic screen identifies a fish line with increased susceptibility to M. marinum infection In order to identify genes that are important for defense against M. marinum infection in zebrafish, we conducted a forward genetic screen using gene-breaking transposon-based mutagenesis. 127 zebrafish lines, carrying a set of random mutations in F3–F5 inbred generations, were screened for altered susceptibility to M. marinum infection. To do this, we infected 3–15 month-old zebrafish from each mutant line with a low-dose (1–72 CFU) of M. marinum and followed the survival of the fish for 12–14 weeks. rag1 hu1999/hu1999 mutant fish, which have no functional T and B cells (Wienholds et al., 2002) were used as a positive control. During the screen, we identified 10 lines with impaired survival and one line with improved survival against M. marinum infection. The most susceptible line, designated mutant463, showed drastically impaired survival (8% at the 14 weeks post infection, wpi endpoint) compared to the WT zebrafish (74%, P< 0.0001) (Fig. 2A). For reference, Fig. 2A also shows the survival rates of two mutant lines as compared to WT; mutant189 (endpoint survival 83%) and mutant229 (endpoint survival 68%). Supplementary Fig. 3A shows the survival of mutant463 fish prior to RNA extraction for the RNA sequencing experiment (survival at the 4 wpi endpoint of 38% in the mutant463 group vs. 100% in the WT group, P= 0.002). This survival was also similar to the mutant463 zebrafish survival kinetics observed during the screen. Buffer injection (PBS) did not cause any mortality among the mutant463 fish (Supplementary Fig. 3A). Infection with 220–328 CFU of S. pneumoniae (T4 serotype) resulted in a minor impairment to mutant463 embryo survival when compared to WT embryo survival (Supplementary Fig. 3B). In our previous studies, high mortality after M. marinum infection associated with high bacterial burden (Harjula et al., 2018a;Myllymäki et al., 2018;Parikka et al., 2012). To clarify whether a defect in tolerance or resistance underlies the decreased survival of the mutant463 fish, we studied the M. marinum burden in the mutants and WT control zebrafish at 7 and 14 dpi (Fig. 2B). At 7 dpi, there was no difference in the bacterial count between mutant463 and WT fish (estimated copy number median 206 vs. 440 CFU, respectively), whereas at 14 dpi there was significantly higher bacterial burden in mutant463 fish compared to WT fish (copy number median 23,076 vs. 701 CFU, P= 0.042, respectively). To further study whether the decreased survival and the increased mycobacterial burden in mutant463 fish is attributed to defects in the innate immune response, we performed a survival assay on WT and mutant463 zebrafish embryos, which rely solely on the innate immune system for protection against infection. We infected 32–198 CFU into the caudal vein of the embryos at 1 dpf and followed their survival for 7 days (Supplementary Fig. 3C). In contrast to the situation in adult fish, mutant463 embryos showed similar survival kinetics to WT embryos. This result suggests that the increased susceptibility of the mutant463 zebrafish line against a mycobacterial infection is due to defective adaptive rather than innate immunity. 3.5. The hypersusceptible mutant463 zebrafish have decreased expression of si:ch211-236p5.3 and unc119b involved in the innate and adaptive immune response To further elucidate the reasons behind increased susceptibility of the mutant463 fish to mycobacterial infection, we analyzed the transcriptome of mutant463 zebrafish at 14 days post M. marinum infection and compared it to WT controls. Out of the 96 upregulated protein coding genes in WT fish, 27 had reduced expression in mutant463 compared to WT fish (Table 4). In addition, out of the 21 non-coding upregulated RNAs, 8 had reduced expression in the mutant463 zebrafish (Supplementary Table 6). None of the typical immune cell marker genes were clearly (≥ 3-fold) downregulated in mutant463 zebrafish compared to WT (Supplementary Fig. 2). However, there were other protein coding genes with a potential immunological role and reduced expression in mutants. For selected genes, the RNA sequencing result was confirmed by qPCR at 7, 14 and 28 dpi (Fig. 3). Among the genes with reduced expression in mutant463 fish were four genes with metabolic roles. One of these, fabp6, showed a 4.5-fold reduction at 14 dpi compared to WT fish in RNA sequencing and a 6.1fold reduction at 28 dpi in qPCR (P< 0.016) (Fig. 3A). In addition to its roles in lipid binding and metabolism, Fabp6 is a bile acid binding protein (Capaldi et al., 2009). Bile acid mediates inflammation through different mechanisms in the liver (Li et al., 2017) and attenuates Fig. 2. Survival and bacterial burden of the mutant463 zebrafish after M. marinum infection. Adult mutant463 zebrafish have impaired survival and elevated bacterial burden after a low-dose M. marinum infection. (A) The survival of WT (wild type) (n= 84), rag1 hu1999/hu1999 (n= 29,a positive control) and three mutant zebrafish lines, mutant189 (n= 29), mutant229 (n= 25) and mutant463 (n= 61), was monitored for 14 weeks post a low-dose (2–35 CFU, colony forming units) M. marinum infection. The graphs represent the data from one (mutant189 and mutant229), two (rag1 hu1999/hu1999 ) or three (WT, mutant463) independent experiments. The survival data is presented as a KaplanMeier survival curve. wpi, weeks post infection (B) The M. marinum burden in the abdominal organ blocks (excluding kidney) of the low-dose (5–9) infected WT (n= 8–12) and mutant463 (n= 8–12) adult zebrafish was measured with qPCR at 7 and 14 dpi (days post infection). The bacterial count is presented as a scatter dot plot and as the median of total bacterial copies (log10). A log-rank (Mantel-Cox) test for panel (A) and a one-tailed Mann-Whitney test for panel (B) was used to perform the statistical comparison of differences. S.-K.E. Harjula, et al. Developmental and Comparative Immunology 103 (2020) 103523 8 inflammation through its receptors in the kidney (Herman-Edelstein et al., 2018). fabp6 as well as another lipid metabolism associated gene stearoyl-CoA desaturase (delta-9-desaturase) (scd) were significantly induced upon M. marinum infection also when measured by qPCR. fabp6 was induced 131.5-fold in RNA sequencing at 14 dpi and in qPCR; 14.9fold at 7 dpi (P= 0.007) and 25.1-fold at 14 dpi (P= 0.015) (Fig. 3A). scd was induced 6.3-fold in RNA sequencing at 14 dpi and 2.1-fold (P= 0.007) in qPCR at 14 dpi (Fig. 3B). Out of the 27 protein coding genes with reduced expression in mutant463 fish, 7 have a documented or predicted role in innate immunity (Table 4). For example, the expression of predicted NACHT, LRR and PYD domains-containing protein 3-like (si:ch211-236p5.3) was almost completely absent in mutant463 zebrafish upon M. marinum infection, whereas it was induced by infection in WT fish (13.3-fold induction in RNA sequencing). The diminished expression of this gene in mutant463 compared to WT fish was also seen in qPCR analysis with larger sample sizes (Fig. 3C) at 14 dpi (1848.4-fold reduction, P< 0.0001) and at 28 dpi (493.7-fold reduction, P= 0.0002). This gene and its paralogs have been predicted to code for a component of the inflammasome and to have a role in intracellular microbial recognition and immune activation. Two of the genes with reduced expression have a role in adaptive immunity (Table 4). While ighv5-5 (3.0-fold reduction) has a predicted role in the B cell receptor signaling pathway and in the positive regulation of B cell activation, unc-119 homolog b (C. elegans) (unc119b) is a zebrafish ortholog for the human gene coding for Unc-119 lipid binding chaperone (UNC119), which has been shown to be essential for T cell activation (Gorska et al., 2004;Stephen et al., 2018). unc119b showed 23.9-fold reduction in RNA sequencing and in qPCR showed the following; 11.4-fold reduction at 7 dpi (P= 0.019), 8.4-fold reduction at 14 dpi (P= 0.003) and 76.0-fold reduction at 28 dpi (P= 0.0002) (Fig. 3D). 4. Discussion The transcriptional responses to mycobacterial infection are an outcome of the dynamic interplay between the host and the pathogen. Thus, it is implicit that the gene expression profile over the duration of M. marinum infection in zebrafish depends on the type of cells analyzed as well as the time-point of the analysis, the host and bacterial strains and the bacterial dosage, particularly when the analysis is carried out at the whole organismal level (Benard et al., 2016;Kenyon et al., 2017; Ojanen et al., 2019;Rougeot et al., 2014;van der Sar et al., 2009). Here, we utilized a zebrafish M. marinum infection model to study response to infection at the transcriptional level, using kidney samples to obtain a focused view of transcriptional changes in hematopoietic cells. Since head kidney is the main hematopoietic organ of zebrafish, it contains immune precursor cells. In our previous studies we have reported the presence of granulomas (Oksanen et al., 2013;Parikka et al., 2012) in adult zebrafish kidney. Therefore our data can give information about the interaction between the mature host immune cells and bacteria in addition to the production of immune cells during infection. As expected, the expression of many genes involved in the innate immune response were induced, such as genes involved in bacterial elimination, immune cell migration and neutrophil degranulation. Activation of an inflammasome is a typical response to intracellular Table 4 The protein coding genes downregulated in mycobacterial infection in mutant463 zebrafish kidney compared to WT (wild type) at 14 dpi (days post infection). The table shows the fold change in expression of the genes upregulated in mycobacterial infection between mutant463 and WT zebrafish (n= 4 in both groups). Fold change represents the fold change between the group medians. The table includes only the genes with at least two samples with ≥20 normalized reads in the WT group and whose expression was reduced at least 3.0-fold (calculated both with the DEseq-tool and from the medians). The fold changes are presented as the quotient of medians when the divisor is zero. Gene symbol Gene name/description Ensembl gene ID Fold change (median) Immune response si:ch211-236p5.3 NACHT, LRR and PYD domains-containing protein 3-like ENSDARG00000086418 −20:0 unc119b unc-119 homolog b (C. elegans) ENSDARG00000044362 −23.9 plaub plasminogen activator, urokinase b ENSDARG00000039145 −12.5 chad chondroadherin ENSDARG00000045071 −7.5 leap2 liver-expressed antimicrobial peptide 2 ENSDARG00000104654 −4.5 si:dkey-21e2.15 mast cell protease 1A ENSDARG00000092788 −4.2 ighv5-5 immunoglobulin heavy variable 5-5 ENSDARG00000096342 −3.0 Metabolism TPO thyroid peroxidase ENSDARG00000033280 −43 CR559945.2 probable N-acetyltransferase CML5 ENSDARG00000106491 −13 hoxb5b homeobox B5b ENSDARG00000054030 −4.6 fabp6 fatty acid binding protein 6, ileal (gastrotropin) ENSDARG00000044566 −4.5 Other role si:ch211-150d5.3 von Willebrand factor A domain-containing protein 7-like ENSDARG00000093384 −24 myhb myosin, heavy chain b ENSDARG00000001993 −22 tmem196a transmembrane protein 196a ENSDARG00000013935 −21.6 mamdc2b MAM domain containing 2b ENSDARG00000073695 −20.8 tspan34 tetraspanin 34 ENSDARG00000103951 −16 tnni2a.1 troponin I type 2a (skeletal, fast), tandem duplicate 1 ENSDARG00000045592 −14.8 myoz3a myozenin 3a ENSDARG00000067701 −13.1 scpp1 secretory calcium-binding phosphoprotein 1 ENSDARG00000090416 −8.0 slco2a1 solute carrier organic anion transporter family, member 2A1 ENSDARG00000061896 −4.5 tcnbb transcobalamin beta b ENSDARG00000091996 −4.4 lum lumican ENSDARG00000045580 −4.3 BX664721.3 endonuclease domain-containing 1 protein-like ENSDARG00000073995 −4.1 krt15 keratin 15 ENSDARG00000036840 −3.5 FO834829.3 ENSDARG00000102718 −3 Unknown role FO704758.1 ENSDARG00000098478 −33.3 CABZ01056516.1 ENSDARG00000102467 −6.4 S.-K.E. Harjula, et al. Developmental and Comparative Immunology 103 (2020) 103523 9