Blunt Snout Bream (Megalobrama amblycephala) MyD88 and TRAF6: Characterisation, Comparative Homology Modelling and Expression

MyD88 and TRAF6 play an essential role in the innate immune response in most animals. This study reports the full-length MaMyD88 and MaTRAF6 genes identified from the blunt snout bream (Megalobrama amblycephala) transcriptome profile. MaMyD88 is 2501 base pairs (bp) long, encoding a putative protein of 284 amino acids (aa), including the N-terminal DEATH domain of 78 aa and the C-terminal TIR domain of 138 aa. MaTRAF6 is 2252 bp long, encoding a putative protein of 542 aa, including the N-terminal low-complexity region, RING domain (40 aa), a coiled-coil region (64 aa) and C-terminal MATH domain (147 aa). Coding regions of MaMyD88 and MaTRAF6 genomic sequences consisted of five and six exons, respectively. Physicochemical and functional characteristics of the proteins were analysed. Alpha helices were dominant in the secondary structure of the proteins. Homology models of the MaMyD88 and MaTRAF6 domains were constructed applying the comparative modelling method. RT-qPCR was used to analyse the expression of MaMyD88 and MaTRAF6 mRNA transcripts in response to Aeromonas hydrophila challenge. Both genes were highly upregulated in the liver, spleen and kidney during the first 24 h after the challenge. While MyD88 and TRAF6 have been reported in various aquatic species, this is the first report and characterisation of these genes in blunt snout bream. This research also provides evidence of the important roles of these two genes in the blunt snout bream innate immune system.


Introduction
The immune system protects an organism against diseases by identifying and eliminating the pathogen [1]. Innate immunity, which is an evolutionarily ancient system present in both invertebrates and vertebrates, is the first line of defence against infection and an internal stimulator for development of antigen-specific acquired immune response and homeostasis [1][2][3][4][5]. The innate immune response is triggered by germline-encoded pattern recognition receptors (PRRs) that are responsible for recognising conserved pathogen-associated molecular patterns of foreign stimuli, such as lipopolysaccharide or peptidoglycan of bacterial cell walls, β-1,3-glucan of fungal cell walls, double-stranded RNA of viruses and endogenous molecules released from damaged cells [5][6][7][8]. To date, four different classes of PRR families have been identified: toll-like receptors (TLRs), C-type lectin receptors (both transmembrane proteins), retinoic acid-inducible gene (RIG)-I-like receptors and NOD-like receptors (both cytoplasmic proteins) [8]. TLRs are responsible for detecting extracellular invading pathogens and intracellular endosomes and lysosomes [9], leading to the activation and recruitment of immune effectors and subsequent stimulation of antimicrobial responses [10]. TLRs are characterised by the N-terminal leucine-rich repeats and a transmembrane region followed by a cytoplasmic Toll/IL-1R homology (TIR) domain [8]. The TIR region is required for initiating downstream signalling by recruiting adaptors that contain a TIR domain: MyD88 (originally identified as myeloid differentiation primary response gene [11]), MyD88-adaptor-like, TIR domain-containing adaptor protein (inducing interferon-β (IFN-β)), TIR-domain-containing adaptor protein and TRIF-related adaptor molecule [12]. While TLR3 and TLR4 are unique in their ability to activate the induction of type I IFNs in MyD88-independent fashion [10,13], most of the TLRs seem to be absolutely dependent on MyD88 for all of their functions [5]. MyD88 stimulates the expression of pro-inflammatory genes, like tumour necrosis factor (TNF) and interleukin-1 (IL-1) through the activation of NF-κB [7,14]. MyD88 consists of a C-terminal TIR and an N-terminal DEATH domain [5,15]. The DEATH domain, a protein interaction module composed of a bundle of six alpha helices, interacts with the corresponding domain in IL-1R-associated kinases (IRAKs), promoting recruitment of downstream immune molecules, including tumour necrosis factor receptor-associated factor 6 (TRAF6) [16]. Activation of TRAF6 by IRAKs recruits the induction of downstream immune molecules that subsequently lead to activation of NF-κB, IRFs and induction of pro-inflammatory cytokines or antiviral genes [8]. TRAF6, the most evolutionarily ancient TRAF family member, is also the only TRAF participating in signal transduction of both the TNF receptor superfamily and the IL-1R/TLR superfamily, which play essential roles in innate immunity, adaptive immunity and bone homeostasis [17].
Both genes have been identified and their functions studied in a broad range of aquatic animals: such as Zhikong scallop (Chlamys farreri) [18], common carp (Cyprinus carpio) [19], miiuy croaker (Miichthys miiuy) [20], whiteleg shrimp (Litopenaeus vannamei) (MyD88) [21], penaeid shrimp (Fenneropenaeus chinensis) [22], zebrafish (Danio rerio) [23], orange-spotted grouper (Epinephelus coioides) [24], whiteleg shrimp (TRAF6) [25] and grass carp (Ctenopharyngodon idella) [26]. However, so far MyD88 and TRAF6 have not been reported in blunt snout bream (Megalobrama amblycephala). Based on the available data from other (related) species, we hypothesize that both genes are present in blunt snout bream and that their structure and functions are conserved. Herein, the transcripts of MaMyD88 and MaTRAF6 were identified from the transcriptomic profile of blunt snout bream and their expression levels after challenge with Aeromonas hydrophila were investigated. Additionally, the entire coding regions of MaMyD88 and MaTRAF6 genes were amplified and sequenced from the genomic DNA. Computational tools were used to analyse the physicochemical characteristics and predict the structure of the encoded proteins. This study provides an insight into the role of MaMyD88 and MaTRAF6 in antibacterial response mechanisms of the blunt snout bream immune system, which can facilitate the utilisation of molecular tools for selective breeding of disease-resistant blunt snout bream broodstock in order to reduce mortality and increase productivity during the cultivation.

