Core ISGs

We identified a core set of 62 genes (“core vertebrate ISGs”, hereinafter core^vert ISGs) that were up-regulated by IFN by all 10 species analysed in this study, with an additional 28 genes up-regulated specifically in the nine mammalian species (“core mammalian ISGs”, hereinafter core^mamm ISGs) (Table 1). The core^vert ISGs represent the ancestral functions of the IFN system and include genes encoding proteins broadly involved in (i) orchestrating antigen presentation, (ii) IFN induction and response, (iii) IFN suppression, (iv) ubiquitination and protein degradation, (v) cell signalling and apoptosis, and (vi) antiviral responses (Table 1, Fig 1E).

Nine of the 62 core^vert ISGs (e.g., various HLA genes, TAP1, ERAP1, etc.) are involved in the generation, trimming, loading, and presentation of MHC-I–restricted antigens, thus providing a direct link between the IFN response and the adaptive immune response via the CD8 T-cell response (Table 1). Additionally, components of the immunoproteasome (e.g., PSMB8 and PSMB9) were specifically up-regulated in the mammalian core (Table 1), reflecting previous studies reporting the absence of the immunoproteasome in birds [10].

We found that various core^vert ISGs are involved in IFN induction and response, including pattern recognition receptors, the key adaptor molecule MyD88, and transcription factors. The classical sensors for RNA PAMPs, including IFIH1/MDA5, DHX58/LGP2, and TLR3 (and RIG-I/DDX58 in the mammals) were among the core^vert ISGs (Table 1, Fig 2). RIG-I is known to be absent in chickens (and other galliformes) but is present and active in other birds, including ducks and geese [11–13]. Furthermore, TRIM25, a ubiquitin E3 ligase responsible for ubiquitinating RIG-I, was also a core^vert ISG [14]. DAI/ZBP1, which was originally classified as a DNA sensor but is now thought to be an RNA sensor [15], was also highly up-regulated by IFN in every mammalian species except humans (Fig 2). No DNA sensors were found among the core^vert ISGs. However, cGAS was found to be an ISG in every mammalian species in this study (Fig 2). By contrast, the DNA sensor AIM2 was only up-regulated in human cells. Interestingly, whilst the basal expression (defined in terms of fragments per kilobase mapped values [FPKM] in the absence of IFN treatment) of RNA sensors was very low, many of the genes in the literature associated with DNA sensing are constitutively transcribed (Fig 2).

Open in a separate window

Fig 2

Basal transcription levels and IFN-induced expression of genes related to PAMP sensing and IFN induction and response.

Boxplots showing differential expression (log2FC) in response to IFN and basal transcription levels (expressed as FPKM) of genes associated with pattern recognition (sensors), downstream signal transduction (adapters), and transcription factors related to either IFN induction or response (transcription factors). Every ortholog for each gene is indicated with a dot coloured according to their presence in the DNA-, RNA-, or both DNA- and RNA-sensing pathways (S1 Data). FPKM, fragments per kilobase mapped values; IFN, interferon; log2FC, log2 fold change; PAMP, pathogen-associated molecular pattern.

With the exception of MyD88, a key adaptor involved in both the RNA- and DNA-sensing pathways [16], we observed limited up-regulation among genes involved downstream of nucleic acid detection (Fig 2). On the other hand, core ISGs included key transcription factors involved in IFN induction and response (IRF1, IRF7, STAT1 and STAT2 are all core^vert ISGs in addition to IRF9 among the core^mamm ISGs) (Fig 2). Importantly, several ISGs that play a role in the suppression of the IFN system, including USP18, USP25, IFI35, and SOCS1 were up-regulated in all species under examination. The encoded proteins of these genes target different points in the IFN response. Thus, negative regulation of the IFN response is multifaceted and a fundamental, ancestral failsafe necessary to avoid excessive/perpetual up-regulation of IFN-induced pathways. Among the core^vert ISGs, we found several genes relating to ubiquitination, such as the ring finger proteins RNF213 and RNF19B (Table 1, Fig 1E), highlighting protein modification as part of the IFN response. Interestingly, N4BP1, originally identified as a target of Nedd4-mediated ubiquitination, has not previously been directly linked to the IFN response.

