Transcriptional and in silico analyses of MIF cytokine and TLR signalling interplay in the LPS inflammatory response of Ciona robusta

The close phylogenetic relationship between Ciona robusta and vertebrates makes it a powerful model for studying innate immunity and the evolution of immune genes. To elucidate the nature and dynamics of the immune response, the molecular mechanisms by which bacterial infection is detected and translated into inflammation and how potential pattern recognition receptors (PRRs) are involved in pathogen recognition in tunicate C. robusta (formerly known as Ciona intestinalis), we applied an approach combining bacterial infections, next-generation sequencing, qRT-PCR, bioinformatics and in silico analyses (criteria of a p-value < 0.05 and FDR < 0.05). A STRING analysis indicated a functional link between components of the Tlr/MyD88-dependent signalling pathway (Tlr2, MyD88, and Irak4) and components of the Nf-κB signalling pathway (Nf-κB, IκBα, and Ikkα) (p-value < 0.05, FDR < 0.05). A qRT-PCR analysis of immune genes selected from transcriptome data revealed Mif as more frequently expressed in the inflammatory response than inflammation mediator or effector molecules (e.g., Il-17s, Tnf-α, Tgf-β, Mmp9, Tlrs, MyD88, Irak4, Nf-κB, and galectins), suggesting close interplay between Mif cytokines and Nf-κB signalling pathway components in the biphasic activation of the inflammatory response. An in silico analyses of the 3′-UTR of Tlr2, MyD88, IκBα, Ikk, and Nf-κB transcripts showed the presence of GAIT elements, which are known to play key roles in the regulation of immune gene-specific translation in humans. These findings provide a new level of understanding of the mechanisms involved in the regulation of the C. robusta inflammatory response induced by LPS and suggest that in C. robusta, as in humans, a complex transcriptional and post-transcriptional control mechanism is involved in the regulation of several inflammatory genes.


