Transcriptomic Analyses Reveal 2 and 4 Family Members of Cytochromes P450 (CYP) Involved in LPS Inflammatory Response in Pharynx of Ciona robusta

Cytochromes P450 (CYP) are enzymes responsible for the biotransformation of most endogenous and exogenous agents. The expression of each CYP is influenced by a unique combination of mechanisms and factors including genetic polymorphisms, induction by xenobiotics, and regulation by cytokines and hormones. In recent years, Ciona robusta, one of the closest living relatives of vertebrates, has become a model in various fields of biology, in particular for studying inflammatory response. Using an in vivo LPS exposure strategy, next-generation sequencing (NGS) and qRT-PCR combined with bioinformatics and in silico analyses, compared whole pharynx transcripts from naïve and LPS-exposed C. robusta, and we provide the first view of cytochrome genes expression and miRNA regulation in the inflammatory response induced by LPS in a hematopoietic organ. In C. robusta, cytochromes belonging to 2B,2C, 2J, 2U, 4B and 4F subfamilies were deregulated and miRNA network interactions suggest that different conserved and species-specific miRNAs are involved in post-transcriptional regulation of cytochrome genes and that there could be an interplay between specific miRNAs regulating both inflammation and cytochrome molecules in the inflammatory response in C. robusta.


Introduction
Cytochrome P450 enzymes (CYPs) consist of a large group of haemoproteins that catalyse a wide range of reactions involved in the oxidative metabolism of xenobiotics, such as drugs, pesticides, and environmental chemicals, or in the metabolism of endogenous compounds, such as steroid hormones, fatty acids, eicosanoids, and vitamins [1][2][3][4]. CYP genes are widely distributed in all domains of life, from bacteria, archaea, and viruses to higher plants and animals and may be classified in class, group, clan, family and subfamily [5]. The CYPs constitute one of the most diverse eukaryotic genes, with a wide complexity within and between species. CYP enzymes use molecular oxygen to modify substrate structure, critical in a huge number of physiological, ecological and toxicological processes. In chordates CYP genes, created by the 2R whole genome duplications, 11 distinct clans have been identified, each containing one or more CYP gene families. CYP constitutes a large gene class descended from a single common ancestor, termed the "cytochrome P450 genesis locus". Duplicated in tandem of progenitor CYP gene gave rise to a set of genes that were precursors of the 11 animal CYP clans: CYP Clans 2, 3, 4, 7, 19, 20, 26, 46, Figure 1. MA plot (a). X axis reports the Log2 mean expression of genes; Y axis reports the log2FC (threshold = 1.5). Red points are upregulated genes, blue points are down-regulated genes. PCA plot (b). Principal Component 1 (PC1) and PC2 of normalized data are reported in Y and X axis. Blue points represent C. robusta exposed to LPS (n = 3) and red points represent untreated ascidians (n = 3). The number of total points in PCA is the number of replicates used for NGS experiment (n = 6).
For functional annotations, unigenes (1227) were aligned to Gene Ontology (GO) terms based on the Protein ANalysis THrough Evolutionary Relationships (PANTHER) classification system for gene ontology annotations (pantherdb.org, 16.0 release); that it is connected to the Clusters of Orthologous Genes (COG) database (geneontology.org, release February 2021). The C. robusta transcripts were classified in the three GO subcategories:  Between genes responsive to LPS, different positively (+) enriched classes belonging to the following GO terms were evidenced: (i) Biological Processes: response to xenobiotic stimulus, response to chemicals, response to stimulus, response to drugs, response to immune system processes and organic acid metabolic processes; (ii) Cellular Components: extracellular space, extracellular region and extracellular region part; (iii) Molecular Function: monosaccharide binding, heme binding, tetrapyrrole binding, monooxygenase activity, amino acid transmembrane transport, oxidoreductase activity.
These results strongly indicated the presence of DE genes linked to Cytochrome P450 (Cyp 450) families and inflammatory response.
Indeed cytochrome P450 enzymes catalyze a wide range of reactions involved in the oxidative metabolism of xenobiotic, or in the metabolism of endogenous compounds; moreover, they use molecular oxygen to modify substrate structure, critical in a huge number of physiological, and toxicological processes.
Finally, the expression of xenobiotic-induced cytochromes is mediated by LPS stimulation, which also increases different cytokines which are well-known pro-inflammatory mediators, thus correlating a potential link between these molecules and immunity response processes.
A deeper investigation on DE genes evidenced different families of cytochromes that resulted in deregulation in NGS experiments. In particular, cytochrome P450 components of 2 and 4 families were deregulated (Table 1).