The core^vert ISGs contained 14 IFN-induced antiviral factors such as MX1, IFIT2, and viperin (Table 1). Interestingly, when we assembled a list of 40 genes that either create a cellular environment hostile to or act directly upon the virus lifecycle (based upon the scientific literature; S2 Table), we noticed that 75% of these antiviral genes were ISGs in at least eight of the 10 species analysed in this study (Fig 3A). In addition, antiviral ISGs were up-regulated to a significantly higher extent than randomly sampled ISGs (P < 0.01 for each species, Fig 3B). Furthermore, some well-studied antiviral ISGs were not up-regulated by IFN in certain species. For example, OAS2 was not up-regulated in the rat, SAMHD1 was not up-regulated in the horse, OASL was not up-regulated in either the cow or the sheep, and the IFITM genes were not up-regulated in the dog (Fig 3C). In general, genes encoding antiviral factors were transcribed to minimal levels in the absence of IFN. Indeed, the median FPKM level for antiviral genes was lower than that of the overall interferome for every antiviral factor except SAT1, SHISA5, and the IFITM genes. Of interest, we noticed particularly high FPKM values for IFITM1 and 3 in the rat and IFITM2 in the microbat (Fig 3C).

Open in a separate window

Fig 3

Properties of antiviral ISGs.

(A) Sinaplot showing the differential expression values (log2FC) of 40 genes previously published as exerting antiviral activity (red dots) (S2 Table) as opposed to the rest of the ISGs (grey dots). ISGs (including antiviral ISGs) were allocated to bins according to the number of species in which they were found to be up- or down-regulated (i.e., the core^vert ISGs are bin 10). The majority of antiviral ISGs (red) were found to be up-regulated in at least eight species (S1 Data). (B) Graph showing the extent of up-regulation of antiviral ISGs compared to nonantiviral ISGs. The mean log2FC of 40 known antiviral ISGs (indicated with an asterisk on the plot) is compared to 100 samplings of 40 randomly selected ISGs from the interferome of each species (box and whiskers). In all cases, antiviral ISGs are up-regulated to a significantly greater extent as compared to nonantiviral ISGs (P < 0.01 for each species). The code used for random sampling and the generation of Fig 3B is available in S2 Data, with the required input files available as S5 Data and S6 Data. (C) Boxplots showing basal transcription levels (expressed as FPKM) and differential expression (log2FC) in response to IFN of known antiviral ISGs as in panels A and B. Every ortholog for each gene is indicated with a dot coloured according to species. The median FPKM value for the entire interferome is indicated with a broken line (S1 Data). FPKM, fragments per kilobase mapped values; IFN, interferon; ISGs, interferon-stimulated genes; log2FC, log2 fold change.

Two core^vert genes, CD47 and IL15RA, encode proteins involved in signalling to components of the adaptive immune system. CD47 is involved in a variety of biological roles, including leukocyte and dendritic cell migration, the development of antigen presenting cells, and immune apoptosis, and it also provides a ‘don’t eat me’ signal [17]. The IL15–IL15Rα axis is well characterised as being important in the promotion of both natural killer cells and a variety of T cell populations, including activated CD8 T cells [18].

Genes previously unrelated to the IFN response