Transcriptomic analysis of Ciona robusta highlights the effects of LPS on transcripts involved in inflammation.
To profile the C. robusta inflammatory response under LPS exposure in vivo, we investigated the differential gene expression in its immunocompetent organ (pharynx) under physiological and challenged conditions (4-h LPS exposure) using NGS transcriptomic analysis. Transcriptomic analysis under physiological conditions revealed 16,504 detected sequences, including protein-coding transcripts, non-proteincoding transcripts, isoforms, novel transcripts and protein-RNA interaction sites (data not shown). All 16,504 genes produced by sequencing were annotated using the Ensembl database (ensemble.org), and for each of them, log-fold changes, p-values and adjusted p-values were calculated (data not shown). With the log-fold change set at 1.5, 1,227 genes were found to be differentially expressed (p-value threshold < 0.05, and adjusted p-value < 0.05). A comparison between the basal gene expression of three specimens showed physiological interindividual variability (Fig. 1a, compared with controls, columns 4, 5 and 6), but a large set of responsive genes (1,227) were found to be reproducibly regulated after 4 h of LPS exposure in vivo (p-value ≥ 1.5-fold higher or lower than controls) (Fig. 1a, columns 1, 2 and 3). Notably, the transcript levels of 447 target genes were found to be downregulated in response to LPS, whereas 780 genes were upregulated (Fig. 1b). For functional annotations, unigenes (1,227) were aligned to Gene Ontology (GO) terms based on the Protein ANalysis THrough Evolutionary Relationships (PANTHER) (pantherdb.org) classification system for gene ontology annotations, classified into three subcategories: (i) GO molecular functions (Fig. 2a); (ii) GO cellular components (Fig. 2b); (iii) GO biological processes (Fig. 2c).
Scientific RepoRtS | (2020) 10:11339 | https://doi.org/10.1038/s41598-020-68339-x www.nature.com/scientificreports/ Genes responsive to LPS were found to be predominantly involved in molecular functions leading to gene expression and regulation (e.g., transcription regulatory region DNA binding, DNA binding, RNA binding), transporter activities (e.g., amino acid transmembrane transporter activity) and monosaccharide-binding activities (Fig. 2a). Gene expression programmes, which establish and maintain specific cell states, are controlled by transcription factors, cofactors, chromatin regulators and noncoding RNAs that interact with regulatory regions. Notably, genes involved in DNA-binding transcription activity were found to be significantly downregulated (3.5-fold decrease) upon LPS exposure, whereas those involved in amino acid transmembrane transporter and monosaccharide-binding activities were found to be upregulated (1.5-and 2.5-fold increase) (Fig. 2a). Concordantly, genes specific to cellular components were expressed at lower levels than the respective control group (e.g., processes involving the Golgi apparatus, nucleoplasm, chromosome, chromatin) (Fig. 2b). Within the biological processes category, the RNA-processing group was found to be strongly reduced (3.8-fold decrease) (Fig. 2c). On the other hand, the response to xenobiotic stimuli and drug processes was found to be moderately increased, whereas organic acid metabolic processes and the immune system processes were found to be slightly increased. The pathway analysis of differentially expressed unigenes with PANTHER revealed different pathways grouped into different categories (e.g., inflammation and immune response, growth factor/growth factor receptor-mediated pathways and intracellular signalling activation, apoptosis and p53 regulation, and metabolism). The first categories listed herein, according to GO analysis results, suggest a functional link among LPS inflammatory processes, cytokines, Tgf-β and Toll receptor signalling pathways (Table 1). STRING-protein-protein interaction networks functional analysis (string-db.org) was used to predict interactions among Tlr/MyD88-dependent signalling pathways and components of Nf-κB signalling by generating network clustered that are visualised according to the k-means algorithm (clustering analysis). The STRING analysis showed a greater functional link between Tlr2 and MyD88 (interaction score of 0.946) than was shown for Tlr1 and MyD88 (interaction score of 0.795) and suggested a strong functional link between Tlrs/Myd88 and Nf-κB, which was confirmed by the identification of two distinct functional clusters: the first included Tlr1, Tlr2, MyD88 and Irak4, and the second included Nf-κB, IκB, IKKα and Rel 1 (Fig. 3). Blue points represent DEGs (up-and downregulated); red points represent non-DEGs. The y-axis indicates − log10 (p-value), and the x-axis indicates log2 (fold change). Genes significantly (p-value < 0.05) up-or downregulated (log2-fold change < 1.5 or log2-fold change > 1.5) are displayed in blue.  19 . Tlr proteins present the "hybrid" functionality of mammalian TLRs: Tlr1 (NP_001159599.2) shows 23% identity and 42% sequence similarity with the structurally and functionally well-characterised human TLR5, whereas Tlr2 (NP_001159600.1) shows 27% identity and 43% sequence similarity with human TLR8. In silico analysis was performed using the RegRNA2.0 database (https ://regrn a.mbc.nctu.edu.tw/html/ predi ction .html) to identify functional RNA motifs and elements in the 3′ (UTR) involved in Tlr mRNA posttranscriptional regulation. We found that the Tlr1 3′-UTR contains a γ-interferon activated inhibitor of translation (GAIT) element, a Mos polyadenylation response element (MOS-PRE) (Fig. 4a). GAIT elements are implicated in several immune-related mRNAs, showing important roles in gene-specific translation control in innate immunity; MOS-PRE is involved in both the translational repression of mRNA stored in immature oocytes and translational induction in maturing oocytes in X. laevis; in contrast, the Tlr2 3′-UTR contains a GAIT element, a Musashi-binding element (MBE), previously identified as being involved in mRNA translational control during cell cycle progression, and an AU rich element (ARE) (Fig. 4a). The C. robusta MyD88 gene has not yet been cloned, and to isolate full-length mRNA, we used a RACE 5′ and 3′ strategy. MyD88 mRNA presents a 5′-UTR of 57 bp, an ORF of 1,210 bp and a 3′-UTR of 593 bp, which contains an MBE and GAIT element (Fig. 4a). It encodes a 403 amino acid polypeptide with a predicted molecular size of 45.169 kDa and a pI of 5.15 ( Supplementary Fig. 1).
In silico analysis of the MyD88 amino acid sequence performed using both the Delta-Blast algorithm and the Simple Modular Architecture Research Tool (SMART) (smart.embl-heidelberg.de) jointly predicted the presence of specific structural domains: a DEATH domain (from 105 to 196 aa) and a Toll/interleukin-1 receptor (TIR) homology domain (from 265 to 399 aa) ( Supplementary Fig. 1).
The C. robusta myD88-deduced amino acid sequences were examined in GenBank through Basic Local Alignment Search Tool (BLAST) (blast.ncbi.nlm.nih.gov) analysis and showed 28% identity and 48% sequence similarity with human MyD88. Molecular Evolutionary Genetics Analysis (MEGA X) software was used to build a phylogenetic tree by the neighbour-joining method (Fig. 4b). Phylogenetic analyses of vertebrate and invertebrate MyD88 proteins supported the idea of a conserved evolution from a common MyD88 ancestral gene among invertebrates, protochordates and vertebrates. The functionally and structurally conserved MyD88 motifs highlighted by CLC workbench 6.4 are shown in the supplementary material ( Supplementary Fig. 2).
STRING-protein-protein interaction networks functional analysis suggested that in C. robusta, the major signalling pathways involved in the activation of Nf-κB were canonical pathways including Ikkα and IκBα proteins. In C. robusta, Nf-κB (NP_001071772.1) shows 53.36% identity with the structurally and functionally characterised NF-κB1 in humans and 41.49% identity with the Dorsal protein of the Rel subfamily present in D.   Figure 5a shows the homodimer molecular model of C. robusta Nf-κB constructed according to the homology-modelling process performed on the basis of the known crystal structure of H. sapiens NF-κB1 (1svc.1.B) (Fig. 5a), generated from the super-imposition of the 2-326 residue sequence, shows 59.81% identity with the template. The phylogenetic tree representing the vertebrate and invertebrate NF-κB family members suggests that C. robusta Nf-κB is an orthologue of vertebrate NF-κB and that, in vertebrates, a single gene was duplicated into two (NF-κB1 and NF-κB2) after the divergence of the tunicate and vertebrate lineages (Fig. 5b).
In humans, the activation of the NF-κB "canonical pathway" activates the protein kinase complex IKK 30,31 ; in C. robusta, only one member of the Ikkα (XP_002125567.1) and IκBα (NP_001071739.1) families was identified. C. robusta Ikkα shows 31.53% identity with human IKKα and 32.62% identity with D. melanogaster Ikk. SMART analysis identified one PKinase-Tyr domain in C. robusta Ikkα (17-261 aa), H. sapiens IKKα and D. melanogaster Ikk (Fig. 5c). Moreover, H. sapiens IKKα presents a PHB and an IKKα/NEMO-binding domain, both of which are involved in the formation of the IKK complex and are not present in C. robusta Ikkα or D. melanogaster Ikk. Figure 5c shows the homodimer molecular model of C. robusta Ikkα based on the known crystal structure of the Mus musculus serine/threonine-protein kinase (4jlc.1.A) generated from the super-imposition of the 12-79    Figure 5e shows the homodimer molecular model of C. robusta IκBα based on the known crystal structure of H. sapiens Tankyrase-1 (3utm.1.A) generated from the super-imposition of the 123-367 residue sequence that shows 23.65% identity with the IkBα protein template. The phylogenetic analysis results support the supposition that C. robusta IκBα (Fig. 5f) is an orthologue of vertebrate IκBα genes and has amino-terminal signal-responsive regions containing conserved serine residues (e.g., Ser 84 and Ser 88 ) that are phosphorylated in an 83-DSPGXXSP-89 amino acid sequence. An in silico analysis was performed based on RegRNA2.0 databases to identify functional RNA motifs and elements in the 3′-UTR of the Nf-κB, Ikkα and IκBα genes with identified GAIT elements (data not shown).