Genomic Organisation Analysis
The entire coding regions of the MaMyD88 and MaTRAF6 genes were amplified and sequenced from the genomic DNA. The analysis revealed that both are composed of a chain of exons, separated by introns.

Physicochemical and Functional Characterisation
Mature MaMyD88 and MaTRAF6 protein sequences without the signal peptide were used as the templates for physicochemical characterisation analyses ( Table 1). The theoretical isoelectric point (pI) of both proteins was lower than 7, indicating they are acidic. As proteins carry a net positive charge below, and negative charge above their pI, this information can be used for the purification of proteins on a polyacrylamide gel by isoelectric focusing. The extinction coefficient (EC) of MaMyD88 and MaTRAF6 measured at 280 nm was 41,660 and 30,880 M −1 ·cm −1 (assuming all pairs of cysteine residues form cysteines) and 40,910 and 28,880 M −1 cm −1 (assuming all cysteine residues are reduced), respectively. The estimated half-life of both proteins is 30 h in mammalian reticulocytes in vitro, >20 h in yeast in vivo and >10 h in Escherichia coli in vivo. Based on the instability index (II) value, which is a measure to evaluate the stability of proteins in a test tube, MaMyD88 protein (II = 33.5) is probably stable (II < 40), while MaTRAF6 (II = 58.6) is probably not stable (II > 40) [27]. High AI (aliphatic index) value of MaMyD88 (87.18) and MaTRAF6 (72.14) indicate high thermostability of these proteins [28]. Low grand average hydropathicity (GRAVY) values of both proteins imply they are hydrophilic in natural conditions. Table 1. Physicochemical characteristics of MaMyD88 and MaTRAF6 proteins: the number of amino acids (No. of aa), molecular weight in Da (Mol. Wt.), isoelectric point (pI), total number of negative (−R) and positive residues (+R), extinction coefficient (EC), instability index (II), aliphatic index (AI) and grand average hydropathicity (GRAVY). EC *-the first value is based on the assumption that all pairs of cysteine residues form cysteines and the second one that all cysteine residues are reduced.
MaMyD88 protein is rich in leucine, aspartic acid, lysine and threonine, while MaTRAF6 is rich in leucine, glutamic acid, and serine. Both proteins are classified as soluble. Twelve cysteine residues were found in MaMyD88 and 32 in MaTRAF6. The most probable pattern of cysteine residue pairing predicted in MaMyD88 is Cys105-Cys221 and Cys125-Cys268, while no cysteine pairing was predicted in MaTRAF6.

