Genome-wide gene expression analysis of amphioxus (Branchiostoma belcheri) following lipopolysaccharide challenge using strand-specific RNA-seq.

Amphioxus is the closest living proxy for exploring the evolutionary origin of the immune system in vertebrates. To understand the immune responses of amphioxus to lipopolysaccharide (LPS), 5 ribosomal RNA (rRNA)-depleted libraries of amphioxus were constructed, including one control (0 h) library and 4 treatment libraries at 6, 12, 24, and 48 h post-injection (hpi) with LPS. The transcriptome of Branchiostoma belcheri was analyzed using strand-specific RNA sequencing technology (RNA-seq). A total of 6161, 6665, 7969, and 6447 differentially expressed genes (DEGs) were detected at 6, 12, 24, and 48 hpi, respectively, compared with expression levels at 0 h. We identified amphioxus genes active during the acute-phase response to LPS at different time points after stimulation. Moreover, to better visualize the resolution phase of the immune process during immune response, we identified 6057 and 5235 DEGs at 48 hpi by comparing with 6 and 24 hpi, respectively. Through real-time quantitative PCR (qRT-PCR) analysis of 12 selected DEGs, we demonstrated the accuracy of the RNA-seq data in this study. Functional enrichment analysis of DEGs demonstrated that most terms were related to defense and immune responses, disease and infection, cell apoptosis, and metabolism and catalysis. Subsequently, we identified 1330, 485, 670, 911, and 1624 time-specific genes (TSGs) at 0, 6, 12, 24, and 48 hpi. Time-specific terms at each of 5 time points were primarily involved in development, immune signaling, signal transduction, DNA repair and stability, and metabolism and catalysis, respectively. As this is the first study to report the transcriptome of an organism with primitive immunity following LPS challenge at multiple time points, it provides gene expression information for further research into the evolution of immunity in vertebrates.

Introduction

Amphioxus, also known as the lancelet, is a living basal chordate (cephalochordates), and it is widely used as a proxy for the common ancestor of cephalochordates and vertebrates., Therefore, amphioxus is a model experimental animal for studying the developmental homology, genome evolution, and comparative immunology of vertebrates.,, Extant amphioxus clades consist of three genera (Branchiostoma, Epigonichthys, and Asymmetron) and a total of ∼31 species, with Branchiostoma encompassing the largest number species (over 28 species). This genus has a wide distribution across the mid-low latitudes of the Atlantic and Pacific oceans., To date, virtually all amphioxus-related studies have been conducted in this genus. For example, there are two species of Branchiostoma with sequenced whole genomes: Branchiostoma floridae and Branchiostoma belcheri, the latter of which has been used in the majority of experiments. These observations not only indicate that Branchiostoma species, especially B. belcheri, have been used as a representative of amphioxus but also that they are useful for investigating scientific questions associated with the evolution of chordates.

As with vertebrates, amphioxus exhibit sensitive immune responses to some microorganisms and compound extracts, such as bacteria (including Vibrio anguillarum, Staphylococcus aureus, and Escherichia coli), lipopolysaccharide (LPS), lipoteichoic acid (LTA), polyhydroxyalkanoates (PHA), and poly (I:C)., Thus, they have been used to explore the evolutionary mechanisms of innate immunity and the origin of adaptive immunity through comparative and functional genomics, transcriptomics, and proteomics.,, Research into the innate immune system has shown that amphioxus shares much of the innate immunity gene repertoire with vertebrates, and systematic genomic analysis has revealed the complexity and diversity of these components., Previous research found that exposure to LPS or microorganisms induced the upregulation of both intercellular and intracellular innate immunity genes in B. belcheri, such as the classical component system protein C1q, MyD88, Toll-like receptors (TLRs), and TNF-associated factor (TRAF). To systematically study the amphioxus expression response under microbial invasion, Huang et al. detected and analyzed the expression profiles of B. belcheri infected by bacteria (S. aureus and Vibrio parahaemolyticus) and control groups using suppression subtractive hybridization (SSH), identifying a set of genes involved in immunity, stress response, digestion, apoptosis, metabolism, and development. Furthermore, the expressed sequence tags (ESTs) approach has also been used to identify immune-related and hepatic cecum-specific genes of Branchiostoma japonicum challenged by LPS, indicating that the hepatic cecum is a key immunological organ., Recently, Liao et al. performed a comparative analysis between B. belcheri at 0 h (control) and 12 h post-injection (hpi) using microarray technology and identified amphioxus miRNAs involved in the acute-phase response to LPS.

Compared with tag- and hybridization-based methods, such as SSH, ESTs, or microarrays, next generation cDNA sequencing technology (RNA-seq) provides the advantages of more complete and updated sequence information, less background noise, and greater dynamic range for gene identification and annotation, discovery of genomic information, and quantification of gene expression levels., In particular, strand-specific RNA-seq can substantially improve the value and quality of data due to direct acquisition of information on the originating strand, allowing for exact identification of antisense transcripts, accurate demarcation of boundaries between adjacent genes, and correction for expression levels of overlapping transcripts. Therefore, strand-specific, massively parallel RNA-seq is a powerful tool for analysis of expression profiling., Additionally, library construction involving the removal of ribosomal RNA (rRNA) avoid shaving to use magnetic beads with polyT for mRNA enrichment, thereby fully retaining mRNAs that may otherwise be eluted. In addition, more genes with low levels of expression have been observed with rRNA-depletion methods.,