Alignment and Phylogenetic Analysis of C. robusta Cytochrome P450 LPS Deregulated Genes
The Cytochrome P450 LPS deregulated genes (Table 1) were aligned with vertebrate and invertebrate Cytochrome P450 2 and 4 family members ( Figure 3). Multiple alignments showed that all the 2 and 4 family members of C. robusta Cytochrome P450 possess the characteristic conserved motif of Cytochrome P450 ( Figure 3; Table 2). The cysteine residue coordinating heme iron result conserved in all the sequences, in a consensus sequence well-conserved in cytochrome P450 (FXXGX(H/R)XCXG), together with the EX1X2R motif located in helix K involved in salt bridge interactions, that are important for its tertiary structure and the correct incorporation of the heme cofactor, (AG)Gx(DE)T(TS) located in I helix and FDDER located in M Coil ( Figure 3; Table 2).
To understand the evolutionary relationship between the known 2 and 4 family members of Cytochrome P450 in vertebrates, invertebrates and C. robusta, we performed a phylogenetic analysis looking for orthologous sequences retrieved in the NCBI database (https://www.ncbi.nlm.nih.gov) (Accessed on 15 June 2021). Orthologous molecules were aligned using the ClustalW algorithm and neighbour-joining trees were constructed. A phylogenetic tree was constructed using the amino-acid sequence of D. melanogaster as an out-group (MEGA X program) ( Figure 4). The Cytochrome P450 in vertebrates and C. robusta formed two distinct major clusters. The first cluster includes the group of vertebrates, invertebrates and C. robusta Cytochrome P450 belong to 2 family members, the second group includes the group of vertebrates, invertebrates and C. robusta cytochrome P450 belong to 4 family members ( Figure 4). The phylogenetic analysis supported the idea of a conserved evolution of Cytochrome P450 from a common ancestral gene among protochordates and vertebrates.

Differential Gene Expression of Cytochrome P450 Genes in Different Tissues
The spatial expression pattern of Cytochrome mRNA (Table 1) in adult C. robusta was investigated by qRT-PCR analysis of total RNA samples from different tissue involved in inflammatory response (pharynx, ovary, stomach and intestine). As shown in Figure 5, Cytochrome P450 expression was higher in the stomach and the intestine, and lower in the ovary and pharynx tissue; 2B10, 2C15, 2J6, 4B1 and 4F6 family members expression was higher in the intestine and the stomach, and lower in the ovary and pharynx tissue; 2U1 family member expression was higher in the ovary, and lower in the intestine, the stomach and pharynx tissue; 2U1-like family member expression was higher in the pharynx and intestine and lower in the stomach and ovary tissue. Values, plotted as mean ± SD, were inferred from four ascidians. * p < 0.05, ** p < 0.005, *** p < 0.001.

Analyses of the Expression of Cytochrome 450 and Cytokines Genes under LPS Exposure
Analyses of the time-course expression of Cytochrome and Cytokines in the pharynx inflammatory response induced by LPS in C. robusta were performed at time points from 0 to 48 h post-LPS challenge by qRT-PCR ( Figure 6). The heatmap shows that cytochrome transcripts were significantly modulated in response to LPS during the 48-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 pro-inflammatory cytokines Mif1, Mif2, Il-17-1, Il-17-2, and Tnf-α and Cyp450 2C15, 2J6, 2C42 and the second comprises Nf-κB and Tgf-β, Il-17-3 and Cyp450 2U1, 2U1-like, 2B10-like, 4B1, 4F6. Specifically, the heatmap highlighted that the inflammatory cytokines Mif1, Il-17-1, Il-17-3, and Tnf-α were upregulated between 1 and 4 h of LPS exposure (p-value < 0.05) and Tnf-α reached its maximum expression level after 2 h of LPS exposure. Notably, Nf-κB and transforming growth factor β (Tgf-β) transcripts displayed a significant increase after 4 h of LPS exposure (p-value < 0.05). On the other hand, after 8 h of exposure, Il-17-1, Il-17-2 levels began to increase at 1 h. Nf-κB and Tgf-β show a second significant increase after 48 h of exposure (p-value < 0.05). Cyp450 2C15, 2J6, 2C42 mRNA were upregulated between 1 and 2 h of LPS exposure, while Cyp450 2U1,2U1-like, 2B10-like, 4B1, 4F6 were upregulated between 1 and 4 h of LPS exposure. These findings suggest an involvement of cytokines in modulating the expression of Cytochromes in response to LPS exposure.  (1-48 h). Time course of gene expression in the pharynx of C. robusta exposed to LPS compared with the gene expression in untreated ascidians. To compute the heatmap was chosen to use the Complete linkage as clustering algorithm, and the Pearson correlation as distance measurement method. Values are represented according to the z-score, which is measured in terms of standard deviations from the mean.

miRNA-Target Interaction Prediction
miRNA-target interaction of deregulated genes belonging to cytochrome (Table 1) and inflammation was investigated, by using the miRNA Target Interactor Predictor (miRNATIP) algorithm to explore how deregulation of Cytochrome P450 enzymes may be driven by non-coding RNA intervention during the inflammation process. miRNAtarget prediction evidenced several miRNAs interacting with both DE genes linked to the Cytochrome P450 family and inflammation. A total of 10,002 interactions were predicted for DE Cytochromes and 3905 interactions were predicted for DE cytokines. Results were then filtered for the lowest free energy values (free energy < −12 Kcal/mol). Table 3 shows a detailed list of miRNA-target interactions, showing just interactions with more than two mRNA targets for each miRNA are showed. A total of 42 of all miRNAs identified by the miRNATIP algorithm were common to the two differentially expressed genes classes, and they belong to 36 different miRNA families. Seven of 36 families were conserved through species and 29 were species-specific. Table 3. Table shows miRNA-mRNA target interactions for deregulated genes linked to cytochrome (blue) and inflammation (red). All interactions showed in the table are divided for conserved and species-specific miRNAs. They are filtered for energy values and are ordered by number of interacting gene for each miRNA. Just results of miRNA-target interactions with more than two interacting targets are showed.

Conserved miRNAs
Target Gene ID
See Table 4 for a detailed list of predicted miRNAs, conserved through species and species-specific, that are common to the Cytochrome P450 and inflammation genes. See supplementary material for all miRNA-target interactions (Supplementary Materials). Table 4. Table shows all miRNAs that are common to deregulated cytochrome genes and cytokines, filtered for low energy values.

Conserved miRNAS Species-Specific miRNAs
cin-miR-5596b-3p cin-miR-5603-3p cin-miR-5607-3p cin-miR-5609-5p cin-miR-5612-5p cin-miR-4219-5p 2.6. Cytochrome and Inflammation miRNA-Target Network Reconstruction miRNA-target interaction network reconstruction allows to easily visualise relevant interactions between different molecules. To build the graph, different steps were required: first, DE genes of Cytochromes and inflammation networks were downloaded into the Search Tool for Recurring Instances of Neighbouring Genes (STRING) database (string-db.org, release 11.5, accessed on 12 August 2021), second, the two different networks produced (cytochromes and inflammation) were imported into Cytoscape app (httos://cytoscape.org, release March 2017) (Accessed on 13 August 2021), to generate graph.ml files, and third, the networks were analysed using R package and R studio (https://www.r-project.org/, release October 2020) (Accessed on 17 August 2021); (https://rstudio. com/products/rstudio/, release October 2020) (Accessed on 17 August 2021), and finally miRNA prediction results were integrated into two different networks ( Figure 7). Interaction network between miRNAs (green), DE Cytochrome genes (blue) and DE inflammatory cytokines (yellow). The network analysis was performed using R package and R studio, and miRNA prediction results were integrated into the two different networks of DE Cytochrome genes and inflammatory cytokines.
The miRNA-target network reconstruction put in evidence that some miRNAs interacting with more of two different targets, some regulating different Cytochrome genes, others regulating both genes of Cytochrome and cytokines (Table 3).
These in silico pieces of evidence let us hypothesize that there could be an interplay between specific miRNAs regulating both the inflammation process and cytochrome molecules.

Transcription Factor Orthologue Identification of Cyp450 Response Elements
After network reconstruction, Cyp450 gene transcriptional regulation by specific transcription factors (TFs) was investigated. Indeed, as known by scientific literature, TFs bind gene response elements and activate mRNA transcription of target genes.
Transcription factors identified in C. robusta genome are ci-Hnf-1 (ncbi gene ID: 778644), ci-Hnf-4 (ncbi gene ID: 778645), ci-AhR (ncbi gene ID: 778536), ci-Rxr (ncbi gene ID: 778746), ci-C/Ebp (ncbi gene ID: 778557) and ci-Vdr-a (ncbi gene ID: 778791). In silico analyses of TF binding sites prediction were finally carried out. Site Tracking and Recognition (SiTaR) tool (https://sbi.hki-jena.de/sitar/, V 0.1) (Accessed on 24 August 2021) was applied to all mentioned TFs and deregulated Cytochrome genes produced by NGS. Result analysis showed that just the ci-Vdr-a transcription factor potentially interacts with most deregulated cytochrome genes analysed (see Supplementary Material for detailed analysis results). Figure 8 shows a schematic representation of the potential transcriptional regulation of Cyp450 genes in response to TFs and non-coding RNAs.

Discussion
Cytochrome P450 (CYP) enzymes are an ancient superfamily of monooxygenase proteins found in all domains of life [35]. Genome sequencing of diverse species showed 12,000 known CYP genes [36]. Cytochromes are expressed in many different tissues of the human body. They are found mostly in intestinal and hepatic tissues. In this study by an in vivo LPS exposure strategy and next-generation sequencing (NGS), we found that in , and amphibians. The CYP2U and CYP2R genes were present in the vertebrate ancestor and are shared across all vertebrate classes; CYP2D and CYP2J are not found in actinopterygian (ray-finned fish) species [37]. Members of this family were found in arthropods and crustacean species, such as CYP2L in lobster [38], and a significant number of CYP2-like genes were found in the sea urchin [39]. In C. robusta LPS Analyses of the time-course expression of Cytochrome P 450 and Cytokines in the pharynx inflammatory response induced by LPS showed that Cyp450 2C15, 2J6, 2C42 were upregulated between 1 and 2 h. In contrast Cyp450 2U1,2U1-like, 2B10-like, 4F3,4F6 were upregulated between 1 and 4 h. While pro-inflammatory cytokine Tnf-α reached its maximum expression level after 2 h of LPS exposure and Nf-κB and transforming growth factor β (Tgf-β) transcripts displayed a significant increase after 4 h of LPS exposure. These findings suggest an involvement of cytokines and Nf-κB in downregulation of the Cytochromes response upon LPS exposure after 4h. In humans, rats, and mice [40][41][42] variations in hepatic CYP expression are common responses to infection and inflammation. Morgan [43] reported that levels of constitutively expressed CYP genes were downregulated in the rat liver when inflammation was induced with lipopolysaccharide (LPS). Furthermore, the changes of pro-inflammatory mediators TNF and IL-1 correlate with the changes in CYP expression and enzymatic activity during infection and inflammation [18][19][20]. Increasing these cytokines appears to strongly affect the pathways regulating inducible CYP expression and LPS markedly alters CYPs' expression levels during inflammation as Cyp2c29 and Cyp 2B10 mRNA in the liver [44]. Lipopolysaccharide (LPS) treatment of rats suppresses Cyp 4F4 and 4F5 expression by 50 and 40%, respectively. In this animal model, a change in the expression of Cyp 4F4 and 4F5 mRNA was observed at 24 h compared to the controls. At time points after 24 h the expression of both isoforms increases dramatically reaching the highest levels at 2 weeks post-injury. These results are consistent with the notion that immediately after injury, concentrations of leukotriene and prostaglandin mediators immediately after injury are elevated by decreased Cyp 4F concentrations. Increased Cyp 4F expression leads to diminished concentrations of leukotriene and prostaglandin mediators and then to recovery and repair [45] In vitro models NF-kB can control the expression of CYP1A1, CYP2B1, CYP2C11, CYP2D5, CYP2E1, and CYP3A7 via interaction with the promoters of these genes, leading to downregulation [46]. In vitro models show that distinct mechanisms are involved in the downregulation of CYP enzyme expression and Cytokine-mediated alteration of gene transcription. De jong [47] suggests that they can be the main regulatory mechanism accountable for changing Cyp450 activity upon inflammation [47]. The mechanisms of downregulation of Cytochromes P450 upon inflammation is an area of intense study. Increasing attention is being given to the potential post-transcriptional mechanisms that could regulate P450 enzymes in inflammation, as the involvement of MicroRNAs. MicroRNAs (miRNAs) inhibit gene expression largely in a posttranscriptional manner by interacting with miRNA response elements (MREs) in the mRNA 3 -UTRs of the target genes [48]. In C. robusta the NGS strategy identified transcripts that were classified as non-coding RNA (ncRNA). Of these 201, 90 were miRNAs. Both conserved (cin-let-7d-5p, cin-mir-153-5p, and cin-mir-92c-5p) and species-specific miRNAs were found involved in post-transcriptional regulation of four pathways (Tgfβ,Wnt, FoxO, and Hh) linked to pharynx tissue homeostasis [33]. Bioinformatics and miRNA-target prediction analysis show that conserved and species-specific miRNAs can be involved in regulating post-transcriptional mechanisms of both gene categories: Cytochrome P450 and cytokines. Particularly cin-let-7 and cin-miR-92 seem vastly involved in post-transcriptional regulation during LPS induced inflammatory processes. In vertebrates, regulation of hepatic CYP levels is involved a wide variety of ligand-activated transcription factors and nuclear receptors. In inflammation, in the regulation of CYP mRNA levels, have involved Transcription factors including the nuclear receptors pregnane X receptor (PXR), the constitutive androstane receptor (CAR), their dimerization partner retinoid X receptor (RXR), the aryl hydrocarbon receptor (AhR) and human nuclear factors (HNFs) that are also responsible for the downregulation of Cytochromes P450 in inflammation. Nuclear Receptors as CAR increase transcription of the human CYP3A4/5, CYP2C9, CYP2C19, and CYP1A2 genes upon drug treatment [49][50][51]. Nuclear receptor PXR and CAR are involved in the mechanism of changes in gene transcription of major Cytochromes in inflammation. Vast evidence shows that inflammation represses PXR levels, leading to downregulation of important CYP enzymes [52,53] and LPS treatment in mice led to functional repression of PXR's dimerization partner RXR [54]. A likely scenario is that the suppressive action of inflammation on PXR expression is mediated through NF-kB activation since Zhou X et al.
showed that NF-kB directly interacts with a functional binding site in the PXR promotor to suppress its transcriptional expression [55]. C. robusta genome reveals only a single putative orthologue to vertebrate nuclear receptors PXR which was most similar to mammalian PXRs [56] with the ability to adapt to different ligands [57]. Furthermore, in silico analysis showed that the nuclear transcription factor ci-Vdr-a belongs to a nuclear hormone receptor family that is homologue to human PXR (ncbi ID: 8856), can potentially interact with many of deregulated Cytochrome genes identified by NGS. Thus, we suggest that ci-Vdr-a could potentially act as its orthologue protein PXR, intervening in the transcriptional regulation of Cytochrome genes during the inflammation process.
Concluding, we assume that in C. robusta ( Figure 8) the modulation of cytochromes belonging to 2B, 2C, 2J, 2U, 4B and 4F subfamilies genes expression during inflammation may proceed through transcriptional downregulation of transcription factors, interference with dimerization/translocation of (nuclear) transcription factors, involvement of NF-kB, and inhibition gene expression in a posttranscriptional manner by interacting with miRNA. The control of Cytochrome P450 gene expression happens in a coordinated manner with the control of gene expression of immune genes as Cytokines and NF-kB according to highly conserved mechanisms during evolution.

Tunicates and LPS Injection
The animal model Ciona robusta was formerly classified as Ciona intestinalis. Molecular studies have confirmed that C. intestinalis constitutes a compilation of species rather than a single speciements [57][58][59][60][61]. C. robusta were collected from Sciacca harbour (Sicily, Italy) and were acclimatized and maintained as reported in Arizza et al. [34]. An LPS solution (Escherichia coli 055:B5, LPS, SIGMA-ALDRICH, Saint-Louis, MI, USA) was prepared in a sterile salt medium (12 mM CaCl2, 11 mM KCl, 26 mM MgCl2, 43 mM Tris HCl, 0.4 M NaCl, pH 8.0). One hundred microliters 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 48 h) and pharynx, ovary, intestine and stomach tissues of naïve were immediately soaked in RNAlater tissue collection solution (AMBION, Austin, TX, USA) and stored at −80 • C. Total RNA extraction was performed using an RNAqueous-Midi kit purification system (AMBION, Austin, TX, USA) as reported in Arizza et al. [34].

RNA Sequencing (RNA-Seq)
The RNA purity and quality of total RNA extracted from the pharynx of C. robusta that were naïve 3 replicates (n = 3) and that were exposed to LPS for 4 h 3 replicates (n = 3) were examined by NanoDrop and Agilent RNA 6000 Nano kits on an Agilent 2100 Bioanalyser (AGILENT, Santa Clara, CA, 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 [34]. All transcripts produced by NGS were annotated by Ensembl database (https://www.ensembl.org/index.html, release August 2020) (Accessed on 3 March 2021). Differential expression between treated (4hrs LPS induction) and untreated genes was performed by BMR genomics [34]. It was performed using edgeR software. It allows to estimate the negative binomial variance parameter globally across all genes. All the data were then normalized by setting the false discovery rate (FDR) to ≤0.05 and the absolute value of the log2 fold change (logFC) to 1.5. A MA plot was performed to visualize the up-regulated and down-regulated genes from RNA-sequencing. The plot visualizes the differences between measurements taken in two samples, by transforming the data onto M (log ratio) and A (mean average) scales, then plotting these values. Results were calculated using ggpubr pakage of CRAN library (https://cran.r-project.org/web/packages/ggpubr/index.html) (Accessed on 30 September 2021).
To visualize the data, the log2 mean expression is calculated, which is the mean of values of treated samples, plus the log2 of the mean of values of controls, all divided by 2. We also provided a PCA plot to show the distribution of LPS-treated and untreated samples. The PCA clearly shows two different cluster samples, one grouping LPS treated samples and one grouping controls. We made a PCA plot of normalized data, showing the first two principal components (PCs) of treated and untreated samples. Results were calculated using PCA tool R library (https://bioconductor.org/packages/release/bioc/vignettes/ PCAtools/inst/doc/PCAtools.html#conduct-principal-component-analysis-pca) (Accessed on 30 September 2021).

Alignments and Phylogenetic Analyses
Multiple alignments of sequences were carried out using CLC (Version 7.0.0). Phylogenetic trees were designed in MEGA X maintaining the bootstrap value of 1000 bootstrap iterations (neighbour-joining method). The accession numbers are listed in Table 5.

qRT-PCR
The differential expression of Cytochrome P450 LPS responsive genes was studied by qRT-PCR using the SYBR-Green method and the specific sets of primers listed in Table 6. qRT-PCR analysis was performed using an Applied Biosystems 7500 Real-time PCR system [34]. Differential expression was determined 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, Waltham, MA, USA). 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 analyzed as described in the Applied Biosystems Use Bulletin N.2 (P/N 4303859). The transcript levels from different tissues were normalized to that of actin to compensate for variations in the amount of RNA input. Relative expression was determined by dividing the normalized value of the target gene in each tissue by the normalized value obtained from the untreated tissue. To examine the time course of the response, LPS-treated ascidians 4 replicates (n = 4) were examined at incremental post-inoculation time points (1, 2, 4, 8, 12, 24, and 48 h). Untreated ascidians (naïve) 4 replicates (n = 4) were used as controls. Table 6. Primers used for qRT-PCR.

Functional Enrichment Analysis
Gene enrichment of C. robusta transcripts produced by NGS was made through The Clusters of Orthologous Genes (COG) database (geneontology.org, release February 2021). The three Gene Ontology (GO) subcategories were investigated, and they are the following: (i) Biological Process (BP); (ii) Molecular Functions (MF); (iii) Cellular Components (CC); (iii) p-value and adjusted p-value thresholds were set to <0.05 to have statistically significant results.
The Protein Analysis Through Evolutionary Relationships (PANTHER GO-slim analysis tool) (pantherdb.org, 16.0 release) System connected to the COG database was used to expand the annotation data of the gene and protein families obtained from GO. The PANTHER "GO-slim" analysis mode was selected to have more reliable and accurate results compared to the GO "GO-complete" annotation mode.

miRNA-mRNA Target Interaction Prediction
RNA-RNA interaction predictions (miRNA-mRNA target) were performed through the miRNA target interaction predictor (miRNATIP) algorithm [63]. The algorithm filters out wrongly predicted interactions, by means of the computation of the binding free energy (free energy threshold <−12Kcal/mol). This last step allows for the strengthening of the power of predicted interactions calculated by the algorithm. Indeed, as it is known, the lower the free energy of two paired RNA strands, the more energy is needed to disrupt this duplex formation, thus deducing that the stronger RNA-RNA binding is when the binding energy is low, and there is a more stable state thermodynamically. To test the best free energy thresholds in the miRNATIP algorithm we did some experiments in C. elegans and in H. sapiens to verify the most favorable ∆G values of binding sites, and they were <−6 and <−7kcal/mol. As some other scientific papers related to RNA-RNA interaction prediction tools used lower energy thresholds to evaluate RNA-RNA stability, we decided then to lower these energy thresholds to <−12kcal/mol energy values also in C. robusta samples [64].

Network Construction and Visualisation
To build interaction networks between DE genes linked to Cytochrome and inflammation and interacting miRNAs, different steps were achieved: first, both DE cytochrome and inflammation genes were downloaded in the STRING database (db) (string-db.org, release 11.5, accessed on 12 August 2021). STRING db, is an open-source platform for reconstructing interaction networks by using direct, indirect and functional interactions between molecules. In the second step, the two different networks previously produced (Cytochrome and inflammation) were imported into the Cytoscape app (https://cytoscape.org, release March 2017) (Accessed on 13 August 2021), to generate graph.ml files.
Cytoscape is an open-source software platform for visualizing molecular interaction networks and biological pathways and integrating these networks with annotations.
In particular, "igraph" package of R was used to analyses networks produced by STRING and to generate the image of integrated networks.
Both orthoDB and REGULATOR collect and documents the transcription factors of metazoan species, between which C. robusta, giving information on the divergence in species evolution. The tools are also linked to other web resources such as InerPro, IntAct, Pfam, STRING, GO annotation, NCBI.
Transcription Factor binding site prediction was also made to predict if orthologue TFs can potentially interact with genomic sequences of deregulated Cytochrome genes. The analysis was performed using the Site Tracking and Recognition (SiTaR) tool (https: //sbi.hki-jena.de/sitar/, V 0.1) (Accessed on 24 August 2021). The method is based on the idea that it can be calculated several motifs with a given nucleotide content and a given number of mismatches in a random sequence. These predicted numbers can then be compared with the real occurrences of motifs in a query sequence. Moreover, this method uses a novel approach compared with position weight matrices (PWM) and Hidden Markow Models, and it is based on the research of motifs on input sequence, where for each motif found, scores are assigned depending on the non-randomness of the motif's occurrence, the number of matching searching motifs and the number of mismatches. The method also has the advantage to significantly reduce the number of false-positive predictions significantly.

Statistical Methods
Fisher's exact test performed statistical assessments of GO term enrichment and pathway analyses 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.
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).