Among the core^vert ISGs, we identified a number of genes with few or no reported associations with the type I IFN response (Table 1). Several of these genes have been studied extensively, but not always in the context of the IFN response. DNAJC13, for example, is reported to be involved in endosome trafficking, and it has been closely linked to Parkinson’s disease [19]. Zinc Finger CCHC-Type Containing 2 (ZCCHC2) has nucleic acid–binding properties and, interestingly, contains a single nucleotide polymorphism (SNP) associated with insect bite hypersensitivity in Exmoor ponies [20]. The fragile X mental retardation (FMR1) gene encodes an RNA-binding protein that plays a role in intracellular RNA transport and in the regulation of translation of target mRNAs. FMR1 was not previously linked to the IFN response, although it has recently been shown to be a proviral factor for influenza virus and previously was shown to induce mild restriction of HIV-1 [21,22]. Cap methyltransferase 1 (CMTR1), also known as ISG95, binds to RNA pol II and is a 2′-O-ribose methyltransferase that participates in the conversion of cap0 to cap1 type transcripts [23]. Interestingly CMTR1 is also known as an important component of IFIT-mediated antiviral activity [24]. We hypothesise that these genes, as core^vert ISGs, likely play fundamental roles in host immunity that are underappreciated or have yet to be fully determined.

Our data also indicate a relatively underappreciated link between local synthesis of early components of the complement system and the type I IFN response [25,26]. C2 was among our core^mamm ISGs. In addition, we found that C1r and C1s, essential components of the C1 complex, were up-regulated by IFN in cells from all species analysed in our study with the exception of the cow. Interestingly, this was also the case for a negative regulator of C1r and C1s (SERPING1/C1-INH).

Interferon-repressed genes and epigenetic controls of the IFN response

Our data enabled an unprecedented opportunity to investigate the interferon-repressed genes (IRGs), which, to date, have received comparatively little attention with regards to their involvement in the innate immune system. Unlike ISGs, the extent of down-regulation of IRGs across all the 10 species used in this study was relatively modest (overall average −0.56 log2FC for IRGs as compared to 1.64 log2FC for ISGs). We found no IRGs shared by all species, although this may reflect the low fold change in expression and/or greater variability in the response of individual genes. This result could also imply that there has been less conservation of the down-regulated genes over time. The most consistently down-regulated genes were FAM117B and KDM5B, which were both significantly down-regulated in eight of the 10 species analysed in this study. Relatively little is known about the function of FAM117B with the exception that it is a risk factor for sarcoidosis [27]. On the other hand, it has been shown that suppression of the KDM5B gene product, a H3K4 demethylase causing transcriptional repression, results in increased expression of IFN-β and other inflammatory cytokines following infection with respiratory syncytial virus (RSV) [28]. We also noticed that, with the exception of the rat, every species analysed down-regulated at least one KDM gene in response to IFN.

It is already established that another form of epigenetic control, acetylation, is also required for robust ISG transcription [29]. ANP32A, a protein involved in acetyltransferase inhibition [30], was down-regulated in five species. Interestingly, ANP32A has been shown to be a host component necessary for influenza virus replication and influences the ability of the virus to replicate in a given animal species [31].

Our data have thus revealed that key epigenetic factors regulating ISG transcription are themselves frequently responsive to IFN treatment. Along these lines, we also investigated the presence of noncoding RNAs (ncRNAs), a class of RNA molecule increasingly recognised as being important in the antiviral response [32,33]. We found that in human cells 75 long intergenic noncoding RNAs (lincRNAs) were differentially expressed (38 up-regulated, 37 down-regulated) in response to IFN (S2 Fig), including NEAT1, a lincRNA that has previously been associated with viral infections [34–36]. As other genomes become increasingly well annotated, it will be possible to resolve a more in-depth understanding of the impact that ncRNAs play in the control of the innate immune system.

Primary cell isolation and cell culture

