Abstract
IL-1R–associated kinases (IRAK) are important regulators in the TLR/IL-1R pathways, but their function appears inconsistent between Drosophila, bony fishes, and vertebrates. This causes a difficulty to understand the IRAK functions. As a step to reveal the evolution of IRAKs, in this study, we performed comparative and functional analysis of IRAKs by exploiting the amphioxus, a pivotal taxon connecting invertebrates and vertebrates. Sequence and phylogenetic analysis indicated three major IRAK lineages: IRAK1/2/3 is a vertebrate-specific lineage, IRAK4 is an ancient lineage conserved between invertebrate and vertebrates, and Pelle is another ancient lineage that is preserved in protostomes and invertebrate deuterostomes but lost in vertebrate deuterostomes. Pelle is closer neither to IRAK4 nor to IRAK1/2/3, hence suggesting no clear functional analogs to IRAK1/2/3 in nonvertebrates. Functional analysis showed that both amphioxus IRAK4 and Pelle could suppress NF-κB activation induced by MyD88 and TRAF6, which are unlike mammalian and Drosophila IRAKs, but, surprisingly, similar to bony fish IRAK4. Also unlike Drosophila IRAKs, no interaction was detected between amphioxus IRAK4 and Pelle, although both of them were shown capable of binding MyD88. These findings, together with previous reports, show that unlike other signal transducers in the TLR/IL-1R pathways, such as MyD88 and TRAF6, the functions of IRAKs are highly variable during evolution and very specialized in different major animal taxa. Indeed, we suggest that the functional variability of IRAKs might confer plasticity to the signal transduction of the TLR/IL-1R pathways, which in return helps the species to evolve against the pathogens.
Introduction
Toll-like receptors are pattern recognition receptors, whereas IL-1R are cytokine receptors. They both share a TLR/IL-1R (TIR) domain at their cytoplasmic tail, and both have important roles in the host-against-pathogen immune responses (1, 2). Both types of receptors receive extracellular signals and induce the classic intracellular MyD88–TRAF6–NF-κB pathway (2, 3). Dysregulation of this pathway may cause chronic inflammation and autoimmune diseases, etc (4, 5).
This pathway is regulated by an important class of kinases, the IL-1R–associated kinases (IRAKs) (6). In mammals, there are four IRAK genes, IRAK1, IRAK2, IRAK3 (or IRAKM), and IRAK4. All four IRAKs have similar protein architecture: an N-terminal death domain, a ProST domain, a kinase domain, and a C-terminal domain (except IRAK4) (6). The death domain is used for interaction with other signal transducers such as MyD88 and other IRAK members (7). The ProST domain is rich in serines, prolines, and threonines. Phosphorylation in this domain plays a role in regulating IRAKs’ own availability (8). IRAK1 is reported to undergo hyperphosphorylation in ProST domain, leading to its release from MyD88 complex and bind to TRAF6 (9). The kinase domain apparently confers kinase activity, which, however, may not be required for some function of IRAKs. For example, IRAK3 lacks kinase activity, and the kinase inactive mutants of IRAK1 and 4 are still able to activate NF-κB (10–12). The C-terminal domain of IRAK1/2/3 is reportedly important for the interaction with TRAF6 (carrying TRAF6 binding sites) (6, 13).
In a typical signal transduction in mammals, TLRs/IL-1Rs use intracellular TIR domain to recruit MyD88 (2, 14, 15). Then, MyD88 recruits IRAKs via the death domains to form a signaling complex (Fig. 1A) (7, 16). IRAK4 is suggested to be the first one to be recruited by MyD88, and its autophosphorylation after recruitment leads to the subsequent recruitment and phosphorylation of IRAK1 (10, 16–19). Hyperphosphorylated IRAK1 is then released from the MyD88 complex and interacts with TRAF6, which in turn activates NF-κB (20, 21). IRAK2 might have a redundant function with IRAK1, but new evidence suggests a different role for IRAK2 in the TLR/IL-1R signaling (6, 10, 22). In contrast, IRAK3 lacks kinase activity but is capable of suppressing the pathway by preventing IRAK1/2 from the dissociation from MyD88 and from the interaction with TRAF6 (23, 24).
However, there are few studies on IRAKs from lower vertebrates and invertebrates. But it appears that the composition and function of the IRAK family could be quite different from mammals (25–32). For example, in lower vertebrates, birds lack IRAK1 and IRAK3, reptiles and amphibians lack IRAK3, and bony fishes lack IRAK2 (33). Interestingly, several reports suggested that bony fish IRAK4 could suppress MyD88-dependent NF-κB activation when transfected into mammalian cells (28, 34, 35), which is in contrast to the activating role of mammalian IRAK4 (Fig. 1A).
As about nonvertebrates, whereas arthropods such as Drosophila have two IRAK homologs (Tube/IRAK4 and Pelle), it appears that poriferans and hemichordates have only IRAK4 (29, 33, 36, 37). It is suggested that the IRAK family derives from a single IRAK4-like kinase in the last common metazoan ancestor (33, 38). Functions of invertebrate Tube/IRAK4 and Pelle have been studied in Drosophila melanogaster. When activated, D. melanogaster Tube and D. melanogaster Pelle are recruited to D. melanogaster MyD88 to form a trimeric complex via the death domains, which is essential for subsequent signal transduction (39–41). Although D. melanogaster Pelle and D. melanogaster MyD88 do not directly contact with each other, they bind to distinct surface of the death domain of Tube (39, 42). Pelle is autophosphorylated and catalyzes the phosphorylation of Cactus, leading to the degradation of Cactus and the activation of Dorsal or Dif (Fig. 1A) (43, 44).
That being said, the function of IRAKs in lower vertebrates and invertebrates is far from being revealed. To gain more understanding of the functional evolution of IRAKs, in this study we exploited an evolutionary model organism called amphioxus and performed corresponding comparative and functional analysis. Cephalochordate amphioxus is the most basal living lineage among three chordate subphyla (cephalochordates, urochordates, and vertebrates), representing a pivotal transitional taxon connecting invertebrates and vertebrates (45–50). There have been extensive studies on the molecular functions of the TLR/IL-1R pathways in amphioxus (51–65), which provide a robust premise for the function study of amphioxus IRAKs. Our study shows that the functions of IRAKs are highly variable between different major animal taxa, suggesting that IRAKs might confer necessary plasticity to the signal transduction of the TLR/IL-1R pathways.
Materials and Methods
Animals and cells
Adult Branchiostoma belcheri were obtained from local fishermen of Zhanjiang, China, and cultured in aquaria supplied with circulating filtered seawater and fed once every 2 d. For immune stimulation, each amphioxus were injected 15 μl LPS (1 mg/ml) in PBS into the coelom, and PBS was injected as negative control. Then, intestine and gill from five fishes were collected and frozen using liquid nitrogen immediately at 0, 2, 8, 24, and 48 h after injection.
HEK293T and Hela cells were grown in DMEM supplemented with 10% FCS at 37°C.
RNA isolation and cDNA synthesis
The total RNA was isolated using TRIzol Reagent (Invitrogen) and isopropanol was used for precipitation. After quality inspection with agarose gel electrophoresis and spectrophotometrically, the extracted RNA was reverse transcribed to synthesize the first strand cDNA using PrimeScript first Strand cDNA Synthesis Kit (Takara). Then, the cDNA was stored at −80°C.
Cloning of B. belcheri IRAK4 and B. belcheri Pelle
Using human IRAKs as baits, we blasted the genome database of B. belcheri (http://genome.bucm.edu.cn/lancelet/index.php) and obtained the consensus sequences of B. belcheri IRAK4 (BbIRAK4) and B. belcheri Pelle (BbPelle). Then two pair gene-specific primers were designed for each gene and nested PCR were performed to amplify the coding sequence of BbIRAK4 and BbPelle. Amplified sequences were inserted into the pGEX-T easy vector (Promega) and verified by sequencing. The primers were showed in Table I.
Bioinformatic analysis
The domain structure was predicted on the main page of the SMART Web site (http://smart.embl-heidelberg.de). BLASTP was performed to analyze the sequence identities. The isoelectric point (pI) and m.w. were speculated on ExPASy Web site (http://www.expasy.org/tools/). Multiple sequence alignments were analyzed using ClustalX 1.83 and were manually corrected using GeneDoc. Neighbor-joining tree was built using MEGA5 with 1000 bootstrap tests and handling gaps by pairwise deletion.
Quantitative real-time PCR
Quantitative real-time PCR (qRT-PCR) was performed to determine the mRNA expression profiles of BbIRAK4 and BbPelle in amphioxus. Healthy adult individuals were chosen for tissue distribution analysis. Tissues (muscle, skin, gill, ovary, hepatic cecum, intestine) were harvested, and total RNA was extracted. qRT-PCR was conducted on a LightCycler 480 instrument using SYBR PrimeScript qRT-PCR kit (Takara). Data were quantified using the 2−ΔΔCt method based on the cycle threshold values of BbIRAK4, BbPelle, and GAPDH. All of the samples were analyzed in three replicates. Primers used for qRT-PCR were listed in Table I.
Construction of the expression vector
For the expression of BbIRAK4 and BbPelle in Hela and HEK293T cells, PCR fragments encoding full-length protein sequences of BbIRAK4 and BbPelle were fused with the cutting site of EcoR I and Xho I restriction enzyme. The PCR fragments and pCMV-hemagglutinin vector (Clontech) were both digested with EcoR I and Xho I restriction enzyme at 37°C for overnight. Then, DNA T4 ligase (Takara) was used to ligate the fragments and vector. The recombinant expression vectors were verified by sequencing and the expression of proteins were confirmed by Western blot. Primers used for expression vector construction are showed in Table I.
Immunofluorescence analysis
Hela cells were digested by trypsin and added to a 24-well plate covered with pretreated coverslips (10 × 10 mm). After 16–24 h, Hela cells were transfected with 500 ng BbIRAK4 or BbPelle recombinant expression vectors by PEI (Polysciences) according to the manufacturer’s instructions. After 24 h, Hela cells were harvested and fixed in 4% formaldehyde solution for 15 min. Then cells were permeabilized by washing three times in 0.05% Tween–20 in PBS (PBST) and blocked nonspecific sites with 5% BSA in PBST at room temperature for 1 h. Next, cells were incubated with anti-hemagglutinin Ab for 1 h and washed three times in PBST at room temperature. Hela cells were further incubated with Alex Fluor 594 Goat Anti-Mouse IgG (H + L) (Proteintech) at 4°C overnight. After triple washing in PBS, cells were labeled with DAPI for 5 min and washed three times in PBS. The coverslips with treated Hela cells were taken out from the 24-well plate and cover on the slides with MOWIOL R4-88 Reagent (Calbiochem). Finally, Hela cells were photographed by Carl Zeiss Axio Imager Z1 microscope.
Luciferase reporter assays
HEK293T cells were digested by trypsin and seeded in 48-well plates. After 16–24 h, cells were transfected with equivalent mixed expression plasmids, which consist of the indicated amount of expression vectors, 50 ng per well of the NF-κB response promoter luciferase reporter plasmid pNF-κB-Luc (StrataGene), and 5 ng per well Renilla luciferase reporter plasmid pRL-TK (Promega) to normalize the data due to different transfection efficiency between wells, and empty vectors to complement the total plasmid quantity of each well to the same. After 24 h, HEK293T cells were harvested and measured by dual-luciferase reporter assay system (Promega). Each experiment was performed at least in triplicate and was repeated at least twice in all cases. The data are shown as the fold change relative to that measured in cells transfected with empty vector.
Coimmunoprecipitation assay
HEK293T cells were digested by trypsin and seeded in six-well plates. After 16–24 h, cells were transfected with 4–6 μg expression plasmids. At 20–36 h posttransfection, whole-cell extracts were prepared by using cell lysis buffer (Cell Signaling Technology) and incubated with primary Abs (Sigma-Aldrich) at 4°C overnight. The next day, Protein G Sepharose (Roche) was added to the mixture and incubated at 4°C for 4–6 h. Then the mixture was washed three times with cell lysis buffer and analyzed by Western blot.
Results
Identification and cloning of two amphioxus IRAK genes
We performed an extensive search in the genome and transcriptome databases and identified two amphioxus IRAK genes. So now it is clear that the genomes of amphioxus, including B. belcheri and B. floridae, only encode two IRAK genes. In this study, we designated them as IRAK4 and Pelle, respectively, according to the results of sequence, phylogenetic and function analyses below. Then we cloned the full-length cDNA of BbIRAK4 and BbPelle with nested RT-PCR (Table I).
Primers . | Primer sequence (5ˊ to 3ˊ) . |
---|---|
Gene-specific primers | |
BbIRAK4-F | 5′-TATTTTTGTGTGGTCTGAGAAC-3′ |
BbIRAK4-R | 5′-TCCCACTCTACTCTCTTAAAATC-3′ |
BbIRAK4-NF | 5′-ATGGGAAGCCTGACAAATATGCCAG-3′ |
BbIRAK4-NR | 5′-TCAAATACAGTCCAAGTTCTTTTTCA-3′ |
BbPelle-F | 5′-TTGGACTACAGGAGACCGGACGGA-3′ |
BbPelle-R | 5′-ATCAAACGCCTTCTGTCCAGTC-3′ |
BbPelle-NF | 5′-ATGGCGGCTGCTGTTCTCGGTAGCA-3′ |
BbPelle-NR | 5′-TCAAACACCACCACGCCAATCAAA-3′ |
Primers for qRT-PCR | |
BbIRAK4-rtF | 5′-ACTGGTGGAGGTCTTACTAC-3′ |
BbIRAK4-rtR | 5′-TGGCTGGTCTGAGATGTT-3′ |
BbPelle-rtF | 5′-CCGAAGTTCTCAAGCCTCTA-3′ |
BbPelle-rtR | 5′-TTGCCGAAGTTCTCATCAG-3′ |
Primers for construction of expression vector | |
BbIRAK4-vF | 5′-CCGGAATTCCGatgggaagcctgacaaatatgccag-3′ |
BbIRAK4-vR | 5′-CCGCTCGAGCGGtcaaatacagtccaagttctttttca-3′ |
BbPelle-vF | 5′-CCGGAATTCCGatggcggctgctgttctcggtagca-3′ |
BbPelle-vR | 5′-CCGCTCGAGCGGtcaaacaccaccacgccaatcaaa-3′ |
Primers . | Primer sequence (5ˊ to 3ˊ) . |
---|---|
Gene-specific primers | |
BbIRAK4-F | 5′-TATTTTTGTGTGGTCTGAGAAC-3′ |
BbIRAK4-R | 5′-TCCCACTCTACTCTCTTAAAATC-3′ |
BbIRAK4-NF | 5′-ATGGGAAGCCTGACAAATATGCCAG-3′ |
BbIRAK4-NR | 5′-TCAAATACAGTCCAAGTTCTTTTTCA-3′ |
BbPelle-F | 5′-TTGGACTACAGGAGACCGGACGGA-3′ |
BbPelle-R | 5′-ATCAAACGCCTTCTGTCCAGTC-3′ |
BbPelle-NF | 5′-ATGGCGGCTGCTGTTCTCGGTAGCA-3′ |
BbPelle-NR | 5′-TCAAACACCACCACGCCAATCAAA-3′ |
Primers for qRT-PCR | |
BbIRAK4-rtF | 5′-ACTGGTGGAGGTCTTACTAC-3′ |
BbIRAK4-rtR | 5′-TGGCTGGTCTGAGATGTT-3′ |
BbPelle-rtF | 5′-CCGAAGTTCTCAAGCCTCTA-3′ |
BbPelle-rtR | 5′-TTGCCGAAGTTCTCATCAG-3′ |
Primers for construction of expression vector | |
BbIRAK4-vF | 5′-CCGGAATTCCGatgggaagcctgacaaatatgccag-3′ |
BbIRAK4-vR | 5′-CCGCTCGAGCGGtcaaatacagtccaagttctttttca-3′ |
BbPelle-vF | 5′-CCGGAATTCCGatggcggctgctgttctcggtagca-3′ |
BbPelle-vR | 5′-CCGCTCGAGCGGtcaaacaccaccacgccaatcaaa-3′ |
F, outer forward; NF, inner forward; NR, inner reverse; R, outer reverse; rtF, qRT-PCR forward; rtR, qRT-PCR reverse; vF, vector forward; vR, vector reverse.
Protein architectural diversity of the IRAK genes
The coding sequence of BbIRAK4 is 1488 bp, encoding 495 aa with a theoretical pI of 5.19 and a calculated molecular mass of 54.6 kDa. It shares ∼84% protein sequence identity with its B. floridae ortholog and ∼40% identity with IRAK4 from vertebrates and other invertebrates. BbIRAK4 contains an N-terminal death domain (residues 1–127), a ProST domain (residues 128–222) and a C-terminal kinase domain (residues 223–495), which is more similar to vertebrate IRAK4 than to Drosophila Tube because Tube lacks the kinase domain (Fig. 1B).
Current knowledge of IRAKs in mammal, bony fish, amphioxus, and Drosophila. (A) The functions of IRAKs in different species. We suggest that BbIRAK4 and BbPelle can repress the activity of BbMyD88 and BbTRAF6, but so far no evidence suggests the direct interaction between BbTRAF6 with BbIRAK4 or BbPelle. (B) The domain architecture of human, zebrafish, amphioxus, and Drosophila IRAKs. Important functional residues were indicated. The C-terminal of Drosophila Tube lack kinase domain, which was indicated as white. In RD kinases, an arginine (R) residue immediately precedes the critical aspartate (D) residue in the kinase activation loop, which is vital for their catalytic activity. And in non-RD kinases, this residue is not arginine (R). The kinase activity of DrIRAK3 is unclear.
Current knowledge of IRAKs in mammal, bony fish, amphioxus, and Drosophila. (A) The functions of IRAKs in different species. We suggest that BbIRAK4 and BbPelle can repress the activity of BbMyD88 and BbTRAF6, but so far no evidence suggests the direct interaction between BbTRAF6 with BbIRAK4 or BbPelle. (B) The domain architecture of human, zebrafish, amphioxus, and Drosophila IRAKs. Important functional residues were indicated. The C-terminal of Drosophila Tube lack kinase domain, which was indicated as white. In RD kinases, an arginine (R) residue immediately precedes the critical aspartate (D) residue in the kinase activation loop, which is vital for their catalytic activity. And in non-RD kinases, this residue is not arginine (R). The kinase activity of DrIRAK3 is unclear.
In contrast, the coding sequence of BbPelle is 1353 bp, encoding 450 aa with predicted molecular mass of 50.5 kDa and theoretical pI of 7.64. It shares ∼86% protein sequence identity with its B. floridae counterpart but shares only ∼35% identity with vertebrate IRAK1/2/3 and other invertebrate Pelle. BbPelle contains a death domain (residues 1–117), a ProST domain (residues 118–221), and a kinase domain (residues 222–450) but lacks the C-terminal region presented in vertebrate IRAK1/2/3 (Fig. 1B). In this sense, BbPelle is more similar to Drosophila Pelle.
Further analysis shows that the kinase domain of amphioxus IRAK4 shares conserved residues with other IRAK4, such as the invariant lysine residues for ATP binding, the threonine/serine residues important for kinase activity, and the tyrosine gatekeeper residue, which are all important for IRAK4 function (Supplemental Fig. 1). As for amphioxus Pelle, its kinase domain is clearly shorter than that of other invertebrate Pelle and vertebrate IRAK1/2/3, which might cause differences in function. Amphioxus Pelle also preserves the invariant lysine residues and the conserved threonine residues, which are critical for the phosphorylation for Drosophila Pelle and vertebrate IRAK1/2 (Supplemental Fig. 1).
Protein kinases can be classified as RD or non-RD kinase based on the residue that immediately precedes the critical aspartate (D) residue in the kinase activation loop region with arginine (R) residue or without (66). Like vertebrate IRAK4, amphioxus IRAK4 belongs to the RD kinase groups, whereas amphioxus Pelle, together with other invertebrate Pelle and vertebrate IRAK1/2/3, are non-RD kinases (Fig. 1B, Supplemental Fig. 1). This indicates that amphioxus IRAK4 could be ortholog of vertebrate IRAK4. Whereas amphioxus Pelle, like Drosophila Pelle (33), could share a certain degree of functional similarity with vertebrate IRAK1/2/3.
Taken together, between human, zebrafish, amphioxus, and Drosophila IRAKs, there are obvious differences in terms of domain (and motif) architectures. We expect that this type of architectural diversity may reflect on their functional diversity.
The phylogenetic origins of the three IRAK lineages
To further understand the sequence evolution of the IRAK family, we collected related sequences from a wide range of species and performed phylogenetic analysis. Two phylogenetic trees were constructed for the kinase domain and the death domain, respectively (Fig. 2). The tree based on kinase domains suggests that IRAK proteins are clearly distinct from other related kinase proteins such as TAK1, TBK1, and IKKs. In contrast the tree based on death domains indicates that the death domain of IRAK proteins is more similar to MyD88 (mainly mediating NF-κB activation) than to FADD and TRADD (mainly mediating apoptosis). Therefore, both trees suggest that the IRAK family is a distinct kinase family.
Phylogenetic trees of IRAKs based on the sequence of kinase domains (A) and death domains (B). The tree was constructed using the neighbor-joining method. Numbers on the lines indicate the percentage bootstrap values for 1000 replicate. Maximum likelihood trees were also constructed with the Posisson model, which were similar to neighbor-joining trees (data not shown). The accession numbers for the used protein sequences are as follows: Homo sapiens IRAK1, NP_001560.2; H. sapiens IRAK2, NP_001561.3; H. sapiens IRAK3, NP_009130.2; H. sapiens IRAK4, NP_001107654.1; H. sapiens TAK1, NP_663304.1; H. sapiens TBK1, NP_037386.1; H. sapiens IKKa, NP_001269.3; H. sapiens IKKb, NP_001547.1; Mus musculus IRAK1, NP_001171444.1; M. musculus IRAK2, NP_751893.3; M. musculus IRAK3, NP_082955.2; M. musculus IRAK4, NP_084202.2; M. musculus TAK1, NP_0033342.1; M. musculus TBK1, NP_062760.3; M. musculus IKKa, NP_031726.2; M. musculus IKKb, O88351.1; DrIRAK1, XP_005166760.1; DrIRAK3, XP_003198307.2; DrIRAK4, NP_956457.1; D. rerio TAK1, NP_001018586.1; D. rerio TBK1, NP_001038213.2; D. rerio IKKa, NP_956611.1; D. rerio IKKb, NP_001116737.1; D. melanogaster Pelle, NP_476971.1; D. melanogaster Tube, NP_001189164.1; D. melanogaster TAK1, NP_524080.1; D. melanogaster IKKb, NP_524751.3; Ciona intestinalis IRAK4, XP_002122012.1; Saccoglossus kowalevskii IRAK4, XP_002739502.1; Caenorhabditis elegans PIK1, NP_001255742.1; Gallus gallus IRAK2, NP_001025776.1; G. gallus IRAK4, NP_001025909.1; Anolis carolinensis IRAK1, XP_008102118.1; A. carolinensis IRAK2, XP_008115205.1; A. carolinensis IRAK4, XP_003221322.2; Xenopus tropicalis IRAK1, NP_001006713.1; X. tropicalis IRAK2, NP_001007501.1; Xenopus tropicalis IRAK4, NP_001116877.1; Eriocheir sinensis Pelle, AKJ26285.1; E. sinensis Tube, AGT21373.1; Scylla paramamosain (SpPelle), AGY49577.1; S. paramamosain (SpTube), AGY49576.1; Camponotus floridanu Pelle, XP_011265220.1; C. floridanu Tube, EFN72687.1.
Phylogenetic trees of IRAKs based on the sequence of kinase domains (A) and death domains (B). The tree was constructed using the neighbor-joining method. Numbers on the lines indicate the percentage bootstrap values for 1000 replicate. Maximum likelihood trees were also constructed with the Posisson model, which were similar to neighbor-joining trees (data not shown). The accession numbers for the used protein sequences are as follows: Homo sapiens IRAK1, NP_001560.2; H. sapiens IRAK2, NP_001561.3; H. sapiens IRAK3, NP_009130.2; H. sapiens IRAK4, NP_001107654.1; H. sapiens TAK1, NP_663304.1; H. sapiens TBK1, NP_037386.1; H. sapiens IKKa, NP_001269.3; H. sapiens IKKb, NP_001547.1; Mus musculus IRAK1, NP_001171444.1; M. musculus IRAK2, NP_751893.3; M. musculus IRAK3, NP_082955.2; M. musculus IRAK4, NP_084202.2; M. musculus TAK1, NP_0033342.1; M. musculus TBK1, NP_062760.3; M. musculus IKKa, NP_031726.2; M. musculus IKKb, O88351.1; DrIRAK1, XP_005166760.1; DrIRAK3, XP_003198307.2; DrIRAK4, NP_956457.1; D. rerio TAK1, NP_001018586.1; D. rerio TBK1, NP_001038213.2; D. rerio IKKa, NP_956611.1; D. rerio IKKb, NP_001116737.1; D. melanogaster Pelle, NP_476971.1; D. melanogaster Tube, NP_001189164.1; D. melanogaster TAK1, NP_524080.1; D. melanogaster IKKb, NP_524751.3; Ciona intestinalis IRAK4, XP_002122012.1; Saccoglossus kowalevskii IRAK4, XP_002739502.1; Caenorhabditis elegans PIK1, NP_001255742.1; Gallus gallus IRAK2, NP_001025776.1; G. gallus IRAK4, NP_001025909.1; Anolis carolinensis IRAK1, XP_008102118.1; A. carolinensis IRAK2, XP_008115205.1; A. carolinensis IRAK4, XP_003221322.2; Xenopus tropicalis IRAK1, NP_001006713.1; X. tropicalis IRAK2, NP_001007501.1; Xenopus tropicalis IRAK4, NP_001116877.1; Eriocheir sinensis Pelle, AKJ26285.1; E. sinensis Tube, AGT21373.1; Scylla paramamosain (SpPelle), AGY49577.1; S. paramamosain (SpTube), AGY49576.1; Camponotus floridanu Pelle, XP_011265220.1; C. floridanu Tube, EFN72687.1.
Both trees suggest that the IRAK family could be separated into three groups: the IRAK4 group, which is conserved between protostomes and deuterostomes; the IRAK1/2/3 group, which is specific to vertebrates; and the Pelle group, which is also conserved between protostomes and invertebrate deuterostomes but lost in vertebrate deuterostomes. The Pelle group tends to represent the intermediate state between IRAK4 and IRAK1/2/3. According to the two phylogenetic trees, amphioxus IRAK4 belongs to the IRAK4 group, whereas amphioxus Pelle is closer to invertebrate Pelle than to vertebrate IRAK1/2/3 (Fig. 2). In line with this, we found that amphioxus Pelle and other invertebrate Pelle have not evolved a distinct C-terminal region, which is presented in vertebrate IRAK1/2/3 (Fig. 1B). It is worth noting that in the kinase-based tree, amphioxus Pelle is clustered with invertebrate Pelle with high confidence, whereas in the death domain based tree, amphioxus Pelle is slightly closer to vertebrate IRAK4s (although not as close as amphioxus IRAK4 to vertebrate IRAK4s) (Fig. 2). This suggests that the death domain of Pelle might evolve slower than its kinase domain. Taken together, the IRAK family of amphioxus and other invertebrates does not have an exact functional counterpart for mammalian IRAK1/2/3.
Taken together, IRAK4, Pelle, and IRAK1/2/3 are three major IRAK lineages in animals. IRAK4 is an ancient lineage conserved between invertebrate and vertebrates. Pelle is also an ancient lineage but lost in vertebrates. IRAK1/2/3 is a vertebrate-specific lineage. Judging from the sequence and phylogenetic positions (Fig. 2), we suspect there might be a chance that IRAK1/2/3 was derived from Pelle and diverged in later evolution.
Expression profile of amphioxus IRAK4 and Pelle
qRT- PCR was used to investigate the expression patterns of BbIRAK4 and BbPelle in different tissues and during immune challenge (Fig. 3). Both BbIRAK4 and BbPelle show relatively high expression in the gill and skin and medium expression in the intestine (Fig. 3A, 3B). These tissues, especially the gill and the skin, are considered the first defense line against microbial pathogens. Indeed, both BbIRAK4 and BbPelle were upregulated in the gill and the intestine after immune challenge with LPS, the common cell wall component of bacteria (Fig. 3C–F). This expression pattern suggests the participation of BbIRAK4 and BbPelle in acute immune reaction. In addition, BbIRAK4 displayed very high expression in the ovary, whereas BbPelle was nearly absent (Fig. 3A, 3B). Because the harvested ovary mainly consisted of unspawned eggs, this observation suggests that BbIRAK4, but not BbPelle, is abundant as maternal transcripts and may have an important role in embryogenesis. In line with this, abundant maternal transcripts of IRAK4 were found in half-smooth tongue sole, zebrafish, and abalone (26, 27, 67). Moreover, the Drosophila Tube, the homolog of IRAK4, has been shown to participate in the formation of the dorsoventral axis during development (66).
Expression profile of BbIRAK4 and BbPelle. (A) Tissue expression of BbIRAK4. (B) Tissue expression of BbPelle. (C and D) Expression profile of BbIRAK4 (C) or BbPelle (D) in amphioxus gill after LPS challenge. (E and F) Expression profile of BbIRAK4 (E) and BbPelle (F) in amphioxus intestine after LPS challenge. Expression level changes were calculated by normalization to the expression of GAPDH. The samples at time 0 h were untreated with either PBS or LPS. Tissues from five fishes were collected for each sample. Two parallel experiments were performed in triplicate, and data are plotted as mean ± SD (n = 3), *p < 0.05, **p < 0.01, ***p < 0.001.
Expression profile of BbIRAK4 and BbPelle. (A) Tissue expression of BbIRAK4. (B) Tissue expression of BbPelle. (C and D) Expression profile of BbIRAK4 (C) or BbPelle (D) in amphioxus gill after LPS challenge. (E and F) Expression profile of BbIRAK4 (E) and BbPelle (F) in amphioxus intestine after LPS challenge. Expression level changes were calculated by normalization to the expression of GAPDH. The samples at time 0 h were untreated with either PBS or LPS. Tissues from five fishes were collected for each sample. Two parallel experiments were performed in triplicate, and data are plotted as mean ± SD (n = 3), *p < 0.05, **p < 0.01, ***p < 0.001.
Amphioxus IRAK4 inhibited MyD88 and TRAF6 by different mechanisms
In agreement with the findings in grouper (35, 68), BbIRAK4 was dispersedly expressed in the cytoplasm (Supplemental Fig. 2). Thus, BbIRAK4 may exert its function in the cytoplasm but not in the nucleus, like mammalian IRAK4. Using luciferase reporter assays, we found that transfecting human HEK293T cell line with high doses (100–300 ng) of BbIRAK4 slightly activated NF-κB signal in an unstable way (Fig. 4A), which we suspected could be due to the minor perturbation conferred by the kinase activity of BbIRAK4. In any case, this warrants further study. In contrast, when B. belcheri MyD88 (BbMyD88) was cotransfected with BbIRAK4 even only 5 ng, the BbMyD88-induced NF-κB activation could be significantly impaired. Morever, a nearly complete blockage could be achieved with a slighted larger dose (50 ng) of BbIRAK4 (Fig. 4B). This observation coincides with the previous reports that bony fish IRAK4 could suppress NF-κB activation in human cells (34, 35). Meanwhile, B. belcheri TRAF6 (BbTRAF6)–mediated NF-κB activation was also suppressed by BbIRAK4, but the effect is limited and not correlated with the dosage of BbIRAK4 (Fig. 4C), hence implying that BbIRAK4 might not act on BbTRAF6 directly. Then we detected the expression levels of TNF-α and IL-8, which are downstream genes of NF-κB after transfection. In line with the reporter assays, BbIRAK4 could intensely downregulate the expression of TNF-α and IL-8 induced by BbMyD88 and also impair their expression induced by BbTRAF6 (Fig. 5A–D).
Regulatory role of BbIRAK4 and BbPelle in NF-κB activity. HEK293T cells were cotransfected with NF-κB transcriptional luciferase reporter plasmid and Renilla luciferase reporter plasmid, together with BbIRAK4/Pelle/MyD88 vectors. (A) BbIRAK4 slightly activated NF-κB signal. (B and C) BbIRAK4 impaired BbMyD88-induced (B) and BbTRAF6-induced (C) NF-κB activation. (D) BbPelle could not activate NF-κB signal. (E and F) BbPelle suppressed BbMyD88-induced (E) and BbTRAF6-induced (F) NF-κB activation. The results are expressed as means ± SD of three samples per treatment, and were confirmed by at least three separate experiments. *p < 0.05, **p < 0.01, ***p < 0.001.
Regulatory role of BbIRAK4 and BbPelle in NF-κB activity. HEK293T cells were cotransfected with NF-κB transcriptional luciferase reporter plasmid and Renilla luciferase reporter plasmid, together with BbIRAK4/Pelle/MyD88 vectors. (A) BbIRAK4 slightly activated NF-κB signal. (B and C) BbIRAK4 impaired BbMyD88-induced (B) and BbTRAF6-induced (C) NF-κB activation. (D) BbPelle could not activate NF-κB signal. (E and F) BbPelle suppressed BbMyD88-induced (E) and BbTRAF6-induced (F) NF-κB activation. The results are expressed as means ± SD of three samples per treatment, and were confirmed by at least three separate experiments. *p < 0.05, **p < 0.01, ***p < 0.001.
Downstream genes were inhibited by BbIRAK4 and BbPelle. HEK293T cells were cotransfected with BbIRAK4/Pelle/MyD88/TRAF6 vectors. The expression level of TNF-α and IL-8 were detected through qRT-PCR assay. (A and B) The expression change of downstream genes induced by BbMyD88, when cotransfected with BbIRAK4 or BbPelle. (C and D) The expression change of downstream genes induced by BbTRAF6, when cotransfected with BbIRAK4 or BbPelle. The results are expressed as means ± SD of three samples per treatment and were confirmed by at least three separate experiments. *p < 0.05, **p < 0.01, ***p < 0.001. DD, death domain.
Downstream genes were inhibited by BbIRAK4 and BbPelle. HEK293T cells were cotransfected with BbIRAK4/Pelle/MyD88/TRAF6 vectors. The expression level of TNF-α and IL-8 were detected through qRT-PCR assay. (A and B) The expression change of downstream genes induced by BbMyD88, when cotransfected with BbIRAK4 or BbPelle. (C and D) The expression change of downstream genes induced by BbTRAF6, when cotransfected with BbIRAK4 or BbPelle. The results are expressed as means ± SD of three samples per treatment and were confirmed by at least three separate experiments. *p < 0.05, **p < 0.01, ***p < 0.001. DD, death domain.
Further coimmunoprecipitation (Co-IP) assays showed that like its mammalian counterparts, BbIRAK4 could use its death domain to directly associate with BbMyD88 (Fig. 6A), suggesting that BbIRAK4 might suppress MyD88 by direct interaction. Immunofluorescence analysis in Hela cells confirmed that BbIRAK4 could colocalize with BbMyD88 in the cytoplasm (Fig. 6E, Supplemental Fig. 3). In contrast, the direct interaction between BbIRAK4 and BbTRAF6 was not observed (Fig. 6C, 6E), suggesting that BbIRAK4 might not affect NF-κB activation through the direct interaction with BbTRAF6.
Functional analysis of BbIRAK4 and BbPelle. For colocalization assay, Hela cells were transfected with BbIRAK4/Pelle/MyD88/TRAF6 vectors. For Co-IP and reporter assays, HEK293T cells were used. (A and B) Co-IP assays showed the interaction between BbMyD88 with BbIRAK4 (A) or BbPelle (B). (C) Co-IP assays showed no interaction between BbTRAF6 with BbIRAK4 or BbPelle. (D) Co-IP assays showed no interaction between BbIRAK4 with BbPelle. (E) Colocalization of BbIRAK4, BbPelle, BbMyD88, and BbTRAF6 in Hela cells. Scale bar, 5 μm. (F) Truncated mutants were constructed. (G and H) The function of BbIRAK4 truncated mutants were analyzed through reporter assays. (I and J) The function of BbPelle truncated mutants were analyzed through reporter assays. Results were confirmed by at least three separate experiments The results of reporter assays are expressed as means ± SD of three samples per treatment. *p < 0.05, **p < 0.01, ***p < 0.001. DD, death domain.
Functional analysis of BbIRAK4 and BbPelle. For colocalization assay, Hela cells were transfected with BbIRAK4/Pelle/MyD88/TRAF6 vectors. For Co-IP and reporter assays, HEK293T cells were used. (A and B) Co-IP assays showed the interaction between BbMyD88 with BbIRAK4 (A) or BbPelle (B). (C) Co-IP assays showed no interaction between BbTRAF6 with BbIRAK4 or BbPelle. (D) Co-IP assays showed no interaction between BbIRAK4 with BbPelle. (E) Colocalization of BbIRAK4, BbPelle, BbMyD88, and BbTRAF6 in Hela cells. Scale bar, 5 μm. (F) Truncated mutants were constructed. (G and H) The function of BbIRAK4 truncated mutants were analyzed through reporter assays. (I and J) The function of BbPelle truncated mutants were analyzed through reporter assays. Results were confirmed by at least three separate experiments The results of reporter assays are expressed as means ± SD of three samples per treatment. *p < 0.05, **p < 0.01, ***p < 0.001. DD, death domain.
Amphioxus Pelle inhibited MyD88 and TRAF6 by different mechanisms
BbPelle was also expressed in the cytoplasm like BbIRAK4 (Supplemental Fig. 2). But, BbPelle could not activate NF-κB even with large dosage (Fig. 4D), which is different to Drosophila Pelle and human IRAK1 (43, 69). BbPelle also impaired the NF-κB activation when cotransfected with BbMyD88, although this required a larger dosage (Fig. 4E). BbPelle also negatively regulated BbTRAF6-mediated NF-κB activation, but different from BbIRAK4, the suppressive effect of BbPelle on BbTRAF6 was correlated with the transfection dosage (Fig. 4F). Further qRT-PCR assays showed that BbPelle impaired the expression of TNF-α induced by BbMyD88 or BbTRAF6 and the expression of IL-8 induced by BbTRAF6 (Fig. 5A–D).
Co-IP assays showed the direct binding between BbPelle and BbMyD88 via the death domains (Fig. 6B) but failed to show the stable interaction between BbPelle and BbTRAF6 (Fig. 6C). The latter observation might be due to the fact that BbPelle lacks the C-terminal region and TRAF6-binding motifs, which are required for the interaction between vertebrate IRAK1/2 and TRAF6 (Fig. 1B). BbPelle also colocalized with BbMyD88 in the cytoplasm of Hela cells (Fig. 6E, Supplemental Fig. 3), in agreement with the Co-IP assay. Taken together, we suggested that BbPelle acts on BbMyD88 to repress its activity. We also suggested that BbPelle might act on both upstream and downstream of BbTRAF6, because BbPelle negatively regulates BbTRAF6 in a dosage-dependent manner without direct contact. As reference, there was a negative feedback regulation role of Drosophila Pelle in the Toll pathway (70). Drosophila Pelle possibly direct phosphorylation of Tube, leading to the disaggregation of membrane-associated Tube clusters (70).
No interaction detected between amphioxus IRAK4 and Pelle
In both Drosophila and vertebrates, Tube (IRAK4 in vertebrates) can interact with both Pelle (IRAK1/2 in vertebrates) and MyD88 to form the so-called Myddosome complex during signal transduction (Fig. 1A). Here, we performed Co-IP assays to test if BbPelle and BbIRAK4 could interact with each other for the regulation of the MyD88-TRAF6 pathway. Surprisingly, although both BbPelle and BbIRAK4 could interact with BbMyD88 (Fig. 6A, 6B), we did not observe the direct interaction between BbPelle and BbIRAK4 (Fig. 6D). However, part of BbPelle and BbIRAK4 could colocalize with each other (Fig. 6E, Supplemental Fig. 3), indicating that they may form complexes through the interaction with MyD88. This suggested that a new functional model for the IRAK proteins in the basal chordate.
Functions of different domains in BbIRAK4 and BbPelle
To investigate which domain was essential for their suppressive activity, we constructed several truncated mutants of BbIRAK4 and BbPelle (Fig. 6F). Two truncated mutants of BbIRAK4 with death domain still strongly suppressed BbMyD88-induced NF-κB activation, whereas the suppressive activity was very weak when cotransfected BbMyD88 with BbIRAK4-△DD (Fig. 6G). These results indicated that the death domain of BbIRAK4 was mainly responsible for its suppressive function on BbMyD88. However, BbIRAK4-△KD (or other truncated mutants) had no (or weak) influence on BbTRAF6-mediated NF-κB activation when cotransfected with BbTRAF6, suggesting that the kinase domain of BbIRAK4 might be responsible for this function (Fig. 6H).
As for BbPelle, both the death domain and the ProST region contributed to the inhibitory effect on BbMyD88, whereas the kinase domain was not required (Fig. 6I). Interestingly, BbPelle-△ProST sharply enhanced both the BbMyD88-induced and the BbTRAF6-induced NF-κB activation (Fig. 6I, 6J), suggesting that the ProST domain was important for its suppressive function. Based on these results, we suggested that BbPelle could bind to MyD88 and then use the ProST region to inhibit the downstream NF-κB activation.
Ex vivo functional analysis of zebrafish IRAKs
We also analyzed the function of zebrafish Danio rerio IRAK (DrIRAK)1/3/4 in HEK293T cells. DrIRAK1 could significantly activate NF-κB signal, whereas DrIRAK3 and DrIRAK4 could not (Fig. 7A). Both DrIRAK3 and DrIRAK4 impaired D. rerio MyD88 (DrMyD88)–induced NF-κB activation with dosage dependent manners, which coincided with the reports of other teleost IRAK4 (Fig. 7B, 7C) (34, 35). When cotransfected with DrIRAK1, DrIRAK3 prominently downregulated DrIRKA1-mediated NF-κB activation (Fig. 7D), whereas DrIRAK4 had no effect on DrIRKA1-mediated NF-κB activation (Fig. 7E). As mammalian IRAK3 prevented IRAK1/2 from the dissociation from MyD88 and negatively regulated the pathway (23, 24), we presumed that DrIRAK3 may perform a similar function as mammalian IRAK3. However, the role of DrIRAK4 was more close to amphioxus IRAK4. We also tried to use ZFL cells, derived from zebrafish liver, to conduct the luciferase reporter assays, but the low transfection efficiency made the result not sufficiently reliable and repeatable (data not shown).
Regulatory role of DrIRAK1/3/4 in NF-κB activity. HEK293T cells were cotransfected with NF-κB transcriptional luciferase reporter plasmid and Rellina luciferase reporter plasmid, together with DrIRAK1/IRAK3/IRAK4/MyD88 vectors. (A) DrIRAK1 could significantly active NF-κB signal. (B) DrIRAK3 impaired DrMyD88-induced NF-κB activation. (C) DrIRAK4 impaired DrMyD88-induced NF-κB activation. (D) DrIRAK3 impaired DrIRAK1-induced NF-κB activation. (E) DrIRAK4 had no influence on DrIRAK1-induced NF-κB activation. Results were confirmed by at least three separate experiments. The results are expressed as means ± SD of three samples per treatment. **p < 0.01, ***p < 0.001.
Regulatory role of DrIRAK1/3/4 in NF-κB activity. HEK293T cells were cotransfected with NF-κB transcriptional luciferase reporter plasmid and Rellina luciferase reporter plasmid, together with DrIRAK1/IRAK3/IRAK4/MyD88 vectors. (A) DrIRAK1 could significantly active NF-κB signal. (B) DrIRAK3 impaired DrMyD88-induced NF-κB activation. (C) DrIRAK4 impaired DrMyD88-induced NF-κB activation. (D) DrIRAK3 impaired DrIRAK1-induced NF-κB activation. (E) DrIRAK4 had no influence on DrIRAK1-induced NF-κB activation. Results were confirmed by at least three separate experiments. The results are expressed as means ± SD of three samples per treatment. **p < 0.01, ***p < 0.001.
Discussion
The IRAK family plays a crucial role in the TLR/IL-1R–MyD88–TRAF6–NF-κB signaling pathway, but their function appears not consistent between Drosophila and vertebrates, and even between bony fishes and mammals (25–30). As amphioxus occupied a unique evolutionary position between vertebrates and invertebrates, the study of amphioxus IRAKs may provide insight into the functional evolution of IRAK family (71). In this study, we showed that amphioxus has only two IRAK genes; one is the ortholog to IRAK4, whereas the other one is the ortholog to invertebrate Pelle, rather than to vertebrate IRAK1/2/3. Phylogenetic trees show that the kinase domain of amphioxus Pelle might evolve faster than its death domain. As we known, the death domain is served as a bridge to link MyD88 and IRAKs for forming Myddosome complex, whereas the kinase domain of IRAKs phosphorylates downstream proteins or itself during signaling transduction (6, 7). Thus, these findings suggested that although the IRAKs have a conserved position in TLR pathway, the specific functions significantly vary between major animal taxa.
We cloned these two IRAK genes from the amphioxus B. belcheri. Both BbIRAK4 and BbPelle may be related to amphioxus immune reaction according to their high expression in gill and skin. The gill of amphioxus belongs to digestive system and is the first defense line of amphioxus immunity, and the skin of amphioxus is also a potential immune organ of amphioxus (71). We observed abundant maternal transcripts of BbIRAK4 in eggs, which has also been reported in half-smooth tongue sole, zebrafish, and abalone (26, 27, 67), suggesting a potentially important role in embryogenesis. For example, it is reported that Drosophila Tube, the homolog of IRAK4, participates in the formation of the dorsoventral axis during development (66). Although BbPelle transcripts were nearly absent in eggs, they indicated an obviously functional difference of amphioxus IRAK4 and Pelle.
BbIRAK4 slightly induced activation of NF-κB, consistent with previous findings of bony fish IRAK4 (34, 35). Amphioxus or bony fish IRAK4 may be not fully compatible with the mammalian system, because zebrafish IRAK4 could significantly stimulate NF-κB activation in ZFL cells (26). BbIRAK4 showed a negative role in the MyD88–TRAF6–NF-κB pathway, in agreement with the previous studies of bony fish IRAK4 (34, 35). Both BbIRAK4 and BbPelle could interact with BbMyD88 through their death domain, suggesting a possibility to form a complex similar to mammal Myddosome complex. Further truncated mutants analysis confirmed that the death domain of BbIRAK4 and BbPelle is essential for suppressing BbMyD88-induced NF-κB activation. However, we could not detect the direct interaction between BbIRAK4 and BbPelle. Thus, the Myddosome-like complex in amphioxus might be IRAK4–MyD88–Pelle but not Pelle–IRAK4–MyD88. This may explain the functional difference of amphioxus IRAKs with mammals. We also observed a suppressive role of BbPelle in MyD88–TRAF6–NF-κB pathway. Also, Drosophila Pelle was suggested to directly phosphorylate Tube, causing the dissociation of Pelle–Tube–MyD88 complex, and downregulate the pathway (70). However, we could not detect the interaction between BbIRAK4/BbPelle and BbTRAF6, which may be due to the lacking of the C-terminal domain and the TRAF6 binding sites. Interestingly, when eliminated the ProST domain of BbPelle, BbPelle-△ProST enhanced the NF-κB activation induced by BbMyD88 or BbTRAF6. The ProST domain is rich in serines, prolines, and threonines, and mammalian IRAK1 was reported to undergo hyperphosphorylation in ProST domain, which was important for its signal transduction (8, 72). Thus, we suspect that hyperphosphorylation in BbPelle ProST domain was essential for its suppressive function. However, BbPelle without kinase domain still downregulated BbMyD88 and BbTRAF6-mediated NF-κB activation, indicating that the ProST domain of BbPelle may be phosphorylated by some other kinases but not autophosphorylation. In HEK293T cells, phosphorylated BbPelle may induce the degradation of BbTRAF6 or directly inhibit some endogenous downstream signal molecules, such as NFKB1 and NFKB2 which also contained death domain, and performed as inhibition function when cotransfected with BbTRAF6. And these require further investigation. Considering that the kinase domain of BbIRAK4 is essential for inhibiting BbTRAF6 induced NF-κB activation, one possible explanation in amphioxus is that BbPelle maybe phosphorylated by BbIRAK4 when they interacted with BbMyD88 at the same time.
We also analyzed the function of zebrafish IRAK1/3/4, indicating that the function of zebrafish IRAK1 and IRAK3 was similar to their mammalian counterparts, whereas zebrafish IRAK4 was more similar to amphioxus IRAK4. Thus, there was no functional counterpart of amphioxus Pelle in vertebrate, and the functional difference of amphioxus and mammalian IRAK4 was evolved after the divergence of bony fishes and amphibian.
In summary, our findings show that the composition and sequence of amphioxus IRAKs are more similar to invertebrates than vertebrates. However, molecular functions of amphioxus IRAKs are quite different from Drosophila and mammals. Together with previous reports, we concluded that although IRAKs are conserved members of the TLR/IL-1R pathways, their functions (unlike those of MyD88 and TRAF6) are highly variable during evolution and clearly specialized in different major animal taxa. Considering their kinase activity, we further suggest that the functional variability of IRAKs might confer important plasticity to the intracellular signal transduction of the TLR/IL-1R pathways, which in return could facilitate the species to evolve and to adapt during the completion against the pathogens.
Footnotes
This work was supported by the National Key R&D Program of China (2018YFD0900503), National Natural Science Foundation Projects 31800729, 31872595, and 1722052, the Marine S&T Fund of Shandong Province (2018SDKJ0302-2), projects from Guangdong (201804010434, 201804020039, 2018A030310157, and 2017B030314021), and by the National Supercomputer Center in Guangzhou and the Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund.
The online version of this article contains supplemental material.
Abbreviations used in this article:
References
Disclosures
The authors have no financial conflicts of interest.