Protein Structure Prediction and Model Validation
In both proteins, alpha helix was predominant among the secondary structure elements, followed by random coil, extended strand and beta turn. The rest of the structure elements were not predicted ( Table 2).  Ramachandran plot analysis results reveal that the models for TIR domain of MaMyD88 and MATH domain of MaTRAF6 have over 90% (92.2% and 90.8%, respectively) of residues in the most favoured region, suggesting a good quality of the homology models (Table 3). Both remaining models had almost all of the residues in most favoured and additional allowed regions combined (MaMyD88 DEATH-96.7%, MaTRAF6 RING-99.3%) indicating acceptable quality. The overall average G-factor of dihedral angles and main-chain covalent forces for the models ranged from −0.15 to 0.18, indicating a very good quality (>−0.5) of the proposed models [32].
LGscore values also indicate "very good" quality (>5.0), except for the MaTRAF6 RING model (2.436), which is "correct" (>1.5). MaxSub validation measure indicates "very good" quality (>0.8) of the MaMyD88 TIR domain, "correct" quality (>0.1) of MaTRAF6 RING domain and "good" quality of the remaining two models (>0.5) [33]. The Z-Scores of the predicted models range from −7.63 to −5.24, which is within the range of the scores typically found for native proteins of the similar size, while plots of single residue energies revealed predominantly negative values, also indicating good quality of the proposed models [34] (data not shown). All of the validation results suggest that the proposed models of MaMyD88 and MaTRAF6 domains can be accepted as relatively accurate.

Expression of MaMyD88 and MaTRAF6 Transcripts after A. hydrophila Infection
RT-qPCR was used to detect the changes in expression of MaMyD88 and MaTRAF6 in liver, spleen and kidney of blunt snout bream after A. hydrophila infection at different time points (4, 12, 24, 72, and 120 h post infection-hpi). In comparison to the control group (normalised to the internal control gene 18S rRNA expression levels), MaMyD88 gene was significantly (p < 0.01) overexpressed in liver (38.6-fold) at 4 hpi, the expression was then reduced at 12 hpi (<5-fold), just to rise again significantly (>30-fold, p < 0.01) at 24 hpi. Expression in kidney followed a similar pattern: 15.9-fold (p < 0.01) at 4, <5-fold at 12 and >10-fold (p < 0.01) at 24 hpi. Expression in spleen was significantly (p < 0.01) upregulated at all three time-points, but the pattern was different: less than 5-fold at 4 hpi, 20.3-fold at 12 and then reverting to <5-fold level at 24 hpi. Expression at 72 and 120 hpi in all tissues was close to, or even below, the normal (control group) levels ( Figure 8A).
The expression of MaTRAF6 in liver followed a pattern somewhat similar to MaMyD88, albeit the over-expression levels were not nearly as high: significantly (p < 0.01) over-expressed (6.4-fold) at 4 hpi, significantly (p < 0.05) underexpressed at 12 hpi and then again over-expressed (2.8-fold) at 24 hpi. Expression then tailed off to the normal level at 72 hpi and to downregulated at 120 hpi (0.2-fold). Expression in kidney and spleen roughly followed a similar, bell curve pattern: peaking at 12 (3.9-fold, p < 0.01) and 24 (5.4-fold, p < 0.05) hpi, respectively, and then tailing off towards the normal level ( Figure 8A). . Each histogram represents the mean ± SE of three replicates. Statistically significant differences from the control group are marked as * p < 0.05 and ** p < 0.01.