Ex vivo skin samples were collected from chickens (n = 3), cows (n = 4), sheep (n = 3), horses (n = 3), a dog (n = 1), and pigs (n = 4) and primary fibroblasts isolated following an explant procedure. Briefly, the hair or feathers were removed from the skin prior to disinfection by soaking in 70% ethanol. After rinsing in PBS, the skin was removed from the underlying tissues, cut into circa 3-mm square explants, and added to cell culture dishes (without media, squamous surface uppermost) for one hour at 37°C before adding Dulbecco’s modified Eagle medium (DMEM) (Gibco) supplemented with 10% fetal bovine serum (FBS) (Gibco), 1% penicillin/streptomycin (p/s) (Sigma) and 100 U/ml nystatin (Sigma). Human primary dermal fibroblasts were purchased from PromoCell (catalogue number C-12302). Rat primary dermal fibroblasts were purchased from the European Collection of Authenticated Cell Cultures (ECACC) (catalogue number 06090769). M. lucifugus (little brown bat, representative of the microbats) primary dermal fibroblast cells, isolated from individuals caught in Oregon, United States of America, were kindly provided by William Kohler [46]. P. vampyrus (the large flying fox fruit bat, representative of the megabats) cells (PVK4) are an immortalised kidney cell line kindly provided by Megan Shaw [45]. The origin of each cell used in this study is summarized in S5 Table.

Human, rat, and dog cells were cultured in fibroblast growth medium 2 (PromoCell) supplemented with 10% FBS and p/s. 293T and PVK4 were all cultured in DMEM supplemented with 10% FBS and p/s. All cell cultures were incubated at 37°C with 5% CO2 in a humidified atmosphere.

Virus infectivity assays

IFN- or mock-treated cells were challenged with infectious VLPs of envelope minus vesicular stomatitis virus (VSV-ΔG-GFP) decorated with a VSV-G envelope (provided in trans during VLP production) in order to confirm the antiviral state of the cells at the time of RNA harvest essentially as already described [47,48]. Briefly, cells were harvested by trypsinization and fixed in 5% formaldehyde. The number of infected cells in IFN- and mock-treated cells was assessed by flow cytometry.

IFN treatment of cells

Parallel sets of cells were plated in multiwell plates 24 or 48 hours prior to IFN treatment and incubated at 37°C. Cells were treated with either 1,000 U/ml Universal interferon (UIFN, PBL InterferonSource), 200 ng/ml canine IFNα (Kingfisher), 1,000 U/ml porcine IFNα1 (Stratech), or 200 ng/ml chicken IFNα (AbD Serotec). Mock treatment was performed in parallel using DMEM lacking IFN. Cells were incubated for the indicated time period at 37°C, washed with PBS, and either lysed in Trizol (Thermo Fisher) for RNA extraction or challenged with VSV-ΔG-GFP to assess the antiviral state. Cells were only further processed for RNAseq analyses when they were in an antiviral state. In this study, cells were considered in an antiviral state when IFN stimulation induced at least 75% inhibition of VSV-ΔG-GFP infectivity (value chosen as average of three independent experiments) (S1 Fig). Pilot experiments were performed for each cell type in order to optimise conditions necessary for cells to reach an antiviral state. With the exception of dog cells, all cells reached an antiviral state after four hours of IFN treatment. The antiviral state in the primary canine cells isolated in these experiments required 24 hours of treatment with canine IFNα (S1 Fig). For comparative purposes, a complete set of experiments ranging from IFN stimulation to interferome analysis was performed in parallel in pig cells stimulated with either UIFN or porcine IFN-α (S3 Fig).

RNA extraction and sequencing

RNA was extracted using Trizol (Thermo Fisher) and RNeasy (Qiagen) protocols. Briefly, chloroform was added to the RNA-containing phase of the Trizol sample and centrifuged. The aqueous phase was then mixed with ethanol and purified using RNeasy columns, incorporating an on-column DNase step (Qiagen) to ensure the complete removal of genomic DNA. Total RNA samples were quantified using the Qubit (Thermo Fisher) and were assessed for integrity using the Bioanalyser pico eukaryotic II RNA chip (Agilent). Only samples with an RNA integrity number (RIN) value > 9 were taken forward for library preparation.