Until now, no studies have used RNA-seq to analyze dynamic variations in the transcriptome of amphioxus over multiple time points post-injection with an antigen. Detailing of the expression profiles of genes involved in the immune response in amphioxus will be valuable for developing a deeper understanding of the evolution of the vertebrate immune system. Based on this, researchers will be able to confirm whether candidate genes possess immune responses, providing preliminarily reference for further experiments. Moreover, by comparing homologous gene expression between amphioxus and vertebrates, researchers will be able to better trace the evolutionary origin of key immune genes, pathways, and even tissues in vertebrates, an approach that has been used in previous studies. The expression profile of amphioxus against an immune stimulator remains largely unclear, causing difficulty in selecting effective immune-related genes and conducting comparative expression analysis between amphioxus and vertebrates. Here, we constructed five transcriptomic libraries of B. belcheri, one at 0 h (control) and four at various time points following LPS injection using rRNA depletion methods. Subsequently, gene expression profiles were analyzed using strand-specific RNA-seq. The aims of this study were to provide a basic view of the ancient immune response in amphioxus and a valuable genomic resource for further studying the origin and evolution of immune-related genes in vertebrates.

Materials and methods

Sample collection

Healthy B. belcheri were obtained from the Beihai Marine Station of Nanjing University (Beihai, China), and individuals of nearly unified sizes were selected. Amphioxus were maintained in acrylic tanks with filtered seawater for a week to empty the contents of their digestive systems and were acclimatized at 24–28°C. Subsequently, we injected LPS in PBS suspension (1mg/ml) into the body coelom of each amphioxus, and the same dosage of PBS was injected into controls according to our previous studies., The injected amphioxus were reared at 26°C, and six individuals were collected at each of 5 time points post-injection: 0 (control), 6, 12, 24, and 48 hpi. All individuals were preserved in TRIzol reagent (Invitrogen, USA) at −80°C for storage before RNA extraction.

RNA extraction and Illumina sequencing

Total RNA was extracted from each individual following previous methods, and RNase-free DNase Set (Qiagen, Germany) was used to remove residual genomic DNA. RNA concentration and quality were detected using a BioPhotometer Plus spectrophotometer (Eppendorf, Germany). Only RNA with an OD260/D280 absorbance between 1.8 and 2.0 and with structural integrity verified by gel electrophoresis was used for further experiments. Furthermore, we used the 2100 Bioanalyzer (Agilent Technologies, USA) to verify the quantity and integrity of the RNA. Total RNA of six individuals collected at each time point was pooled equally. Library construction and sequencing were performed at the Beijing Genomics institute (BGI-Shenzhen, China) under standard procedures. Briefly, we removed rRNA from the total RNA. The remaining RNAs were fragmented into short fragments (around 200–500nt), and then first-strand cDNA was synthesized from fragments by random hexamer primers with dUTP substituted for dTTP during the synthesis of the second strand. Ends of short fragments were processed by adding a single adenine after purification and connected with adapters. Finally, the second strand was degraded using UNG (uracil-N-glycosylase). The paired-end cDNA library was sequenced on the HiSeq 2000 platform (Illumina, USA) that generated paired-end reads of 90 bp.

Sequence filtering, mapping, assembly, and annotation