Discussion
In this study, MyD88 and TRAF6 genes were identified and characterised in blunt snout bream for the first time. The genes were identified both from the transcriptome profile and genomic DNA. To test whether the functions of these genes are conserved in blunt snout bream, phylogenetic analysis was conducted, functional domains of amino acid sequences were predicted and characterised, secondary and tertiary structure of the proteins were predicted using computational tools and mRNA expression patterns in response to pathogen challenge were analysed. Genomic organisation analysis may also be helpful for further studies associated with identifying and developing gene polymorphisms related to disease susceptibility/resistance in blunt snout bream. Phylogenetic analyses suggest that both proteins are very similar to homologs from other cyprinids. There is even more divergence between different common carp MyD88 sequences, than between the common carp ADC45019 MyD88 and MaMyD88 ( Figure 5). While no differences were found between cDNA and genomic MaTRAF6 sequences, genomic MaMyD88 translates into a larger putative protein than cDNA due to insertion at the position 204 ( 204 QVCL 207 ). However, alignment with other sequences revealed the absence of this insertion (Figure 2), suggesting that the cDNA sequence is the correct protein template.
The presence of conserved functional domains in MaMyD88 protein suggests that it probably has the same function (interaction with TLRs and IL-1R) as its homologs in other animals [26,28]. The DEATH domain plays important functions in DEATH signal transduction, regulation of apoptosis and inflammatory responses [35], whereas the TIR domain is important in activating the innate immune response by the TLR/IL-1R superfamily mediated pathway [36]. Three highly conserved box regions found within the TIR domain ( Figure 1A) are in accordance with results reported in frogs, mammals and fish [28]: boxes 1 and 2 mediate the coupling to inflammation signalling pathways and box 3 is primarily related to the control of subcellular location of the receptor, possibly through interactions with cytoskeletal elements [37]. Furthermore, the presence of mRNA instability motif in the 3' UTR is characteristic of some fish inflammatory mediator-coding genes and is believed to be responsible for destabilising mRNA [27,38].
MaTRAF6 primary and secondary structure is comparable to those reported in other fish species, such as zebrafish [23], common carp [19] and grass carp [26]. However, zinc fingers were seemingly absent from the MaTRAF6 amino acid sequence. The RING and/or zinc finger regions of TRAFs are critical for signalling, probably due to their association with downstream kinase molecules such as MEKK1 [31]. When it was compared to the most similar, grass carp homolog (Figure 6), where two zinc fingers were reported: 152-204 and 205-262 aa [26], two aa mutations were found in the first zinc finger. At the position 195-glycine (G/GGG) in grass carp was changed to glutamic acid (E/GAA) in MaTRAF6 and at the position 201-lysine (K/AAG) was changed to threonine (T/ACC). However, the second zinc finger appears identical in both species, suggesting that the reported differences could have actually been caused by different parameters used by the employed software, and not by different functional properties of the proteins. To further corroborate this, a common carp TRAF6 sequence, where no zinc fingers were reported [19], was also compared to the grass carp: only one mutation was found in the putative zinc finger one region (position 201) and none in the zinc finger two. Comparative analysis revealed that position 201 is the least conserved position in the cyprinid zinc finger one domain, suggesting its low relevance for the protein function. Similarly, only one zinc finger was reported in zebrafish, spanning aa 205 to 260 [23]. Only one mutation was found at the position 256, where threonine (T) was changed to proline (P/CCG) in MaTRAF. Furthermore, MaTRAF6 RING and MATH domains share very high similarity with the homologous domains in other teleosts (Figure 4). The coiled-coil region plays a key position in the auto-ubiquitination and activation of the NF-κB signalling pathway [39], while the MATH domain is responsible for binding TRAF6 to upstream molecules [17]. While seven exons and six introns were found in the genomic sequences of TRAF6 in orange-spotted grouper [24] and zebrafish [23], identical numbers of exons were reported in common carp [19] and blunt snout bream MyD88 (5) and TRAF6 (6) genes. All these results suggest highly conserved functions of the homologs.
As protein structure often reflects its functions [40], physicochemical analysis and prediction of protein structures can provide basic insights into the functions of proteins at the molecular level. In theory, secondary protein structure is directly related to the hydrogen bonding of α-helix and β-sheet, while random coils usually indicate the absence of a regular secondary structure [41]. The observed patterns of cysteine pairing in MaMyD88 indicate that the protein contains disulphide bonds that are essential in the folding of the protein and responsible for stabilisation of protein structure [42]. Analysis of three-dimensional structure of proteins has previously revealed the structural differences between receptor recognition by TRAF6 and other TRAFs, providing insight into the mechanisms by which TRAF6 mediates several signalling cascades [17]. As there is no available data about tertiary structure of these proteins in blunt snout bream, or even in closely related fish species, this study provides the basis for the further research of functional properties of these genes in fish.
Regulation patterns of both MaMyD88 and MaTRAF6 in response to a challenge by A. hydrophila reported here suggest a response of the host's immune system to bacterial infection, and thus are consistent with the hypothesis of their conserved function in the immune system. Similar patterns were also observed in other studied fish: miiuy croaker MyD88 was significantly upregulated at different time points (from 12 to 72 hpi) in liver, kidney and spleen after a challenge with Vibrio anguillarum [20]; Yan et al. [28] reported that MyD88 transcription was induced in orange-spotted grouper spleens eight hours after a challenge with Singapore grouper iridovirus; it was also highly expressed in all the examined tissues in large yellow croaker (Pseudosciaena crocea) in response to formalin-inactivated Vibrio parahaemolyticus challenge [43].
In agreement with results presented here, TRAF6 expression in aquatic animals was also previously found to be highly responsive to pathogen challenges [16,[24][25][26]. In grass carp, TRAF6 was up-regulated in spleen and head kidney six hours after a challenge with Ichthyophthirius multifiliis, while the expression was significantly down-regulated in head kidney at 72 hpi [26]. Similar expression patterns were observed in some crustaceans as well: in whiteleg shrimp, TRAF6 was upregulated in gills and hepatopancreas at 3 hpi with white spot syndrome virus (WSSV), and then down-regulated and maintained at a low level until the last tested time-point (24 hpi) [25]; similarly, TRAF6 was upregulated in tiger shrimp (Penaeus monodon) haemocytes and lymphoid organ at the late stages of WSSV and poly I:C infection [16]. These results strongly suggest that MaMyD88 and MaTRAF6 play key roles in the blunt snout bream innate immunity. A good understanding of the blunt snout bream immune system is necessary for the prevention of diseases relevant for the farming of this species, however, further studies of the specific functions of these two genes are needed in order to fully confirm these assumptions.

Fish and Challenge Experiment
Healthy blunt snout bream specimens (average body weight: 27.3 ± 6.8 g) used in this study were collected from the farm in Tuanfeng, Huanggang, Hubei, China. Fish were maintained in a 1 m 3 tank with aeration at about 28 °C for 2 weeks in the laboratory at the College of Fisheries, Huazhong Agricultural University. Fish were fed commercial pelleted feed twice a day. For the analysis of mRNA expression, fish were challenged with A. hydrophila. Fish were divided into two groups: 87 specimens in the control group and 151 specimens in the experimental group. Fish from both groups were intraperitoneally injected equal volumes (0.1 mL) of sterile physiological saline (control) and a bacterial suspension at 1.8 × 10 6 cfu/mL (experimental). Six specimens from each group were sampled at 4, 12, 24, 72 and 120 hpi anesthetised in MS-222 (Sigma-Aldrich, Saint Louis, MO, USA) at 100 mg/L concentration and immediately killed. Liver, spleen and kidney were separately sampled from all specimens, immediately flash-frozen in liquid nitrogen and stored at −80 °C.

Total RNA Preparation and cDNA Synthesis
The total RNA was extracted from each sample with RNAisoPlus Reagent (Takara Bio Inc., Dalian, China), according to the manufacturer's instructions. Quality and quantity of the extracted RNA were estimated using electrophoresis in 1% agarose gel and Nanodrop 2000 spectrophotometry (Thermo Scientific, Wilmington, DE, USA). Equal amounts of the total RNA of six specimens sampled at each time-point were pooled. cDNA library was synthesised using PrimeScript ® RT reagent Kit with gDNA Eraser (Takara Bio Inc., Dalian, China) following the manufacturer's instructions, serially diluted 10-fold and used as the template for RT-qPCR.

Cloning of Full-Length cDNAs and Bioinformatics Analyses
The full-length cDNAs of MaMyD88 and MaTRAF6 were obtained from the blunt snout bream transcriptome profile, constructed using Solexa/Illumina technology. The unigenes were firstly obtained via the de novo assembling of short reads from the transcriptome of blood, liver, gill, intestine, spleen and kidney. Then, the transcripts were annotated and identified through BLAST homology search against the GenBank database [44]. The full-length cDNA sequences of MaMyD88 and MaTRAF6 were obtained using SMART™ RACE cDNA Amplification Kit (Takara Bio Inc., Dalian, China), following the manufacturer's instructions. The specific primers are described in Table 4. Thermal conditions for RACE-PCR were as follows: 94 °C for 5 min; 30 cycles of 94 °C for 30 s, 65 °C for 30 s, 72 °C for 2 min; extension at 72 °C for 5 min. The 3'-RACE and 5'-RACE products were ligated into the pGEM-T Easy vector (Promega, Madison, WI, USA) for sequencing. The full-length cDNA sequences were assembled using the SeqMan software. The putative amino acid sequences were predicted using the NCBI's Open Reading Frame Finder [45]; signal peptide was predicted by using the SignalP server [46]; the poly-A tail was predicted by using The GENSCAN Web Server at MIT [47] and the domain structures using SMART program [48]. Homologous MyD88 and TRAF6 amino acid sequences from other species were obtained from the NCBI database. Multiple-sequence alignment was performed using the ClustalW2 server [49]. Neighbour-joining phylogenetic trees were constructed using MEGA 5.2 [50], and the topological stability of the trees was evaluated by 1000 bootstrap replications.

Protein Physicochemical and Functional Characterisation
Expasy's ProtParam prediction server [51] was employed to determine the physicochemical properties of MaMyD88 (284 aa) and MaTRAF6 (542 aa) polypeptide chains: molecular weight, amino acid composition, theoretical isoelectric point, total number of positive and negative residues, extinction coefficient, instability index, aliphatic index and grand average of hydropathicity. SOSUI [52] was used to identify the types of protein (soluble or membrane) and CYS_REC [53] to predict the presence of disulphide bonds and their bonding patterns.

Protein Structure Prediction
Secondary structure of MaMyD88 and MaTRAF6 was predicted using the SOPMA server [54], with default parameters (window width-17; similarity threshold-8; number of states-4), while three-dimensional homology models were constructed using the SWISS-MODEL server [55]. The modelled structures were selected on the basis of sequence identity with the Protein Data Bank (PDB) templates [56]. Stereochemical quality and accuracy of the predicted models were analysed using PROCHECK's Ramachandran plot analysis [57,58], ERRAT [59], ProQ [34] and ProSA [35].

Quantitative Real-Time PCR (RT-qPCR) and Statistics
RT-qPCR was performed with SYBR ® Premix Ex Taq™ (Takara Bio Inc., Dalian, China) in a Rotor-Gene Q real-time PCR cycler (Qiagen, Dusseldorf, Germany). The total reaction volume of 10 mL contained 5 µL of SYBR ® Premix Ex Taq II (2×), 0.4 µL of each primer (10 µM), 0.8 µL of cDNA template and 3.4 µL of dH2O. The primers (Table 4) were designed using Primer Premier 5 software (Premier Biosoft, Palo Alto, CA, USA) and synthesised by Sangon Biotech Co., Ltd. (Shanghai, China). 18S rRNA was used as the reference gene [60]. Thermal conditions were as follows: denaturation at 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s, annealing temperature at 55 °C-MaMyD88/53 °C-MaTRAF6/60 °C-18S rRNA for 20 s and elongation at 72 °C for 15 s. All reactions were performed in triplicate. The data were analysed using the Rotor-Gene Q series software 1.7 (build 94) (Qiagen). The Ct values were obtained using a threshold value of 0.05 for all three genes. Relative gene expression was quantified by the comparative Ct method, expressed as 2 −ΔΔCt [61]. The relative expression of the target genes was normalised to an endogenous reference (18S rRNA), where ΔCt was calculated as CtTest − Ct18S rRNA, and ΔΔCt was calculated as ΔCtTest − ΔCtControl. The RT-qPCR data were analysed statistically with Microsoft Excel and with one-way analysis of variance (ANOVA) in the SPSS 16.0 software. Differences were considered statistically significant at p < 0.05 and p < 0.01.