Libraries of mock- and IFN-treated cell RNA were assembled using equal masses of total RNA. The External RNA Controls Consortium (ERCC) spike control was added to the total RNA sample in order to assess library quality following sequencing. RNA samples were first enriched for mRNA by selecting poly(A) RNA using the Dynabeads mRNA DIRECT Micro purification kit (Thermo Fisher). The eluted RNA was used to generate RNAseq libraries using the Ion Total RNA-seq Kit v2 (Thermo Fisher) following the manufacturer’s instructions, with the exception that RNA samples were sheared for just 1.5 minutes. Amplified and purified libraries were checked for quality and quantity using the Agilent Tapestation (D1000 tape) and Qubit (hsDNA assay). Libraries were run on the Ion Proton (Thermo Fisher) according to the manufacturer’s instructions.

Data QC

Raw data were trimmed and assessed for quality using FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We performed multidimensional scaling (MDS) of normalised counts per million data using EdgeR (Bioconductor) in order to assess the impact of IFN and assessed biological covariance as a further control for data quality.

In order to assess the presence of cell culture contaminants on the transcriptomic data, we used Kraken to assign taxonomic labels to the reads using the MiniKraken database, which contains all complete bacterial, archaeal, and viral genomes in RefSeq [49]. To assess mycoplasma levels in the cell cultures, we used Bowtie2 to map the data to six different strains of mycoplasma known as frequent contaminants of cell cultures (S4 Fig).

Data processing and differential expression analysis

Reads were aligned to host genomes (S6 Table) using a two-step procedure. A first round of mapping used TopHat2, followed by a second round of mapping using Bowtie2 in an attempt to map the remaining unmapped reads [50]. HTSeq-Count [51] was used to count reads mapping to genes annotated in.gtf files. Genes with <1 read mapping in at least half of the samples were removed prior to differential gene expression analysis using the EdgeR package (Bioconductor) [52,53]. FDR values were calculated using the Benjamini–Hochberg method. MDS and statistical analyses of the data were performed in R.

Generation of orthologous gene clusters

We utilised the Ensembl Compara database [54], combined with our expression data, to generate a table of orthologs with the following associated data for each species: Ensembl ID, Gene name, log2FC, and FDR following IFN treatment. The Compara database provides a thorough account of gene orthology based upon whole genomes available in Ensembl and thus provided us with a standardised approach by which to define phylogenetically the clusters of orthologous genes relative to the chicken taken as an outgroup in the orthology assignment. However, certain gene families relating to innate immunity (for example, the IFITM genes) have undergone lineage-specific expansion, potentially resulting in genes not being annotated and clustered within the database [55]. To account for the misannotation and absence of genes in Compara, we improved the table by manually checking (and annotating, if necessary) genes initially found to be ‘missing’ in either one, two, or three of the 10 species analysed in this study. For this subset of genes, we searched for the presence of an as-yet-unannotated ortholog in the Ensembl genomes using blastn. In cases where a clear ortholog was detected, we included this gene within the appropriate orthogroup. In total, we identified an additional 18 genes that were added to the final table. In cases where a gene was not annotated in the Ensembl genome but a sequence (or predicted sequence) for the homologous species was available in NCBI, we mapped the RNAseq reads to the gene of the homologous species using Bowtie2. The number of reads mapping from each sample was then counted. In cases where an ortholog was present in NCBI from a closely related species, a relaxed Bowtie2 algorithm was first used to map reads to the sequence. The consensus sequence of the resulting contig was then used to count the reads in individual samples using Bowtie2 as above. Differential expression values were then determined using EdgeR having appended the results to the HTseq file. In total, 64 orthologs were added to the table as a result of the orthologous sequence mapping approach. Finally, we modified the .gtf file of sheep to reannotate STAT4 (to become STAT1 and STAT4) and SOCS1, and manually annotated the ZCCHC2 gene in the cow .gtf file.

Verifying species-unique genes

