Abstract The house fly Musca domestica, a cosmopolitan dipteran insect, is a significant vector for human and animal bacterial pathogens, but little is known about its immune response to these pathogens. To address this issue, we inoculated the larvae with a mixture of Escherichia coli and Staphylococcus aureus and profiled the transcriptome 6, 24, and 48 h thereafter. Many genes known to controlling innate immunity in insects were induced following infection, including genes encoding pattern recognition proteins (PGRPs), various components of the Toll and IMD signaling pathways and of the proPO-activating and redox systems, and multiple antimicrobial peptides. Interestingly, we also uncovered a large set of novel immune response genes including two broad-spectrum antimicrobial peptides (muscin and domesticin), which might have evolved to adapt to house-fly's unique ecological environments. Finally, genes mediating oxidative phosphorylation were repressed at 48 h post-infection, suggesting disruption of energy homeostasis and mitochondrial function at the late stages of infection. Collectively, our data reveal dynamic changes in gene expression following bacterial infection in the house fly, paving the way for future in-depth analysis of M. domestica's immune system. Introduction Although lacking acquired immune systems, insects have efficient and potent innate immune systems to discriminate and combat foreign invaders successfully [35][1], [36][2]. It is generally acknowledged that the insect immune system involves cellular and humoral immune reactions against microbial infections that maintain close networks with each other and occur first in the epidermis, gut and tracheal respiratory organs and then in the hemocoel [37][2]. One characteristic of insect immunity is rapid activation of immune genes upon microbial infection, which produces effectors such as antimicrobial peptides. Insects have evolved sensitive mechanisms for recognition of pathogens including bacteria, fungi, parasites and viruses, which subsequently trigger cellular immune [38][3], [39][4] and humoral immune reactions [40][5], [41][6] via signal transduction pathways. Four signal transduction pathways, Toll, IMD, JNK and JAK/STAT, are regarded as the main pathways regulating the immune response of insects [42][7]. Genes encoding effectors are activated through signaling cascades and a set of these molecules are produced in specific tissues and secreted into the hemolymph [43][1]. House flies Musca domestica are endemic, and are the carriers of more than 100 harmful pathogens that can have severe consequences for human and animal health. Unfortunately, controlling the human diseases transmitted by house flies has not been successful due to the lack of knowledge of the basic molecular mechanism of this species [44][8]. Adaptation to distinct ecological environments might result in the evolution of specific immunity of house flies. Therefore, comparing the innate immune systems of Musca with those of the species that face different ecological pressures and pathogens such as Drosophila and Anopheles can be very informative, and thus offer clues on how house flies can flourish in close contact with many pathogens [45][8]. Recently, next generation sequencing technologies, such as the 454 Life Sciences (Roche) pyrosequencing platform, the Illumina Genome Analyzer, and the Applied Biosystems Solid platform provide rapid and high-throughput methods of identifying differentially expressed genes and their expression profiles [46][9], [47][10]. Identification and characterization of the host genetic factors released in response to pathogens is essential for understanding of innate immunity of M. domestica. However, information on the host genes involved in antibacterial defense is still limited. In this study, we performed transcriptome analysis and digital gene expression profile analysis of M. domestica challenged with Escherichia coli and Staphylococcus aureus, using high-throughput sequencing methods (Illumina Solexa Sequencing). The aims of this study were to uncover some information about the house fly immune response and discover new genes involved in bacterial infection in order to better understand the bacteria-host interaction. At the same time, the high-throughput sequencing in this study will identify a large number of transcripts that are comparable to the available transcripts in other species, and provide strong support for the genomic analysis of M. domestica. Materials and Methods Fly maintenance and bacterial challenge experiments The laboratory colony of M. domestica used in this study was a gift from Miss Fengqin He, Institute of Zoology, Chinese Academy of Sciences. Musca larvae were raised on artificial diet consisting of bran and water until pupariation. After eclosion, adult flies were fed with water, sugar, and milk powder. Specimens at all life stages were kept in a temperature-controlled room at 25±1°C, 70±5% relative humidity, and a photoperiod cycle of LD12:12. Septic injury was produced by pricking the abdomen of the 2^nd instar larvae with a needle previously dipped into a concentrated mixed bacteria suspension of E. coli and S. aureus [48][11]. The bacterial challenged larvae were maintained at 25°C on fresh medium for 6 h, 24 h, and 48 h before RNA extraction. RNA isolation Total RNA was isolated from the following developmental stages: eggs, 1^st instar larvae, 2^nd instar larvae, 3^rd instar larvae, pupae and newly-emerged adults (within 3 days of eclosion) in a 1∶1 female∶male ratio. RNA was extracted using the RNAiso Plus Kit (TaKaRa) following the manufacturer's instructions. Construction of the cDNA library and Illumina sequencing for transcriptome analysis Briefly, 20 mg total RNA (a mixture of RNA from eggs, 1^st instar larvae, 2^nd instar larvae, 3^rd instar larvae, pupae, adults, and bacterial challenged 2^nd instar larvae at equal ratios) was used to construct a cDNA library. Poly (A) mRNA was purified from total RNA using oligo (dT) magnetic beads. Fragmentation buffer was added for resizing mRNA to short fragments. Taking these short fragments as templates, random hexamer-primer was used to synthesize the first-strand cDNA. The second-strand cDNA was synthesized using buffer, dNTPs, RNaseH and DNA polymerase I, respectively. Short fragments were purified with the QiaQuick PCR extraction kit and resolved with EB buffer for end reparation and adding poly (A). After that, the short fragments were connected with sequencing adapters. And, after the agarose gel electrophoresis, the suitable fragments were selected for the PCR amplification as templates. Finally, the library could be sequenced using Illumina HiSeq™ 2000. Bioinformatics analysis Transcriptome de novo assembly was carried out with short reads assembling program Trinity [49][12]. Trinity first combines reads with a certain length of overlap to form longer fragments, which are called contigs. The reads are then mapped back to contigs. The paired-end reads enable the software to detect contigs from the same transcript as well as the distances between these contigs. Next, Trinity connects the contigs, and gets the sequences that cannot be further extended on either end. Such sequences are defined as unigenes. The unigenes were lengthened by the increased sequence depth using normalized library with the 454 platform [50][13]. Finally, BLASTx alignment (E-value<10^−5) between unigenes and protein databases like nr, Swiss-Prot, KEGG and COG (Clusters of Orthologous Groups of proteins) were performed, and the best aligning results were used to decide sequence direction of unigenes. When different databases gave conflicting results, we prioritized them in the following order: nr, then Swiss-Prot, KEGG, and COG. When unigene did not align with any of the entries in these databases, ESTScan [51][14] was used to predict the unigene's coding regions and to determine its sequence direction. Candidates of immunity-related genes from the fruit fly were used to query the M. domestica transcriptome. The nucleotide CDS and protein amino acid sequence for each of the identified Drosophila melanogaster immune genes were downloaded from the flybase ([52]http://flybase.org/). The M. domestica transcriptome was searched for homologues to these sequences using CLC Main Workbench 5.5 with the tBLASTn program with a cutoff E-value of 10^−5. Tentative matches were checked by searching the nr NCBI database using BLASTn for gene prediction errors. We collected the unigenes that are homologous to the important known innate-immunity-relevant genes, and aligned them to the reference gene sequences of D. melanogaster for counting the gene variant numbers. At first, we included the reference gene variants of a certain D. melanogaster gene together with the homologous unigenes from our data to do the alignment, and then we assigned the unigenes to their most homologous reference gene variants respectively. As for those none-full-length unigenes, sometimes we couldn't tell which reference gene variant they belong to because of their lacking of the specific sequence region, so that we assigned them to the first reference gene variant showed in the alignment result. Finally we counted the least variant number of the unigenes under each reference gene variant, and added up the numbers as the least variant gene number of the gene in the house fly. When we counted the least variant number of the unignenes, in order to avoid counting the none-overlapping fragments from the same variant, we only took into account the unigenes that overlap with each other and cover the highest phylogenetic diversity site in the alignment. Gene expression library preparation and sequencing Total RNA was extracted separately from immune-challenged larvae at 6, 24 and 48 h post-infection following the manufacturer's instruction, as described above. RNA extracted from non-challenged larvae at each matching time point was taken as a control. Ten larvae were collected for RNA extraction from each group. Next, a gene expression library was prepared using an Illumina gene expression sample prep kit. Briefly, mRNA was enriched using the oligo (dT) magnetic beads from the total RNA, and cDNAs were synthesized as described above for each sample. The library products were then ready for sequencing analysis via Illumina HiSeq™ 2000 using paired-end technology in a single run. Six libraries from control and challenged groups were constructed. Analysis and annotation of gene expression tags The 50 bp length reads were gained through base calling and were filtered to get the clean reads by removing the dirty reads with adaptors, reads in which unknown bases are more than 10%, and low quality reads (the percentage of the low quality bases of quality value ≤5 is more than 50% in a read). Then, clean reads were mapped to the above transcriptome reference database using SOAPaligner/soap2 [53][15], allowing no more than 2 mismatches. Gene expression levels were evaluated with RPKM values [54][16], and the value was substituted by 0.001 if the gene has no expression. NOIseq method that can work in absence of replication was used for gene expression analysis [55][17]. If there was more than one transcript for a given gene, the longest transcript was used to calculate its expression level and coverage. To identify differentially expressed genes between two samples, the false discovery rate (FDR) method was used to determine the threshold of P-value in multiple tests [56][18]. We use “FDR≤0.001 and the absolute value of log[2]Ratio ≥1” as the threshold to judge the significance of gene expression difference. In the gene expression profiling analysis, Gene Ontology (GO) enrichment analysis of functional significance was applied using the hypergeometric test to map all differentially expressed genes to terms in the GO database, looking for significantly enriched GO terms in differentially expressed genes, and comparing them to the transcriptome database. For the pathway enrichment analysis, we mapped all differentially expressed genes to terms in the KEGG database and looked for significantly enriched metabolic pathways or signal transduction pathways in differentially expressed genes comparing with the whole transcriptome database. Then, all genes expression was subjected to KEGG enrichment analysis compared to the transcriptome background using a hypergeometric test. Annotated KEGG pathways with a Q value≤0.05 were considered as significantly enriched terms in differentially expressed genes. qPCR validation The accuracy of the genes expression approach has been validated by qPCR analysis. Total RNA was extracted from bacterial challenged and non-challenged larvae as described for gene expression library preparation, and then was reverse-transcribed according to the protocol provided with M-MLV reverse transcriptase (Promega, USA). Ten differentially expressed genes were randomly chosen to verify gene expression sequencing results using four replicates. The β-actin gene was used as a constitutive expression control for normalization. The primers are shown in [57]Table S1. qPCR was performed following the above-mentioned methods with modified annealing temperature for each pair of primers [58][19]. The relative quantification (comparative method) was calculated using the ΔΔCt method [59][20]. All samples were normalized to the ΔCt value of a reference gene to obtain a ΔΔCt value (ΔCt target − ΔCt reference). The final relative expression was calculated using the following formula: F = 2^−ΔΔCt. The data obtained from qPCR were analyzed for statistical significance using Graph-Pad Prism [60][21]. The significance at P<0.05 was analyzed using one-way ANOVA. The qPCR results were then compared with gene expression data to detect the expression correlation of each gene. Peptide synthesis and antibacterial assays The peptides muscin and domesticin, were chemically synthesized using the solid-phase method on the Applied Biosystems model 430A peptide synthesizer by Shanghai mocell biotech Co., Ltd. (Shanghai, China). Synthesized peptides were purified by high performance liquid chromatograph (HPLC), and verified by mass spectrometry and amino acid composition analysis. Each synthesized peptide was diluted with sterile deionized water to different concentrations, and applied directly to antibacterial assays. Twelve bacterial strains used in the tests were a gift from Shunyi Zhu, Institute of Zoology, Chinese Academy of Sciences (Beijing, China). Minimum inhibitory concentrations (MICs) were determined in duplicate by the liquid growth inhibition assay [61][22]. Briefly, 10 µl of each dilution (sterile deionized water as a control) were incubated in sterile microtiter plates (96 wells) with 100 µl of a suspension of midlogarithmic phase culture of bacteria at a starting optical density of OD[600] = 0.001, or with 80 µl of fungal spores (final concentration 10^4 spores/ml) and 10 µl of water. Poor-broth nutrient medium (1% bactotryptone, 0.5% NaCl w/v, pH 7.5) was used for standard bacterial cultures. Anti-fungal assays were performed in potatoes dextrose broth (Difco) at half-strength supplemented with 10 mg/ml tetracycline. Bacteria were grown during 24 h under vigorous shaking at 30°C. Fungi were grown at 25°C in the dark without shaking for 48 h in a moist chamber. Microbial growth was controlled by measurement of the optical density at D[600] after a-24 h incubation. Inhibition of filamentous fungi growth was observed at microscopic level after 24 h and measured at 600 nm after 48 h. The MIC values are expressed as the range between the highest concentration of the peptide where bacterial growth was observed and the lowest concentration that caused 100% of inhibition bacterial growth [62][23]. Results and Discussion Sequencing and sequence assembly To obtain the global gene expression profile of the house fly, RNA samples from M. domestica of six developmental stages (eggs, 1^st instar larvae, 2^nd instar larvae, 3^rd instar larvae, pupae, adults) were prepared, equally-mixed and then sequenced by the Illumina platform in a single run which generated 70,335,268 raw reads, 66,049,270 clean reads and 5,944,434,300 total clean nucleotides ([63]Table 1). The average read size, Q20 percentage (sequencing 1% error rate,) and GC percentage were 90 bp, 96.25%, and 53.69%, respectively. These short reads were assembled into 116,687 contigs with a mean length of 319 bp. Then, after gap filling of contigs using paired-end reads from the transcriptome sequencing data of 454, we obtained 47,086 unigenes. The mean size of these unigenes was 757 bp and the N50 was 1,186 bp. Here N50 was the length of the smallest contig in the set that contained the fewest or largest contigs whose combined length represents at least 50% of the assembly. Of these unigenes, 10,933 (23.2%) were larger than 1,000 bp ([64]Figure S1). Due to the short length of the sequencing read, unigenes obtained from Illumina platform were short in general. In present study, the unigenes were lengthened by the increased sequence depth using our previous normalized library with 454 platform [65][13]. Table 1. Summary of the M. domestica transcriptome. Parameters Number Total clean reads 66,049,270 Total Nucleotides (nt) 5,944,434,300 Total number of contigs 116, 687 Total number of unigenes 47,086 Mean size of unigenes (bp) 757 N50 of unigenes (bp) 1,186 Sequences matched known genes with E-value<10^−5 27,021 [66]Open in a new tab Unigene functional annotation Unigene sequences were annotated by searching the Nr of NCBI, Swiss-Prot, KEGG and COG protein database using BLASTx with a cut-off E-value of 10^−5. A total of 27,021 distinct sequences (57.39% of unigenes) matched known genes ([67]Table S2). The remaining 20,065 unigenes (42.61%) could not be done and needed more genetic data to annotate. It is noteworthy that a large part of the unigenes in M. domestica transcriptome database is with unknown functions. Based on sequence homology, 23,212 unigenes (85.9%) were annotated and divided into 25 different COG categories ([68]Figure 1). The general function category that contained 3,880 unigenes (16.7%) was the largest, followed by carbohydrate transport and metabolism (1,959, 8.4%), transcription (1,935, 8.3%), translation, ribosomal structure, and biogenesis (1,656, 7.1%), eplication, recombination and repair (1,498, 6.5%), and post-translational modification, protein turnover, and chaperones (1,329, 5.7%). Only two unigenes belonged to nuclear structure, which was the smallest group. Figure 1. Classification of the clusters of orthologous groups (COG) for the M. domestica transcriptome. [69]Figure 1 [70]Open in a new tab 8,549 unigenes (31.6% of the total annotated unigenes) were divided into 25 specific categories. The GO categories recovered from Blast2GO analyses were summarized by the proportion of unique sequences annotated in each GO level 2 category ([71]Figure 2). We categorized 7,063 unigenes (26.1% of total) into 48 function groups. Cell, cellular process, cell part, binding and metabolic process were the five largest groups, containing 3,812, 3,751, 3,389, 2,919, and 2,914 unigenes, respectively. Figure 2. Classification of the gene ontology (GO) for the M. domestica transcriptome. [72]Figure 2 [73]Open in a new tab 7,063 unigenes (26.1% of total) were categorized into 48 function groups. To identify the biological pathways that are active in M. domestica, functional classification and metabolic pathway assignment were performed by the KEGG annotation system. A total of 16,440 unigenes (60.8%) were classified into 239 KEGG pathways. Metabolic pathways contained 2,274 unigenes (13.8%) was significantly larger than other pathways, such as pathways in cancer (619, 3.77%), focal adhesion (565, 3.44%), and regulation of actin cytoskeleton (545, 3.32%). In addition, some immune related pathways, including phagosome, lysosome, ECM-receptor interaction, Fc gamma R-mediated phagocytosis, leuko-cyte trans-endothelial migration, complement and coagulation cascades, and many signaling transduction pathways such as MAPK, VEGF, JAK-STAT, PPAR, and Toll-like receptor signaling pathway were predicted in the KEGG database ([74]Table S3). These annotations provide a valuable resource for investigating specific processes, functions and pathways and allow for the identification of novel genes involved in the pathways of immunity. Annotation of immune-relevant genes Based on genome-wide analysis, many immune-related genes have been identified from D. melanogaster (265), Anopheles gambiae (304), Apis mellifera (138), and Bombyx mori (220) [75][2]. To gain deep insight into the molecular biology of immune systems in M. domestica, the immune-relevant genes were discovered by referencing to those identified in other insect species. Approximately 279 genes were found to be homologous to known immune-relevant genes ([76]Table 2), including the most important elements of innate immunity, such as pattern recognition receptors, and other immune-related genes (such as PPO, NOS, caspase, dicer2, argonaute2). It is noteworthy that many unigenes obtained by next generation sequencing should be different fragments or even allelic or splice variants of the same gene [77][24]. In this study, the number of gene sequences predicted to encode immune genes would be therefore overestimated over the actual number of genes belonging to each of the currently characterized gene families. The putative genes uncovered provide the basis for further understanding of the physiological functions of candidate genes in M. domestica immune responses. Table 2. Gene counts belonging to immunity-related gene families. Gene family Gene numbers Drosophila melanogaster [78]^a Anopheles gambiae [79]^a Apis mellifera [80]^a Bombyx mori [81]^a Musca domestica [82]^b Recognition PGRP 13 7 4 12 16 βGRP 3 7 2 4 3 fibrinogen-related protein 14 61 2 3 19 scavenger receptor 21 21 13 18 13 C-type lectin 34 25 10 21 16 hemocytin 1 0 1 1 1 galectin 6 8 2 4 7 TEP 6 15 3 3 36 nimrod 10 4 4 4 8 draper 1 1 1 1 2 eater 1 1 0 0 0 dscam 1 1 1 1 1 Modulation clip serine protease 37 41 18 15 6 serpin 30 17 5 26 40 Toll pathway spätzle 6 6 2 3 1 toll 9 10 5 14 14 MyD88 1 1 1 1 1 tollip 1 2 1 2 1 tube 1 1 1 1 1 pellino 1 1 1 1 1 pelle 1 1 1 1 1 TRAF2 1 1 1 1 0 ECSIT 1 1 1 1 2 cactus 1 1 3 1 2 Dif/Dorsal 2 1 2 1 2 IMD pathway IMD 1 1 1 1 1 DREDD 1 1 1 1 1 TAK1 1 1 1 1 1 FADD 1 1 1 1 1 TAB2 1 1 1 1 2 IAP2 1 1 1 1 1 IKK 2 2 2 2 1 UBC13 1 1 1 1 0 relish 1 1 2 1 1 JNK pathway HEM 1 1 1 1 1 JNK 1 1 1 1 3 FOS 1 1 1 1 1 JUN 1 1 1 1 1 JAK/STAT pathway UPD3 1 0 0 0 0 PIAS 1 1 1 1 1 SOCS 1 1 1 1 1 domeless 1 1 1 1 4 hopscotch 1 1 1 0 2 STAT 1 2 1 1 1 Antimicrobial peptide ceropin 4 4 0 13 3 attacin 4 1 0 2 9 diptericin 2 0 0 0 3 defensin 1 4 2 1 1 gloverin 0 0 0 4 0 moricin 0 0 0 9 0 lebocin 0 0 0 1 0 domesticin [83]^* 1 0 0 0 1 muscin [84]^* 0 0 0 0 1 Melanization PPO 3 9 1 2 7 DDC 1 7 1 1 2 DCE 2 0 18 16 4 TH 1 1 1 1 2 punch 1 0 1 0 1 Other effectors NOS 1 1 1 2 1 POI 2 1 0 1 0 lysozyme 11 8 2 4 13 Other immune molecules caspase 6 11 4 4 5 dicer2 1 1 1 0 6 argonaute2 1 1 1 1 2 Total 265 304 138 220 279 [85]Open in a new tab ^a Number of gene sequences from data of Christophides et al. (2002) [86][57], Evans et al. (2006) [87][7], Tanaka et al. (2008) [88][2]. ^b Number of gene sequences obtained in this study that annotated with NCBI nr database. * Two novel antimicrobial peptides identified and named in present study. Gene expression library sequencing Having established the house fly transcriptome database, we next determined how bacterial infection altered the transcriptome. We sequenced RNA from 2^nd instar larvae at 6 h, 24 h, 48 h following inoculation with a mixture of E. coli and S. aureus; RNA from uninfected larvae served as the control (designated as M6, M24, M48 and CK6, CK24, CK48, respectively). After filtering the dirty tags, each sample generated 3,510,021∼3,672,652 reads. Among these clean reads, 78.51∼86.51% were mapped to unigenes in the transcriptome database ([89]Table 3). The percentage of clean reads ranged from 96.78% to 99.29%, reflecting the high quality of the sequencing. About 4∼5% of genes were covered between 90∼100% in each library. Table 3. Statistics of gene expression sequencing. Summary M6 CK6 M24 CK24 M48 CK48 Total reads 3,661,408 3,510,021 3,513,957 3,661,248 3,560,749 3,672,652 Total base pairs 179,408,992 171,991,029 172,183,893 179,401,152 174,476,701 179,959,948 Total mapped reads 3,154,425 3,015,374 3,040,076 2,874,614 2,971,222 3,041,448 Perfect match 2,421,048 2,312,301 2,287,588 2,002,766 2,242,858 2,159,328 ≤2 bp mismatch 733,377 703,073 752,488 871,848 728,364 882,120 Unique match 2,216,969 2,097,216 2,069,330 2,134,543 2,166,611 2,164,870 Multi-position match 937,456 920,759 970,746 740,071 804,611 876,578 Total unmapped reads 506,983 494,647 473,881 786,634 589,527 631,204 [90]Open in a new tab Effects of bacterial infection on gene expression The “FDR≤0.001” and the absolute value of “log[2] Ratio≥1 or ≤−1” were used as the threshold to identify and compare differentially expressed genes between larvae non-challenged and challenged at different stages. Numbers of differentially expressed genes for each comparison were shown in [91]Figure 3. 572, 3194, and 3544 genes were significantly affected by bacterial infection at 6, 24 and 48 h post-infection, respectively ([92]Table S4), They include genes involved in insect immunity (ECM-receptor interaction, phagosome, complement and coagulation cascades, PPAR signaling pathway) and protein, carbohydrate and lipid metabolism. Besides, genes reflecting mitochondrial oxidative stress and endoplasmic reticulum (ER) stress were up-regulated. Thus, infection caused extensive changes in the M. domestica physiological status. Figure 3. Differentially expressed genes of each library compared with CK. Figure 3 [93]Open in a new tab Unigenes changed at transcriptional level following bacterial challenge were identified by filtering of the one-fold up- and down-regulated ones with FDR≤0.001. Importantly, some of these changes are temporal-specific. There were 12, 33, 43 pathways significantly enriched in the 3 comparisons respectively (Q value≤0.05). The results are shown in [94]Table S5. For example, at 6 h after bacterial infection, metabolic pathways were preferentially altered, with the genes affected including those controlling fatty acid, galactose and amino acid metabolism, ribosome biogenesis, phagocytosis, neuroactive ligand-receptor interaction, collecting duct acid secretion, terpenoid-quinone biosynthesis and adipocytokine signaling. In contrast, at later stages, particularly at 48 h post-bacterial challenge, genes controlling oxidative phosphorylation were strongly down-regulated, suggesting profound mitochondrial dysfunction at these stages. Differentially expressed genes involved in immune response Among the genes affected by infection, 509 were presumably involved in innate immune responses ([95]Table S6), including pattern recognition, Toll and IMD signaling, direct antimicrobial defense, proPO-activating cascade and redox, as described below. Pattern recognition proteins Invading microbes are detected by pattern recognition proteins (PRPs), which bind conserved pathogen-associated molecular patterns (PAMPs) shared by broad classes of microorganism. Several types of PRRs have been reported in invertebrates, such as peptidoglycan recognition proteins (PGRPs), β-1,3-Glucan recognition protein (βGRP), C-type lectins (CTLs), scavenger receptors (SRs), and thioester-containing proteins (TEPs) [96][25]. PGRP, first discovered in haemolymph of B. mori, binds bacterial peptidoglycan and triggers the prophenoloxidase cascade, culminating in the activation proPO and spätzle [97][26], [98][27]. PGRPs are considered to be the largest and most versatile family of pattern recognition molecules for microbial products in insects [99][28]. Indeed, we found that there are multiple putative PGRP unigenes in the house fly, each induced at all three time points post-infection, except that Unigene29467 was repressed at 6 h post-challenge. Of note, recent studies in Drosophila indicates that amidase PGRPs negatively regulates the IMD pathway by degrading PGN [100][29], [101][30], suggesting some late-expressed PGRPs in house fly may act to dampen immune response. In contrast, βGRP was hardly responsive to infection, with only one gene (Unigene30021) mildly (1.54x) up-regulated at only one time point (6 h post-challenge). βGRPs were first identified in B. mori and crayfish Pacifastacus leniusculus, for their ability to bind β-glucans, and role in β-glucan-induced activation of phenoloxidase [102][31]. These proteins are found only in invertebrates and contain a protein domain that is similar to several bacterial glucanases. Some of the proteins have broader, some narrower defense specificities, and all are associated with the regulation of immune signaling pathways [103][32]. C-type lectins (CTLs) are a large family of PRRs found in almost all metazoans and mainly exert their functions depending on a common structural motif, the carbohydrate recognition domain [104][33]. The members of CTL family are abundant in shrimp and have various functions in innate immunity, including phagocytosis, melanization, respiratory burst, agglutination, antibacterial and anti-viral responses [105][34]. A total of 21 CTL unigenes were identified in the M. domestica transcriptome data, and most of them were massively up-regulated in larvae during various stages of infection. Thus, CLT may be important for antibacterial defense. Scavenger receptors (SRs) recognize different PAMPs, including LPS, lipoteichoic acid (LTA) and yeast zymosan/β-glucan, and act as phagocytic receptors mediating non-opsonic phagocytosis of pathogens [106][35], [107][36]. We identified 73 unigenes 7 of them up-regulated at 48 h (the last stage post-infection) but none at earlier stages, raising the possibility that SRs might function to clean damaged biomolecules and even cell debris at late stage of infection. Thioester-containing proteins (TEPs) are structurally related to the complement C3/alpha (2)-macroglobulin family [108][37], [109][38]. One of the best characterized TEPs in invertebrate is the mosquito A. gambiae TEP1, where upon bacterial infection, TEP1 was cleaved to release the C-terminal fragment that promotes phagocytosis of bacteria [110][39]. We found 40 putative unigenes, 6 of them up-regulated at 6 h and 48 h post-challenge, consistent with their roles in innate immunity. Toll and IMD signaling pathways Toll and IMD pathways activate antimicrobial peptide genes and regulate the host humoral response [111][40]. We identified many components in these two pathways, including spätzle, toll, MyD88, tollip, tube, pellino, pelle, TRAF2, ECSIT, cactus, Dif/Dorsal, IMD, DREDD, TAK1, FADD, TAB2, IAP2, IKK and relish ([112]Table 2). A few unigenes encoding some of these proteins (e.g., relish, cactus, and toll) were up-regulated following bacterial infection with different expression patterns ([113]Table S6). Antimicrobial peptides A hallmark of the insect defense is the induction and secretion of antimicrobial peptides that accumulate in the hemolymph where they oppose invading microorganisms [114][40]. Twenty unigenes encoding 4 known antimicrobial peptides were found to be up-regulated drastically in the infected M. domestica larvae at all three time points post-bacterial infection, reinforcing their crucial roles in innate immunity ([115]Table S6), including cecropin, attacin, defensin and diptericin. We also found two novel dramatically up-regulated unigenes (Unigene29893 and Unigene13085), which we named muscin (GenBank: [116]KF514663) and domesticin (GenBank: [117]KF514664) ([118]Figure S2 and [119]Figure S3). The deduced peptides, containing putative signal peptides and rich in positively charged and hydrophobic amino acids, were synthesized and their antimicrobial activities tested using liquid growth inhibition assay. The two peptides were both active against all the tested bacteria ([120]Table 4). Muscin and domesticin were most effective against Serratia marcescens and Micrococcus luteus, respectively (MIC 1.25–2.5 and 0.1–0.2 µM), and weaker activities were observed against Salmonella typhimurium (MIC 50–100 and 25–50 µM for muscin and domesticin respectively). However, the peptides lacked any detectable activity against fungi Neurospora crassa even at 100 µM. No homolog of muscin was found through searching the GenBank database,while domesticin shows some similarity to a few peptides of Drosophila, including IM18 peptide, an immune-induced molecule with unknown function ([121]Figure S4) [122][41]. On the other hand, we failed to identify the homologus genes of gloverin, moricin, lebocin, antimicrobial peptides in B. mori. We predicted that the rapid evolution of antimicrobial peptide genes in insects is likely to be the result of host-pathogen co-evolution, indicating a specific role of antimicrobial peptides against a restricted subset of pathogens. Comparative analysis of gene repertoires of the antimicrobial peptides from different insects suggests that evolution of antimicrobial peptides followed independent scenarios as a result of specific adaptation to distinct ecological environments. Table 4. Antibacterial activities of muscin and domesticin. Microorganism MIC (µM) muscin domesticin Gram-positive bacteria Bacillus megaterium 2.5–5 0.2–0.4 Bacillus sp. DM-1 12.5–25 1.25–2.5 Bacillus subtilis 2.5–5 5–10 Micrococcus luteus 25–50 0.1–0.2 Staphylococcus aureus 12.5–25 2.5–5 Gram-negative bacteria Agrobacterium tumefaciens 25–50 5–10 Escherichia coli ATCC 25922 12.5–25 12.5–25 Pseudomonas aeruginosa 25–50 2.5–5 Serratia marcescens 1.25–2.5 0.5–1 Salmonella typhimurium CCTCCAB 94007 50–100 25–50 Xanthomonas oryzae 12.5–25 1.25–2.5 Fungi Neurospora crassa >100 >100 [123]Open in a new tab MICs are expressed as the interval a–b, where a is the highest concentration tested at which microorganisms are growing and b the lowest concentration that causes 100% growth inhibition. Antimicrobial peptide genes were dominant in up-regulated transcripts at all three time points post-bacterial infection. These results indicated that various antimicrobial peptides could respond rapidly to invaded pathogens and keep high-level expression during all these detected stages. The up-regulation of antimicrobial peptides expression after infection is a common defense strategy by hosts to rapidly destroy invaders. Differential expression patterns among different antimicrobial peptides may represent a combined strategy to form a defense network against diverse microbial pathogens. Lysozymes, present in many organisms, cleave the β-1, 4-glycosidic linkage between N-acetylmuramic acid and N-acetylglucosamine found in certain bacterial cell walls [124][42]. Their functions have been widely described in species that ingest or harbour bacteria throughout their life cycles, such as M. domestica and D. melanogaster, and also characterized as digestive enzymes [125][42]–[126][45]. Other reports have suggested a major role of lysozymes in immune responses to pathogens [127][42], [128][46]. We found that two of the four lysozyme unigenes up-regulated and two other down-regulated post-infection, suggesting Musca lysozymes may function as defense molecules as well as digestive enzymes. proPO-activating system An immediate immune response in insects is the cell-mediated melanization reaction observed at the site of cuticular injury or on the surface of parasites invading the hemocoel. Melanization requires the activation of proPO (prophenoloxidase), an enzyme that catalyzes the oxidation of mono- and diphenols to orthoquinones, which polymerize nonenzy-matically to melanin [129][40]. The proPO-activating system mainly includes many genes such as serine proteinases and their inhibitors (serpins), proPO-activating enzyme (PPA), proPO and its active form, phenoloxidase (PO) [130][47], [131][48]. After stimulated by injury or PAMPs, a serine proteinase cascade is first activated [132][49], which leads to the cleavage of the pro-form of the prophenoloxidase-activating enzyme (pro-PPA) into active PPA. Then PPA further activates the proPO into the active enzyme, PO, through proteolysis of its propeptide [133][47]. In the present study, many differentially expressed genes were annotated to be tentative members of the proPO-activating system ([134]Table S6). These genes were mainly a kind of serine proteinases, and their inhibitors. The present gene expression data revealed that most members in the serine proteinase cascade and proPO system were responsive to bacterial infection in M. domestica. The expression level of proPO and PPA transcripts were increased significantly at 24 h or 48 h post-bacterial challenge. The expressions of most serine proteases were increased in bacterial-challenged larvae, indicating a positive response of the serine proteinase cascade in the immune defense. However, the profiles of serpins were also up-regulated unexpectedly during this process, which seemed incompatible with their roles in regulation of the proPO-activating system. A similar consequence was shown in shrimp Fenneropenaeus chinensis [135][50]. It might be taken as a negative feedback mechanism to avoid damage of host tissues and cells by excess reactive components generated by PO. Redox system Reactive oxygen species (ROS) is referred to as a class of radical or non-radical oxygen-containing molecules that have high oxidative reactivity with lipids, proteins, and nucleic acids [136][51]. Generation of ROS in infected organisms not only helps to kill pathogens but also acts on the host cells themselves, including altering the intracellular redox balance and functioning as signaling molecules involved with the regulation of immunomodulatory genes. ROS production is regarded as an immediate acute-phase oxidative defense in response to pathogen assault or cellular stress such as phagocytosis and melanotic encapsulation [137][52]. Our analysis showed that the transcript levels of some cytochrome P450 and xanthine oxidase genes were strongly induced in the process of bacterial infection, especially at 48 h. Moreover, some unigenes encoding oxidoreductases, such as dehydrogenase/reductase SDR family member 11-like, 15-hydroxyprostaglandin dehydrogenaseor-related dehydrogenase, FAD-NAD binding oxidoreductase, peroxidasin-like, NADH dehydrogenase subunit 5 etc., were up-regulated significantly in this process. However, few unigenes encoding antioxidant enzymes were also seen up-regulations in this study. Previous study has revealed that M. domestica SOD genes could be induced mainly at 72 h post E. coli or S. aureus challenge [138][19]. We suppose that comprehensive up-regulated expressions of antioxidant enzyme genes would appear in the later stage. Furthermore, increasing sequencing depth in the future might increase both the number of immune related genes, and provide more data on their differential expression pre- and post-infection. There must be a fine redox balance maintained by innate immune system, which is critical for insect to survive from the war between pathogen and host. The underlying molecular mechanisms, however, still remain elusive. Our further research will focus on mitochondrion, an organelle which is considered as the central platform for innate immunity [139][53]. In addition to those known immune-related genes, we are surprised to find that hexamerins are up-regulated strongly at 24 and 48 h post-challenge ([140]Table S4). In general, hexamerin is thought to act as storage protein which is used as a source of amino acids and energy for protein synthesis during metamorphosis. But some reports indicated that insect hexamerin may be involved in innate immunity [141][54]–[142][56]. As humoral procoagulant, hexamerin can bind to the surfaces of bacteria invaded, and its subunits are confirmed as the major constituent of the aggregates in Drosophila hemolymph [143][54]. So we speculate that hexamerins might be involved in the innate immune response of the house fly. qPCR To validate the gene expression result, qPCR analysis was performed using gene-specific primers for ten changed unigenes selected at random. Results are shown in [144]Table 5 and are compared with the gene expression data. Although the qPCR data supports the trends of our gene expression results, the levels of change tend to be different between methods for each gene. We attribute differences in the level of change to the sensitivity of biases occurred between qPCR and gene expression. In general, the results of gene expression profiling are reliable. Table 5. Comparisons of relative gene expression fold between gene expression data and qPCR results. Name gene ID Fold change at 6/24/48 h (gene expression data) qPCR fold change at 6/24/48 h attacin Unigene12874 15.2 4.7 3.9 210.8±59.0** 52.6±15.8** 25.1±4.2** defensin CL7841.Contig1 15.9 3.4 6.0 220.3±60.5** 45.1±12.9** 63.2±16.1** diptericin Unigene13439 12.8 9.8 5.7 180.5±51.4** 112.3±30.8** 62.9±18.5** muscin Unigene29893 8.3 - 1.9 55.0±16.2** 16.0±4.3** 26.0±7.6** prophenoloxidase Unigene28739 - - 10.9 42.3±12.1** 51.6±10.9** 116.7±32.2** PGRP-SD CL486.Contig1 4.2 1.4 1.6 7.8±2.4** 1.5±0.3** 1.1±0.3* HSP67B2 CL8219.Contig1 - 1.1 1.2 1.0±0.3 0.9±0.2 1.0±0.2 eiger CL358.Contig2 - 1.5 3.5 2.5±0.7** 3.9±1.1** 4.8±1.1** metallothionein Unigene13124 - 1.2 - 1.0±0.3 1.0±0.2 3.1±0.9** SOD CL6208.Contig1 - - 1.9 1.9±0.5 1.1±0.3 2.1±0.5* [145]Open in a new tab Significant differences between the challenged and control group are indicated with * at P<0.05 or ** at P<0.01. Conclusions In this study, bacteria-challenged M. domestica transcriptome profiles were investigated and the substantial amount of transcripts was recovered, which provided a strong support for the future genomic research on M. domestica, especially on in-depth genome annotation in insects. Universally identified immune-candidate genes, infection markers, and putative signaling pathways were found in M. domestica and especially a considerable amount of immune-relevant genes and pathways in the house fly showed significant similarity to Drosophila, Anopheles, Apis, and Bombyx, suggesting that mechanisms underlying the innate immunity in insects might be conserved in invertebrates. After the bacterial challenge, pattern recognition proteins, especially the PGRPs, significantly increased. Antimicrobial peptides, including ceropin, attacin, diptericin, defensin, were also found increased. Moreover, component genes in Toll and IMD signaling pathways, as well as the genes involved in proPO-activating system and redox system, were strongly induced in the process of bacterial infection. In addition, a large set of novel immune response genes that have never been linked previously to immune responses in other insects indicated that the immune system of insect might be much more complex than previously reported, and house fly-specific immune events might have happened during evolution as a result of specific adaptation to distinct ecological environments. Particularly, we realized that regulations of whole-body energy and metabolic homeostasis were disrupted, and there was a strong suppression of oxidative phosphorylation during antibacterial defense reaction. We firmly believe that unclear repair and rebalance mechanisms at the later stage of infection should be crucial for insect to survive from the pathogen-host battle. Two novel antimicrobial peptides, muscin and domesticin, were found in challenged house flies, and they both showed broad spectrum of bactericidal activities. Domesticin is homologous to some peptides of Drosophila, including an immune-induced molecule IM18 peptide, while we failed to identify any known homologous peptide of muscin. This kind of newly found antimicrobial peptides could be a part of the explanation to how the house fly was able to flourish in the septic environments. Although RNA-seq technology reduces the need of technical replication within our experiments and the fact that we used the NOISeq method can empirically model the noise in count data, the absence of biological replication within our experiments constrains the possibility of more detailed analyses. We feel safe to draw the draft conclusion above out of the multiple points data, including pre-challenge and post-challenge 6 h, 24 h, and 48 h. However, studies with more intensive time intervals and biological replications are needed to provide deep insight into the immunogenetics of the house fly, which also may contribute to a better understanding of the evolutionary history of innate immunity from insects to vertebrates. Supporting Information Figure S1 Length distribution of unigenes. (DOC) [146]Click here for additional data file.^ (24KB, doc) Figure S2 The nucleotide and deduced amino acid sequences of M. domestica antimicrobial peptide muscin . (DOC) [147]Click here for additional data file.^ (25KB, doc) Figure S3 The nucleotide and deduced amino acid sequences of M. domestica antimicrobial peptide domsticin . (DOC) [148]Click here for additional data file.^ (26KB, doc) Figure S4 Multiple alignment of domesticin with its homologs from Drosophila by CLC Sequence Viewer. (PDF) [149]Click here for additional data file.^ (260.2KB, pdf) Table S1 Primers used for qPCR analysis. (DOC) [150]Click here for additional data file.^ (47.5KB, doc) Table S2 Top hits obtained by BLASTx for the unigenes. (XLS) [151]Click here for additional data file.^ (25.5MB, xls) Table S3 KO annotation of unigenes. (XLS) [152]Click here for additional data file.^ (42KB, xls) Table S4 The data of all the differentially expressed genes. (XLS) [153]Click here for additional data file.^ (2.6MB, xls) Table S5 KEGG enrichment analysis of all the differentially expressed genes. (XLSX) [154]Click here for additional data file.^ (60.7KB, xlsx) Table S6 Detail information of 509 putative immune-related differentially expressed genes. (XLS) [155]Click here for additional data file.^ (149KB, xls) Acknowledgments