Analyses of the expression of Tlr/MyD88-dependent and Nf-κB canonical signalling pathways and immune-related genes under LPS exposure. Analyses of the time-course expression of Tlr/
MyD88-dependent and Nf-κB canonical signalling pathway immune-related genes in the pharynx inflammatory response induced by LPS in C. robusta were performed at time points from 0 to 72 h post-LPS challenge by qRT-PCR (Fig. 6). The heatmap shows that a large portion of the transcripts was significantly modulated in response to LPS during the 72-h period of LPS exposure (p-value < 0.05). Based on the expression patterns of the transcripts, two major clusters were highlighted: the first includes proinflammatory cytokines and components of MyD88/Nf-κB signalling, and the second comprises the Mif2 cytokine IκBα and a few inflammation effectors (e.g., lectins, lysozymes, and Po

Discussion
In this study, we found that the innate immune signalling pathway activated by LPS in C. robusta is evolutionarily conserved and involves both the Tlr and Mif signalling pathways, findings in agreement with reports highlighting the key immune signalling pathways activated against invading microbial pathogens and other potential threats to a mammalian host 46 . In general, these findings are not surprising, since this invertebrate occupies a key phylogenetic position in chordate evolution and is considered to be a member of the sister group of vertebrates 5-7 . In our model, the Mif signalling pathway seems to be activated as the first line of defence against the bacterial www.nature.com/scientificreports/ endotoxin LPS and appears to be a major regulator of Il-17s, Tgf-β and Mmp9 transcription, as these levels were found to be upregulated after 1 h of exposure. In contrast, the Nf-κB signalling pathway seems to be recruited a few hours later to activate Tgf-β through the involvement of Tlr2, MyD88 and Irak4. In humans, LPS induces TLR4 activation through binding with three proteins (the lipopolysaccharide-binding protein LBP, the cluster differentiation antigen CD14, and the myeloid differentiation protein receptor MD-2) 47 . In humans, MD-2 is a component of the MD-2-related lipid-recognition (ML) superfamily, which contains a large set of genes encoding proteins such as MD-1, MD-2, and Niemann-Pick type C2 (NPC2) protein. The ascidian C. robusta possesses two Tlr-related genes presenting "hybrid" biological and immunological functionality of mammalian TLRs 19 , but LPS significantly activates only one of them, Tlr2; on the other hand, no MD-2 orthologue was found in the genome of C. robusta. The subtractive hybridization strategy allowed the identification of Npc2 mRNA in C. robusta, an homologue of ML superfamily components, preferentially expressed in haemocytes inside the vessel lumen of the pharynx 48 . This finding supports the hypothesis that the recognition of LPS by TLR4 through MD-2 binding may have been gained early during the evolution of vertebrates and that C. robusta may respond to LPS through a complex with alternative LPS sensors (i.e., the Ncp2 protein) developed before the differentiation of the MD-2 protein in vertebrates. LPS activation of Tlr2 likely regulates the Rel/Nf-κB signal transduction pathway, and we hypothesise that Tlr1 in is involved in the signalling that activates Mif1 in C. robusta. Notably, the molecular mechanisms underlying Tlr ligand recognition and signal transduction are distinct from those of vertebrate TLRs 19 . The qRT-PCR results suggest an intriguing possibility of biphasic action of the Nf-κB and Mif signalling pathways during the 72 h of exposure to LPS, suggesting their different roles in regulating immune cell behaviour, respectively, with a pro-inflammatory action period (from 0 to 12 h) and an anti-inflammatory action period at which time homeostasis is restored (from 12 to 72 h).   Recently, a few studies have indicated that, in humans, intracellular MIF is involved in the induction of NF-κB activity 49 , reinforcing the idea of a close relationship between Mifs, Tlrs and Nf-κB in C. robusta. In agreement, Roger et al. 50 found that MIF upregulates Tlr4 expression in mouse macrophages and activates the MAPK and NF-κB signalling pathways in immortalised cell lines and mouse macrophages 46 . In humans, the upregulation of MIF during inflammation and under other stress conditions is primarily regulated at the level of MIF release rather than transcriptional induction 51 , which explains its rapid activation. In C. robusta, Mif2 is not regulated by LPS, and its 3′-UTR mRNA has a CPE element, confirming that it can be stored in the cytoplasm for rapid activation by a cytoplasmic polyadenylation mechanism; instead, Mif1 is upregulated after LPS stimulation and does not have a 3′-UTR CPE 45 . Moreover, the in silico analysis of the 3′-UTR of Tlr1 revealed a MOS-PRE that includes the 5′ sequence of the CPE 52 . Both the CPE and MOS-PRE enable translational repression or activation in response to specific stimuli; therefore, we can hypothesise that, following LPS stimulation, Mif1 and Tlr1 transcripts gain poly(A) tails, which activate their translation and allow the rapid and abundant expression of proteins for their prompt response.
In humans, MIF can regulate TNF and other cytokines by affecting the expression of TLR4, p53, ERK, c-Jun-N-terminal kinase, p38, and MAPK phosphatase-1 53 . On the other hand, the stimulation of dendritic cells with TLR4 ligands elicits high levels of MIF production, and secreted MIF acts as an autocrine/paracrine enhancer of TNF production. Our findings support the notion that there is interplay between Mifs and Nf-κB pathway components in C. robusta, meaning that they have an effect on each other. Specifically, we speculate that to actively maintain the immune response and increase its efficiency, Mif2 and Tlr1 are involved in a rapid response (within a few minutes) and stimulate the transcription of Mmp9 and the pro-inflammatory cytokines Mif1, Tgf-β, Il-17 and Tnf-α, which, in turn, activate Nf-κB through Tlr2, eliciting renewed Mif2 activation until immune response resolution. In silico analyses of the cis-regulatory elements of the 3′-UTR in C. robusta based on RegRNA 2.0 revealed a GAIT element in Tlr2, MyD88, IκB, Ikkα, and Nf-κB. In humans, GAIT elements are implicated in several immune-related mRNAs, showing an important role in gene-specific translation control in innate immunity 54 . In C. robusta, GAIT elements have already been identified in Cap 55 and Tnf-α 56

and in
Mif mRNAs 45 . Notably, in C. robusta, Mifs also seem to regulate Mmp9, which has several important immune functions, such as extracellular matrix degradation and Tnf-α and Tgf-β activation 57 . In humans, TGF-β is a crucial enforcer of immune homeostasis and tolerance, inhibiting the expansion and function of many immune system components. It can be viewed as immunosuppressive and plays important roles both during the initial phase of immune injury and during the later remodelling phase 58 . In agreement with these findings, Tgf-β was www.nature.com/scientificreports/ overexpressed in C. robusta after 4 h of LPS exposure and later, from 12 to 72 h, thus regulating the resolution of the inflammatory response. Finally, we propose a novel schematic model representing the putative interplay between the Mif and Nf-κB pathway components (Fig. 7). To the best of our knowledge, this is the first study to show evidence of interplay between the Mif and Nf-κB pathways using a marine invertebrate model through a wide-ranging approach. Further studies in this direction are needed to cover knowledge gaps regarding the hierarchically organised set of molecular, cellular and organismal networks involved in universal immune interactions with pathogens and subsequent intracellular signal transduction.

Methods
Tunicates and LPS injection. Molecular studies have led to the hypothesis that Ciona intestinalis constitutes a compilation of species rather than a single species [10][11][12][13] . Accordingly, our model organism, originally from the Mediterranean Sea and formerly classified as C. intestinalis, corresponds to C. robusta [10][11][12][13] .
C. robusta specimens were collected from Sciacca Harbour (Sicily, Italy) and were acclimatised and maintained under controlled temperature conditions (15 °C) in tanks supplied with flow-through oxygenated seawater. The animals were fed every 2 days with Coraliquid marine invertebrate food (Sera Heinsberg, Germany). An LPS solution (Escherichia coli 055:B5, LPS, SIGMA-ALDRICH, Germany) was prepared in sterile salt medium (12 mM CaCl 2 , 11 mM KCl, 26 mM MgCl 2 , 43 mM Tris HCl, 0.4 M NaCl, pH 8.0). One hundred microlitres of the LPS-containing suspension was injected into the tunic matrix surrounding the pharynx wall (median body region) at a final LPS concentration of 100 μg. C. robusta not exposed to LPS (naïve) were used as controls. Fragments of pharynx tissue (200 mg) explanted at various times (from 1 to 72 h) were immediately soaked in RNAlater tissue collection solution (AMBION, Austin, TX) and stored at − 80 °C. Total RNA extraction was performed using an RNAqueous-Midi kit purification system (AMBION, Austin, TX).
Transcriptomics. The RNA purity and quality of total RNA extracted from the pharynx of C. robusta that were naïve (N = 3) and that were exposed to LPS for 4 h (N = 3) were examined by NanoDrop and Agilent RNA 6,000 Nano kits on an Agilent 2,100 Bioanalyzer (AGILENT, USA), respectively. High-quality RNA samples (A260/A280 = 1.9-2.1, RIN ≥ 7) were used for cDNA library construction. RNA sequencing (RNA-Seq) was performed by BMR Genomics (Padua, Italy) on an Illumina platform in a single-end format 75 bp (1 × 75 bp) containing ~ 40 million ± 10% of reads/sample. An analysis of differentially expressed genes was performed by BMR Genomics (Padua, Italy) using edgeR software, which estimates the negative binomial variance parameter globally across all genes. All data from the edgeR analysis were normalised by setting the false discovery rate (FDR) to ≤ 0.05 and the absolute value of the log fold change (logFC) to 1.5. A heatmap and a volcano plot were generated to visualise the results indicating the genes differentially expressed between the exposed samples and controls. The heatmap was produced using a heatmapper tool (https ://www.heatm apper .ca). Complete linkage clustering was applied, and Spearman's rank correlation was used as the method of distance measurement. In the volcano plot, the log 2 CPM, which represents the log base 2 of the count per million (CPM) for each transcript, was calculated by the CPM function. Negative log values represent transcripts presenting a CPM lower than 1, and positive log values represent those presenting a CPM ≥ 1.

Gene functional enrichment analysis and immune network construction. The Clusters of
Orthologous Genes (COG) database (geneontology.org) was used to perform a Gene Ontology (GO) enrichment analysis of a set of C. robusta transcripts produced by NGS. The ontology analysis was divided into three subcategories: (i) molecular functions (MF); (ii) cellular components (CC); (iii) and biological pathways (BP). The Protein Analysis Through Evolutionary Relationships (PANTHER GO-slim analysis tool) System connected to the COG database was utilised to expand the annotation data of the gene and protein families obtained from GO. A GO analysis was performed by selecting the PANTHER "GO-slim" analysis mode, as it is more reliable and accurate than the GO annotation mode "GO-complete". Different parameters used to measure the functional role of the list of genes were taken into account and reported as fold enrichment. In addition to p-value evaluation (p-value < 0.05), false discovery rate (FDR < 0.05) control procedures were also considered. The COG database was also used to perform pathway analysis of the NGS data. The Search Tool for the Retrieval of Interacting Genes/Protein (STRING) database, which allows the visualisation of complex networks (through clustering analysis), was used to retrieve the predicted interactions for the identified immune-related proteins. Levels of significance were set at a p-value < 0.05; the log2-fold change was </> 1.5. Significant data (proteins significantly modulated) were clustered according to the k-means algorithm, an unsupervised clustering algorithm based on an adjacency matrix.
Cloning and bioinformatics analyses of MyD88. The Ensembl automatic annotation of the genome sequence (https ://www.ensem bl.org/) was used to identify the following sequence: MyD88 (ENSC-ING00000017616). MyD88 cDNA was obtained using the GeneRacer kit (INVITROGEN, USA), as reported in Vizzini et al. 45 . In short, 5′-and 3′-RACE was directed by the listed primers reported in supplementary Table 1; the overlapping RACE was cloned into a pCRII vector (TA cloning kit, INVITROGEN, USA) and sequenced. cDNA was amplified by polymerase chain reaction (PCR), and the amplicon was cloned into a pCRII vector (TA cloning kit, INVITROGEN, USA) and sequenced. The full-length MyD88 cDNA was analysed with the ExPASy translation tool (https ://web.expas y.org/trans late/) to obtain the open reading frame (ORF) regions, the leader and trailer sequences (UTRs) and the nucleotide sequence, which was translated into a protein sequence. The protein sequence was located using the Basic Local Alignment Search Tool (BLAST, NCBI) to identify known protein sequences that are homologous to MyD88, and the physical/chemical parameters (e.g., molecular mass qRT-PCR. The differential expression of 27 LPS immune responsive genes was studied by qRT-PCR using the SYBR-Green method and the specific sets of primers listed in Table 3. qRT-PCR analysis was performed using an APPLIED BIOSYSTEMS 7500 Real-time PCR system 45 . Differential expression was determine in a 25 µl PCR mixture containing 2 μl of cDNA converted from 250 ng of total RNA, 300 nM primer (forward and reverse), and 12.5 μl of Power SYBR-Green PCR MasterMix (APPLIED BIOSYSTEMS). Amplification specificity was tested by real-time PCR melting analysis. To obtain sample quantification, the 2 −ΔΔCt method was used, and the relative changes in gene expression were analysed as described in the APPLIED BIOSYSTEMS Use Bulletin N.2 (P/N 4303859). The transcript levels from different tissues were normalised to that of actin to compensate for variations in the amount of RNA input. Relative expression was determined by dividing the normalised value of the target gene in each tissue by the normalised value obtained from the untreated tissue. To examine the time course of the response, LPS-treated ascidians (N = 4) were examined at incremental post-inoculation time points (1,2,4,8,12,24,48, and 72 h). Untreated ascidians (naïve) (N = 4) were used as controls.
Statistical methods. Statistical assessments of GO term enrichment and pathway analyses were performed by Fisher's exact test in combination with a robust false discovery rate (FDR) correction for multiple testing. The row p-value and FDR threshold were set as < 0.05.
Scientific RepoRtS | (2020) 10:11339 | https://doi.org/10.1038/s41598-020-68339-x www.nature.com/scientificreports/ Minitab 17 statistical software was used for the qRT-PCR data analysis. Statistical differences were estimated by one-way ANOVA, and the significance of differences among groups was determined by Tukey's t-test. The level of significance was set at a p-value ≤ 0.05. The data are presented as the means ± SD (N = 4).