The raw data were cleaned by filtering out adaptor sequences, low-quality reads with unknown bases (N bases) >10%, and sequences whose proportion of low-quality bases (Q value ≤ 10) was more than 50% using SOAPnuke software (https://github.com/BGI-flexlab/SOAPnuke). Clean paired-end reads were mapped to the B. belcheri genome (v18h27.r3) using TopHat, allowing no more than five mismatches per read in the alignment. The expression level of each gene was determined by Cufflinks assembler based on reads per kilobase per million mappable read (FPKM) values. Using transcripts from the reference genome, we annotated gene functions using the NCBI non-redundant protein sequence (Nr), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases.

Identification of differentially expressed and time-specific genes

We retained expressed genes (FPKM > 0.1) for further analysis., The fold-change of gene expression was calculated using the formula log2(FPKM-control/FPKM-treatment). The differentially expressed genes (DEGs) and their significance levels (P-values) between library pairs (comparison groups, CGs; CG1, control versus 6 hpi; CG2, control vs. 12 hpi; CG3, control vs. 24 hpi; CG4, control vs. 48 hpi; CG5, 6 hpi vs. 48 hpi; CG6, 24 hpi vs. 48 hpi) were detected by the DEGseq package based upon FPKM values. To avoid false-positive results, we performed a multiple testing correction on significance levels by controlling the false discovery rate (FDR) according to the Benjamini and Hochberg (BH) method using the p.adjust module of the R package. Finally, if the ratio of gene expression was ≥ 2 (absolute value of log2ratio ≥ 1) and the FDR was < 0.001, DEGs were considered reliable targets for further analyses. We used the index τ to estimate the time specificity of gene expression according to previous methods using a custom Perl script., Genes with τ > 0.85 were considered time-specific genes (TSGs)., Subsequently, Jensen-Shannon (JS) scores were used to evaluate time points of specific expression for each TSG based on previous descriptions of the JS algorithm.,

Enrichment of functions and pathways

We conducted GO enrichment analyses based on the gene annotation results and performed Fisher's exact tests using the Blast2GO pipeline. In addition, KEGG pathway enrichment analysis was performed using KOBAS 2.0 software. The significance levels of GO terms and KEGG pathways were corrected by FDR, and only terms and pathways with an FDR < 0.05 were considered as significantly enriched targets.

Quantitative real-time PCR

We conducted quantitative real-time PCR (qRT-PCR) analysis of 12 randomly selected DEGs (significantly different expression in all 4 CGs) using the same RNA samples that were used in sequencing. First-strand cDNA was synthesized using the RevertAid™ First Strand cDNA Synthesis Kit (Thermo Scientific, USA) following the product manual. Beacon Designer 7 was used to design gene-specific primers for qRT-PCR (Supplementary Table S1). The cDNA products (100 ng/μl) were amplified by SYBR® Premix Ex Taq II Kit (Takara, Japan) on the ABI PRISM 7300 Real-Time PCR System (Applied Biosystems, USA). We conducted each qRT-PCR in triplicate, and then relative gene expression levels were calculated using the 2−△△CT method and EF-1A as a reference gene. Consistency of fold-changes between RNA-seq and qRT-PCR were estimated by calculating Pearson's correlation coefficient (r) and its significance (P-value) using IBM SPSS statistics 22.

Results

Genome-guided mapping

Strand-specific RNA-seq was performed for B. belcheri 0 (control), 6, 12, 24, and 48 hpi with LPS. We obtained an equal number of reads (∼95 million clean reads) for each library. Using the B. belcheri genomic sequence (http://genome.bucm.edu.cn/lancelet/), 77–80% of the clean reads in each of the five libraries were mapped to the whole B. belcheri genome, and the 42–44% of the clean reads were mapped to the 30,392 protein-coding genes. We detected 22,518 (74.1%), 21,171 (70.4%), 21,770 (71.6%), 22,027 (72.5%), and 21,501 (70.7%) expressed genes in the 0, 6, 12, 24, and 48 hpi libraries, respectively, for downstream analysis.

Detection of DEGs

A total of 6161, 6665, 7969, and 6447 DEGs were identified in CG1, CG2, CG3, and CG4, respectively (Supplementary Table S2). Of these genes, 2,058 were upregulated and 4,103 were downregulated at 6 hpi compared with control expression levels (CG1); 3,507 and 3,158 were up- and downregulated, respectively, in CG2; 4,879 and 3,090 were up-and downregulated, respectively, in CG3; and, 3,358 and 3,089 were up- and downregulated, respectively, in CG4 (Fig. 1A). As the greatest number of DEGs was observed in CG3, the top 10 up- and downregulated genes from this CG are presented in Table 1. The top 20 up- and downregulated genes for all CGs can be seen in Supplementary Table S3. The most upregulated gene was tumor necrosis factor receptor superfamily member 16 (NGFR), which was observed in 3 CGs (CG2, 3, and 4). Some genes encoding well-known proteins involved in immune response were upregulated only in a single CG, such as complement C1s subcomponent (C1S), coagulation factors, immunity-related GTPase family M protein (IRGM), E3 ubiquitin-protein ligase, neurogenic locus Notch protein (Notch), mucin, and interferon-induced protein. Additionally, we frequently detected genes encoding serine/threonine-protein kinase (Nek7), mitofusin-2 (MFN2), complement component C6 (C6), E3 ubiquitin-protein ligase, SCO-spondin (SSPO), and adhesion G protein-coupled receptors (GPCRs) among the top 20 most downregulated genes in multiple CGs. Some genes were annotated as encoding hypothetical proteins and were thereby assigned unidentified or unknown functions. Furthermore, 1,654 DEGs were shared among CG1–4, among which 502 genes were consistently upregulated compared with controlsand 970 were downregulated (Fig. 1B); interestingly, the remaining 182 DEGs (gene set 1, GS1) exhibited divergent expression changes once at least among CGs, indicating that GS1 exhibited the most acute response to LPS immunity. In addition, to better visualize the resolution phase of the immune response, we also conducted comparison of gene expression levels between time points post-injection with LPS (during immune response, CG5 and CG6) (Supplementary Table S2), 4,207 and 1,850 DEGs (a total of 6057) were up- and downregulated at 48 hpi compared with 6 hpi expression levels (CG5), respectively; 1,926 and 3,309 (a total of 5235) were up-and downregulated, respectively, in CG6. In the top 20 up- and downregulated genes (Supplementary Table S3), 3 genes showed upregulated expression in both CG5 and CG6, including transmembrane emp24 domain-containing protein 9 (TMED9), translation initiation factor IF-2 (INFB), and mucin-4 (MUC4); downregulated expression of 3 genes was shared by these 2 CGs, including death domain-containing protein (CRADD), ribonuclease inhibitor (Rnh1), and coagulation factor XIII B chain (F13B). We also found that some typical immune-related genes were CG-specific between CG5 and CG6, such as caspase-6 (CASP6), MFN2, Complement component C1q receptor (CD93), NF-kappa-B inhibitor β (NFKBIB), and so on.

Figure 1.

A. Summary of the number of up- and downregulated DEGs. B. Venn diagram of DEGs among four comparison groups (CGs). The number of common DEGs is indicated in black, the number of commonly up-and downregulated DEGs is indicated in orange and purple, respectively. C. Venn diagram of GO terms enriched among DEGs in the 4 CGs. D. Venn diagram of KEGG pathways enriched among DEGs in the 4 CGs. CG1, control vs. 6 hours post infection (hpi); CG2, control vs. 12 hpi; CG3, control vs. 24 hpi; CG4, control vs. 48 hpi.

Table 1.

Top 10 most upregulated and downregulated genes between 0 and 24 hours post-LPS injection.

GeneID	log2Ratio(0h/24h)	Annotation	Gene name	FDR
The top 10 most upregulated genes
112430F	-18.22032425	Tumor necrosis factor receptor superfamily member 16	Ngfr	2.76E-11
026750F	-16.0790908	Fibronectin (Fragments)	Fn1	0.000192
248780F	-15.69322215	Trypsin delta/gamma	DeltaTry	0
203350F	-15.58255283	Histone H2A.1	H2A.1	0
002980R	-15.30114853	Transmembrane matrix receptor MUP-4	Mup-4	0
256850R	-14.8359887	Histone H3-like centromeric protein CNA1	CNA1	0
003260R	-14.4435152	Zona pellucida sperm-binding protein 1	Zp1	0
088970R	-14.22843715	CDKN2AIP N-terminal-like protein	Cdkn2aipnl	7.69E-07
018230F	-14.1855803	Pancreatic secretory granule membrane major glycoprotein GP2	Gp2	0
248630F	-13.92488095	Histone H2B.6	Th123	0
The top 10 most downregulated genes
102380F	9.337660986	Transmembrane emp24 domain-containing protein 9	Tmed9	2.05E-114
100790R	9.521525178	SCO-spondin	Sspo	4.02E-64
230910R	9.597579418	Complement component C6	C6	8.14E-67
098080R	9.668176687	Hypothetical protein	/	2.02E-69
034440R	9.860613824	E3 ubiquitin-protein ligase TRIM32	Trim32	5.24E-77
139330R	9.895920693	Dromaiocalcin-1	Dcal1	0
112040R	9.897505202	3-α-hydroxysteroid sulfotransferase	Std1	0
293820R	10.05203989	Kinesin light chain 3	Klc3	2.45E-85
245710F	10.43252017	Hypothetical protein	/	1.68E-104
133300R	10.46494165	Pro-Pol polyprotein	Pol	0

GO and KEGG enrichment analysis of DEGs

GO enrichment analysis was performed for the DEGs of each CG (Supplementary Table S4). A total of 167, 252, 501, and 150 terms were significantly (FDR < 0.05) overrepresented in CG1–4, respectively. Enriched GO terms included 67 (48.20%) biologic process terms, 22 (15.83%) cellular component terms, and 50 (35.97%) molecular functions terms in CG1; 155 (61.51%) biologic process terms, 34 (13.49%) cellular component terms, and 63 (25.00%) molecular functions terms in CG2; 334 (66.67%) biologic process terms, 61 (12.18%) cellular component terms, and 106 (21.16%) molecular function terms for CG3; and 86 (57.33%) biologic process terms, 17 (11.33%) cellular component terms, and 47 (31.33%) molecular functions terms in CG4. As shown in the results (Supplementary Table S4), 76 terms shared by CG1–4 were primarily related to molecule binding, enzyme and protein activity, metabolic process, cell death, response to stimulus, and cell region. Among the remaining 537 GO terms were some evidently immune-related categories that were shared among some CGs or showed specific distribution in a single CG (Fig. 1C), such as defense response (GO: 0006952), inflammatory response (GO: 0006954), response to wounding (GO: 0009611), cell surface receptor signaling pathway (GO: 0007166), immune system development (GO: 0002520), humoral immune response (GO: 0006959), response to bacterium (GO:0009617), response to molecules of bacterial origin (GO: 0002237), immune response (GO:0006955), and T cell differentiation in the thymus (GO:0033077). In addition, we detected GO categories involved in apoptotic processes, such as cell death (GO: 0008219), death (GO: 0016265), apoptotic process (GO: 0006915), regulation of cell death (GO: 0010941), and programmed cell death (GO: 0012501). For GS1, 182 DEGs were significantly enriched in only 2 cellular component terms: extracellular region part (GO: 0044421, FDR = 0.0086) and extracellular region (GO: 0005576, FDR = 0.0112). Additionally, a total of 153, 105 terms were significantly enriched for DEGs in CG5 and CG6, respectively. Obtained GO terms by enrichment analysis included 86 (56.21%) biologic process terms, 25 (16.34%) cellular component terms, and 42 (27.45%) molecular functions terms in CG5; 56 (53.33%) biologic process terms, 16 (15.24%) cellular component terms, and 33 (31.43%) molecular functions terms in CG6 (Supplementary Table S4). Concretely, 20 biologic process terms shared by CG5–6 were primarily associated with response to stimulus, defense response, metabolic process, cellular process, and response to stress. Among the remaining biologic process terms, we detected some evidently apoptosis-related terms that exhibited CG-specific distribution in CG5 or CG6, such as cell death (GO: 0008219), death (GO: 0016265), apoptotic process (GO: 0006915), programmed cell death (GO: 0012501).

We conducted KEGG pathway enrichment analysis of DEGs to detect significantly altered pathways in CGs. In total, 37, 15, 11, and 15 pathways were significantly enriched in CG1–4, respectively (Supplementary Table S5). These pathways were primarily involved in innate immunity, infection, cancer and disease, metabolism, and cell apoptosis, such as the NF-κB signaling pathway (ko04064), cytokine-cytokine receptor interaction (ko04060), RIG-I-like receptor signaling pathway (ko04622), TGF-β signaling pathway (ko04350), cell cycle (ko04110), endocytosis (ko04144),antigen processing and presentation (ko04612), Salmonella infection (ko05132), pathways in cancer (ko05200), metabolic pathways (ko01100), and apoptosis (ko04210). Moreover, one pathway, the NOD-like receptor signaling pathway (ko04621), was shared by all these 4 CGs. Conversely, the number of CG-specific pathways were 23, 6, 9, and 2 for CG1–4, respectively (Fig. 1D). In CG1, specific pathways were primarily related to infection and invasion of microbes, extracellular innate immunity process, and cancer and immunity pathways; more detailed information is listed in Table 2. CG2-specific pathways were involved in intracellular immune signaling, function of immune cell, and cancer. CG3-specific pathways included terms associated with nucleic acid heredity and repair, as well as metabolism and cancer. Two pathways related to fat metabolism were specific to CG4 (Table 2). For two CGs during immune response, in total, 20 and 17 pathways were obtained in CG5–6, respectively (Supplementary Table S5). Seven pathways shared by CG5 and CG6 were primarily related to innate immunity, and cell apoptosis, including the complement and coagulation cascades (ko04610), NOD-like receptor signaling pathway (ko04621), NF-kappa B signaling pathway (ko04064), apoptosis (ko04210), cytokine-cytokine receptor interaction (ko04060), TNF signaling pathway (ko04668), notch signaling pathway (ko04330). The number of CG-specific pathways were 13 and 10 for CG5–6, respectively. However, unlike the pathways shared by these two CGs, these specific terms were primarily involved in cancer and disease, metabolism, cell cycle.

Table 2.

The shared and specific pathways enriched by DEGs among all 4 comparison groups (CG1–4).

Pathway	Pathway ID	FDR
The shared pathways across all CGs
NOD-like receptor signaling pathway	ko04621	2.04E-04
The CG1-specific pathways
Focal adhesion	ko04510	3.60E-07
Dilated cardiomyopathy	ko05414	1.75E-06
ECM-receptor interaction	ko04512	1.18E-05
Arrhythmogenic right ventricular cardiomyopathy (ARVC)	ko05412	9.58E-05
Autoimmune thyroid disease	ko05320	4.46E-04
Prion diseases	ko05020	7.59E-04
RIG-I-like receptor signaling pathway	ko04622	1.05E-03
Haematopoietic cell lineage	ko04640	2.18E-03
Salmonella infection	ko05132	4.01E-03
p53 signaling pathway	ko04115	4.01E-03
Pancreatic secretion	ko04972	5.58E-03
Tyrosine metabolism	ko00350	6.54E-03
Melanoma	ko05218	1.01E-02
Cardiac muscle contraction	ko04260	1.97E-02
Regulation of actin cytoskeleton	ko04810	1.97E-02
Drug metabolism – cytochrome P450	ko00982	2.63E-02
Bacterial invasion of epithelial cells	ko05100	2.71E-02
Phototransduction – fly	ko04745	2.93E-02
Bladder cancer	ko05219	3.25E-02
Pathogenic Escherichia coli infection	ko05130	3.51E-02
Toxoplasmosis	ko05145	3.86E-02
Complement and coagulation cascades	ko04610	4.07E-02
Carbohydrate digestion and absorption	ko04973	4.07E-02
The CG2-specific pathways
Apoptosis	ko04210	1.61E-03
Natural killer cell mediated cytotoxicity	ko04650	1.01E-02
Endocytosis	ko04144	1.17E-02
Antigen processing and presentation	ko04612	1.47E-02
Thyroid cancer	ko05216	3.42E-02
Colorectal cancer	ko05210	4.12E-02
The CG3-specific pathways
DNA replication	ko03030	1.25E-08
Cell cycle	ko04110	2.16E-08
Nucleotide excision repair	ko03420	1.30E-02
Metabolic pathways	ko01100	1.30E-02
Progesterone-mediated oocyte maturation	ko04914	2.98E-02
Glycosylphosphatidylinositol(GPI)-anchor biosynthesis	ko00563	2.98E-02
Proteasome	ko03050	2.98E-02
Mismatch repair	ko03430	3.64E-02
Non-small cell lung cancer	ko05223	4.18E-02
The CG4-specific pathways
Fat digestion and absorption	ko04975	2.86E-02
Biosynthesis of unsaturated fatty acids	ko01040	3.15E-02

Identification of TSGs in B. belcheri

Results showed that 16.52% (5,020 genes) of all B.belcheri genes exhibited time-specific expression (τ > 0.85) in this study. Moreover, population distributions of maximal JS scores were similar with those of previous studies (Fig. 2A),, suggesting our analysis was reliable. As a result, the 48 hpi sample harbored the largest number of TSGs (1,624 genes), followed by 0 (1,330 genes), 24 (911 genes), 12 (670 genes) and 6 (485 genes) hpi (Fig. 2B).

Figure 2.

A. Time specificity of genes. Distribution of maximal tissue specificity scores calculated for each time-specific gene (τ > 0.85, TSGs) across all 5 time points. B. Summary of the number of TSGs at each time point.

GO and KEGG enrichment analysis of TSGs

To explore the functions of the TSGs at each time point, we performed GO enrichment analysis for each TSG set. TSGs were assigned to the largest number of GO terms (174 terms) at 48 hpi, followed by 24 (76 terms), 0 (69 terms), 12 (44 terms), and 6 (8 terms) hpi (Supplementary Table S6). Based on these enriched GO terms, we obtained 14, 0, 3, 45, and 85 time-specific categories for 0, 6, 12, 24, and 48 hpi, respectively (Supplementary Figure S1). Meanwhile, some useful time-specific GO information was obtained at each time points (with exception of 6 hpi without significantly enriched terms) (Supplementary Table S7). For example, for 0 h, time-specific GO terms were involved in reproductive development and biomacromolecule biogenesis and organization, including reproductive system development (GO: 0061458), reproductive structure development (GO: 0048608), cellular component biogenesis (GO: 0044085), protein–DNA complex subunit organization (GO: 0071824), and nucleosome organization (GO: 0034728). There were three time-specific GO terms associated with regulation and signal transduction at 12 hpi, including transferase activity, transferring sulfur-containing groups (GO: 0016782), biologic regulation (GO: 0065007), and signal transduction activity (GO: 0004871). For 24 hpi, we found that most time-specific GO terms were related to cell cycle and death, DNA replication and repair, and chromosome stability, including cell cycle phase (GO: 0022403), programmed cell death (GO: 0012501), DNA replication (GO: 0006260), response to DNA damage stimulus (GO: 0006974), and chromosome organization (GO: 0051276). Time-specific GO terms at 48 hpi were involved primarily in biosynthetic, metabolic, and catabolic processes, such as regulation of biosynthetic process (GO: 0009889), cellular macromolecule metabolic process (GO: 0044260), and heterocycle catabolic process (GO: 0046700).

KEGG enrichment analysis of the TSGs at each time point was conducted to further uncover time-specific pathways. TSGs were significantly enriched in the largest number of pathways at 6 hpi (13 terms), followed by 24 (10 terms), 48 (9 terms), 12 (5 terms), and 0 (3 terms) hpi (Supplementary Table S8). Furthermore, some enriched pathways were shared by multiple time points, such as neuroactive ligand-receptor interaction (ko04080) (shared by 6, 12, 24, 48 hpi); microRNAs in cancer (ko05206) and cell adhesion molecules (CAMs) (ko04514) (shared by 6, 12, 48 hpi); and ABC transporters (ko02010), Fanconi anemia pathway (ko03460), and primary immunodeficiency (ko05340) (shared by 6 and 12, 24, or 48 hpi, respectively). In particular, 3, 7, 1, 8, and 5 pathways were specific to 0, 6, 12, 24, and 48 hpi, respectively (Table 3). The 3 time-specific pathways at 0 hpi were vasopressin-regulated water reabsorption (ko04962), progesterone-mediated oocyte maturation (ko04914), and phospholipase D signaling pathway (ko04072). Interestingly, although time-specific pathways were primarily related to immunity and defense responses at 6, 12, 24, and 48 hpi, at each time point post-LPS treatment, the specific pathways participated in different immune steps. First, at 6 hpi, specific pathways were involved in extracellular process, immune cell surface signaling, and phagocytosis, including ECM-receptor interaction (ko04512), T cell receptor signaling pathway (ko04660), and Fc gamma R-mediated phagocytosis (ko04666). At 12 hpi, only one specific pathway was identified, complement and coagulation cascades (ko04610). Next, at 24 hpi, most specific pathways played key roles in cell cycle and DNA replication and repair, including cell cycle (ko04110), DNA replication (ko03030), and mismatch repair (ko03430). Finally, at 48 hpi, we detected time-specific pathways involved in organics synthesis and metabolism, including glycosaminoglycan degradation (ko00531) and thyroid hormone synthesis (ko04918).

Table 3.

Time-specific pathways enriched among TSGs at all 5 time points. TSGs, time-specific genes at each time point; hpi, hours post-injection.

Pathway	Pathway ID	FDR
0 h (control)
Vasopressin-regulated water reabsorption	ko04962	1.08E-09
Progesterone-mediated oocyte maturation	ko04914	1.96E-03
Phospholipase D signaling pathway	ko04072	2.93E-03
6 hpi
Prion diseases	ko05020	1.571E-10
Dorso-ventral axis formation	ko04320	1.174E-06
Notch signaling pathway	ko04330	1.266E-06
Thyroid hormone signaling pathway	ko04919	3.719E-05
T cell receptor signaling pathway	ko04660	0.0007369
Fc gamma R-mediated phagocytosis	ko04666	0.0060352
ECM-receptor interaction	ko04512	0.0464513
12 hpi
Complement and coagulation cascades	ko04610	1.25E-02
24 hpi
DNA replication	ko03030	2.105E-15
Cell cycle	ko04110	5.712E-08
Mismatch repair	ko03430	0.0004709
Pyrimidine metabolism	ko00240	0.0011964
Homologous recombination	ko03440	0.0032069
Oocyte meiosis	ko04114	0.0057862
Nucleotide excision repair	ko03420	0.0114305
Systemic lupus erythematosus	ko05322	0.049154
48 hpi
Olfactory transduction	ko04740	0.0037022
Glycosaminoglycan degradation	ko00531	0.004054
Thyroid hormone synthesis	ko04918	0.0116272

Validation of gene expression profiles using qRT-PCR

To validate the accuracy of the expression profiles obtained from RNA-seq, we measured the expression of 12 randomly selected genes by qRT-PCR across all time points using the same samples originally used in RNA-seq. A linear regression analysis of the fold-changes in the gene expression ratios between the RNA-seq and qRT-PCR data showed a significant correlation for each CG (Pearson's correlation coefficients, r = 0.88 for CG1, r = 0.89 for CG2, r = 0.90 for CG3, r = 0.89 for CG4 (Fig. 3).

Figure 3.

Confirmation of the RNA-seq results via qRT-PCR of 12 randomly selected DEGs in the CG1 (A), CG2 (B), CG3 (C), and CG4 (D). Each point represents the fold-change in expression relative to 0 h (control), as measured by RNA-seq (x-axis) and qRT-PCR (y-axis).

Discussion

Libraries constructed by removing rRNA have shown higher rates of intron-mapped reads than those of polyA libraries in previous studies. RNA-seq of total RNA requires higher read depths than that needed for rRNA-depleted libraries to achieve equal coverage of exons. Here, we obtained ∼3.3Gb clean data that could be mapped to the CDS regions of the B. belcheri genome (∼32Mb CDS sequences); thus, read coverage showed a high depth (∼100 ×). Moreover, due to the availability of whole-genome sequence information for B. belcheri, we were able to efficiently extract reads that mapped to exons to avoid problems caused by intron-mapped reads, and we removed intron-mapped reads for identification of novel genes and redefinition of transcript structure in the future. With the exception of reads mapped to CDS regions, the other reads obtained here primarily derived from non-coding transcripts, such as long non-coding RNAs and small nucleolar RNAs. These non-coding resources may be further used in the identification of immune-related non-coding RNAs in amphioxus, as reported in fish and shellfish., Therefore, strand-specific RNA-seq using rRNA-depletion represents a sensitive and unbiased approach for providing complete genetic information.,, To better understand the mechanisms of the molecular response to LPS challenge and to uncover key genes, functional classes, and complex pathways involved in immune response, we analyzed genes in B. belcheri that were differentially expressed between different time points (including controls vs. different time points after LPS injection, CG1–4; comparison between different time points after LPS injection, CG5–6) using RNA-seq technology. Subsequently, we also identified a set of TSGs for each time point, allowing us to obtain time-specific functional clusters and pathways to further understand the specific molecular bases of the immune response in amphioxus at different stages. The qRT-PCR results of 12 randomly selected DEGs showed expression patterns consistent with the RNA-seq experiment, indicating that our sequencing results were reliable.

NGFR, a known nerve growth factor receptor and a member of the tumor necrosis factor receptor (TNFR) superfamily, stimulates neuronal cells to survive and differentiate. Neurotrophins might signal a cell to die via apoptosis by activating NGFR, considered a marker of perineural invasion in malignant melanomas in mammals. LPS-induced neuro-inflammation was confirmed by experiments in a rat model, and we also detected several positive sites in the dorsal nerve area of adult B. belcheri challenged by LPS (24 hpi) using immunohistochemistry (data not shown). Here, widely different expression patterns of NGFR were found across 4 CGs (CGs 2–5). Notch, a neurogenic locus protein, controls oligodendrocyte maturation and affects the activation of microglia under inflammatory conditions. In this study, Notch was upregulated at 12 hpi compared with control levels, and gene encoding neurogenic locus notch homolog protein 1 (Notch1a) also showed different expression in CG6. SSPO, an evolutionarily conserved gene in the central nervous systems of chordates, has a neuroprotective effect by promoting neuronal differentiation and growth in spinal cord injuries., Our results showed that the downregulation of SSPO could be triggered by LPS infection in amphioxus. Therefore, we propose that the differential expression of these key neuro-related genes may be involved in neuro-inflammation, implying that a rudimentary, active neuro-immune system probably existed in basal chordates and acted as an ancient immune system in the ancestor of vertebrates.

Many of the top 20 DEGs in each CG were related to the immune response of amphioxus to LPS, such as genes encoding complement proteins, coagulation factors, E3 ubiquitin-protein ligase, mucin, interferon-induced protein, Nek7, MFN2, TMED9, and adhesion GPCRs.,,, This observation suggests that these genes play a key role in immunity and can be considered molecular indicators in immune research of cephalochordates under LPS stress. Interestingly, West et al. found an important function of mitochondria in innate immunity. MFN2, a molecule embedded in the mitochondrial outer membrane, participates in RIG-I-like receptor signaling (a pathway that was significantly enriched in our study, ko04622) in antiviral immunity. In addition to its antiviral function, we noted that MFN2 displayed an acute expression response to LPS here. Mitochondrial fission is key organelle response to immune challenge. Mitochondrial fission factor (MFF) promoted this fission process by interacting with dynamin-related protein 1 (Drp1) in mammalian cells, and we found different expression of MFF in CG6. Furthermore, adhesion GPCRs could affect mitochondrial function and have been shown to play an important pathological role in zebrafish via knockdown of adhesion GPCRs. This finding suggests that mitochondrial immunity may have been important in the ancestor of vertebrates and that this type of immunity may have already appeared by ∼520 million years ago.

By comparing the results of GO enrichment of DEGs among CGs, we found that genes associated with molecule binding, enzyme and protein activity, metabolic process, cell death, response to stimulus, and cell region were implicated at all immune time points within 48 h after treatment (CG1–4), and genes involving response to stimulus, defense response, metabolic process, cellular process, and response to stress showed significantly different expression between time points during immune response (CG5–6), suggesting that these genes may be essential in the immune response of amphioxus to LPS challenge. Meanwhile, enriched biologic process terms that directly participated in immune signaling pathways were detected in at least one CG in CG1–4 (excluding all CG1–4), such as defense response (GO: 0006952) and inflammatory response (GO: 0006954). Similarly, we also detected directly immune-related biologic process terms in at least one CG in CG5–6, such as response to wounding (GO: 0009611) and inflammatory response (GO: 0006954). The expression patterns of these related genes were the same as those in vertebrates during the acute immune response,, implying a certain evolutionary conservation in immunity between vertebrates and amphioxus. In addition, apoptosis plays a critical role in the immune response to remove injured or necrotic cells or cells with potential danger. Here, GO terms involved in apoptosis were enriched in each of CG1–5, suggesting that apoptosis-related genes play protective roles in host defense against LPS infection in B. belcheri, as reported in the hemocytes of Litopenaeus vannamei exposed to nitrite stress. GS1 genes were significantly enriched in GO terms involved in extracellular region, suggesting that immune function against LPS infection in amphioxus depends not only on the complexity and diversity of sequences but also on the high expression plasticity of extracellular genes.

The results of KEGG enrichment analysis of DEGs revealed the well-known signaling pathways of innate sensors, including the RIG-I-like receptor signaling pathway (ko04622) and NOD-like receptor signaling pathway (ko04621). Some adaptors in these 2 pathways trigger the NF-κB signaling pathway (ko04064) to directly activate NF-κB. In addition, other many pathways also have important functions in innate immunity, such as cytokine-cytokine receptor interaction (ko04060), the TGF-β signaling pathway (ko04350), cell cycle (ko04110), endocytosis (ko04144), antigen processing and presentation (ko04612), Salmonella infection (ko05132), apoptosis (ko04210), and metabolic pathways (ko01100)., Due to the absence of an adaptive immune system, amphioxus strongly relies on its innate immune system to protect the host from pathogens., Indeed, many of the best-known pathways associated with innate immunity were enriched in this study. Notably, some genes were involved in pathways in cancer (ko05200), bladder cancer (ko05219) and microRNAs in cancer (ko05206), as well as Huang et al. annotated some genes involved in cancer in the B. belcheri genome. Our results indicate that the origin of functional cancer genes may be extremely ancient within chordates. Pathways associated with infection and invasion of microbes and extracellular innate immunity processes were specifically identified in CG1, while pathways involved in intracellular immune signaling and immune cell function, nucleic acid heredity and repair, and metabolism and cancer were specifically observed in CG2, CG3, and CG4, respectively. Additionally, for two CGs during immune response, we found very few pathways involving innate immune signaling in the list of CG-specific terms between CG5 and CG6, and many pathways related to cancer and disease, metabolism, cell cycle were overrepresented. These specific KEGG terms identified in each CG indicate that the dynamic changes in expression of gene involving immunity, metabolism and disease primarily occur in response to LPS infection in amphioxus.

Gene expression exhibits spatiotemporal specificity in normal tissues and organisms challenged by LPS, as shown in previous studies., Here, we identified many TSGs, which further extended our knowledge of the ways in which gene expression contributes to amphioxus immunity through high spatiotemporal plasticity at different time points. The present study showed that some TSGs enriched in certain GO and KEGG terms were shared across multiple time points. Furthermore, we obtained some time-specific terms for each time point. GO and KEGG terms involved in reproductive development and biomacromolecule biogenesis and organization were specifically detected at 0 h. Previous studies have shown that organismal LPS exposure can alter reproductive development in rat., In addition, DEGs were significantly enriched in GO terms involved in biomacromolecule biogenesis in L. vannamei under nitrite stress. LPS immune challenge thus may affect the reproduction and development of amphioxus. While no time-specific GO terms were detected at 6 hpi, we identified some well-known pathways that participate in the function of cell surface receptors and phagocytosis. At 12 hpi, time-specific GO terms were related to transferase activity, signal transduction, and biologic regulation, and these functions play key roles in “complement and coagulation cascades” pathway, indicating that immune signal transduction and cascades may be extremely active at this stage. At 24 hpi, almost all specific GO and KEGG terms were related to DNA repair and stability and the cell cycle, suggesting that TSGs associated with response to DNA damage, replication, and stability may be important in amphioxus at this time point under LPS challenge., Finally, time-specific GO terms and pathways at 48 hpi were primarily related to biosynthetic, metabolic, and catabolic processes. It has been widely reported that the genes belonging to these GO categories participate in immune response,, as various related catabolic and metabolic genes are required to meet metabolic demands for cell growth, proliferation, and differentiation in cancer and immunity.,

In summary, to investigate the immune responses of B. belcheri to LPS, we performed a gene expression analysis at the transcriptional level of 5 rRNA-depleted libraries using strand-specific RNA-seq. The potential roles of DEGs identified here were analyzed using alignment and GO and KEGG enrichment analysis. The results demonstrate that the acute immune response to LPS challenge in amphioxus shares certain similarities to that in vertebrates. In addition, the identification and GO and KEGG enrichment analysis of TSGs demonstrated a dynamic time-specific process of immunity. Notable, due to a large unknown in immune response, the functional relevance of some genes in the present study could not be explained. Besides, because of complexity in immune system, DEGs, GO terms and pathways obtained in this study includes not only members directly participated in immune response, but also indirect those (such as typical genes and terms in metabolism, development, regulation, and so on). To our knowledge, this is the first study to explore the transcriptome dynamics at multiple time points at a genome-wide level after LPS injection in amphioxus. Our experiment not only provides useful information to further our understanding of the complex and diverse immune systems of cephalochordates, but our data may also aid in revealing a deeper understanding of the origin of the vertebrate immune system in the future.