Using the human genome as a gold standard for annotation, we next assessed the authenticity of each “species-specific” ISG (i.e., an ISG present only in one of the 10 species analysed in this study) based upon its current annotation within the genomes. We first generated a subset of species-specific genes by applying an arbitrary cutoff of their differential expression of ≥2 log2FC. This cutoff resulted in a total of 102 genes across eight species. A total of 33 genes (32%) were in this category as a result of misannotation. In particular, the current little brown bat and pig genomes appear to be currently less well annotated compared to other genomes, as a large proportion (73% and 79%, respectively) of seemingly species-unique genes were in reality misannotations.

In silico genome screening

We used the database-integrated genome screening (DIGS) tool [37] to systematically screen the genomes of 111 mammalian species for sequences disclosing homology to 79 of the 90 core^mamm ISGs. The peptide sequences of the human copies of these 79 genes (obtained from BioMart [56]) were used as probes for tBLASTn-based searches of each species genome. Sequences disclosing above-threshold similarity to peptide probes were extracted and classified by BLASTx-based comparison to a reference library. This library contained, for each ISG, peptide sequences of human paralogs and orthologs from selected additional mammal species. The DIGS tool captures screening results in a relational database, wherein they could be interrogated using structured query language. The number of significant matches for each gene was determined using a gene-specific bitscore cutoff. These counts were normalised to the median value across the mammalian genomes screened to account for variation in exon numbers. Because DIGS is based upon sequence-similarity screening, high counts for a particular gene do not necessarily reflect bona fide gene expansion. Where no matches to a given gene were identified and no ortholog had been annotated in Ensembl, we attempted to confirm that deletion had occurred by viewing the corresponding genomic region in the UCSC and Genomicus genome browsers [57,58] and by comparing alignments of the orthologous genomic regions derived from species with and without the deleted gene. Deletions could not always be confirmed mainly due to low coverage or relatively poor quality assembly of the available genomes.

dN/dS ratios for vertebrate interferomes

Clusters comprising one-to-one orthologs present in the interferome of each species were extracted and filtered to check for the presence of a gene in species X matching a human Ensembl ID. The human Ensembl ID was then used to query BioMart to extract the corresponding pairwise dN and dS values for species X against the human. Values for differential expression (log2FC) and FDR values for species X were merged with the dN/dS values and histograms of both significant (ISG) and nonsignificant (non-ISG) clusters plotted. dN and dS values were not available in Ensembl for the Large flying fox or the chicken. The overall distributions of dN/dS values for ISGs and non-ISGs were compared using the Kruskal–Wallis rank sum test and Wilcoxon rank sum test.

Webserver

In order to allow mining of our data by the wider research community, we created a web-based public interactive data server, accessible at http://isg.data.cvr.ac.uk. The server hosts a database containing the orthologous ISG clusters studied in this paper. This web tool allows users to search for and download orthologous clusters and the associated experimental results. The search may be based on various search criteria:

1. Search clusters by HUGO Gene Nomenclature Committee (HGNC) gene name or Ensembl ID.

2. Select clusters based on the presence and expression characteristics in particular species.

3. Select clusters by presence or expression characteristics in a specific number of species.

It should be noted that genes can have different aliases to those in Ensembl, and these must be checked if a gene that is not initially present, e.g., MDA5, is present as IFIH1. These aliases are stated in Ensembl.

URLs

The webserver for querying the dataset is available at http://isg.data.cvr.ac.uk/.

The DIGS blast framework is available at http://giffordlabcvr.github.io/DIGS-tool/.

The authors wish to thank the University of Glasgow School of Veterinary Medicine and Pablo Murcia for sources of clinical material, William Kohler, Michael Whitt and Megan Shaw for useful reagents. We additionally thank Megan MacLeod, Paul Garside, Ian Goodfellow, Rick Randall, Simon Babayan, and Chris Boutell for helpful discussions.