Gene Expression Profiles from Disease Discordant Twins Suggest Shared Antiviral Pathways and Viral Exposures among Multiple Systemic Autoimmune Diseases

<div><p>Viral agents are of interest as possible autoimmune triggers due to prior reported associations and widely studied molecular mechanisms of antiviral immune responses in autoimmunity. Here we examined new viral candidates for the initiation and/or promotion of systemic autoimmune diseases (SAID), as well as possible related signaling pathways shared in the pathogenesis of those disorders. RNA isolated from peripheral blood samples from 33 twins discordant for SAID and 33 matched, unrelated healthy controls was analyzed using a custom viral-human gene microarray. Paired comparisons were made among three study groups—probands with SAID, their unaffected twins, and matched, unrelated healthy controls—using statistical and molecular pathway analyses. Probands and unaffected twins differed significantly in the expression of 537 human genes, and 107 of those were associated with viral infections. These 537 differentially expressed human genes participate in overlapping networks of several canonical, biologic pathways relating to antiviral responses and inflammation. Moreover, certain viral genes were expressed at higher levels in probands compared to either unaffected twins or unrelated, healthy controls. Interestingly, viral gene expression levels in unaffected twins appeared intermediate between those of probands and the matched, unrelated healthy controls. Of the viruses with overexpressed viral genes, herpes simplex virus-2 (HSV-2) was the only human viral pathogen identified using four distinct oligonucleotide probes corresponding to three HSV-2 genes associated with different stages of viral infection. Although the effects from immunosuppressive therapy on viral gene expression remain unclear, this exploratory study suggests a new approach to evaluate shared viral agents and antiviral immune responses that may be involved in the development of SAID.</p></div

Four distinct viral oligo probes from three HSV-2 genes that were differentially expressed between probands and unaffected twins, as well as between probands and matched, unrelated healthy controls^*.

<p>*Paired comparisons were made between probands (P) and unaffected twins (U) or matched, unrelated healthy controls (C). Statistically significant viral probes (fold change, FC >1.5, <i>q</i> <0.05) were selected.</p><p>Four distinct viral oligo probes from three HSV-2 genes that were differentially expressed between probands and unaffected twins, as well as between probands and matched, unrelated healthy controls^{<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0142486#t003fn001" target="_blank">*</a>}.</p

Comparisons of differential gene expression values determined by relative quantitative-polymerase chain reaction and microarray analyses^*.

<p>* Fold change (FC) values indicating relative expression levels of two HSV-2 genes (<i>UL36</i> and <i>RS1</i>) between two study groups: probands (P) and matched, unrelated healthy controls (C). qPCR, quantitative polymerase chain reaction; qRT-PCR, quantitative reverse-transcription PCR; Positive, the affected probands with probe values above the proposed normal range (mean plus three standard deviation of control group); Negative, the affected probands with probe values within the proposed normal range.</p><p>Comparisons of differential gene expression values determined by relative quantitative-polymerase chain reaction and microarray analyses^{<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0142486#t004fn001" target="_blank">*</a>}.</p

Dot plots of three HSV-2 gene expression levels in SAID affected probands, unaffected twins, and matched, unrelated healthy controls.

<p>The individual normalized signal data for <i>UL19</i> (A, probe 1 in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0142486#pone.0142486.t003" target="_blank">Table 3</a>), <i>UL36</i> (B) and <i>RS1</i> (C) were plotted according to the three study groups: probands, unaffected twins, and matched, unrelated healthy controls. Statistically significant <i>p</i> values between any two groups are shown (*) and were calculated using one-way ANOVA with corrections for multiple comparisons. The proposed “normal range” values were derived from the mean plus three standard deviations (SD, dash line) of healthy control expression values.</p

Ingenuity pathway analysis of differentially expressed human genes between disease discordant twins.

<p>The analysis compared 789 oligo probes corresponding to 537 human genes significantly differentially expressed (fold change >1.5, <i>q</i> <0.05) between twins discordant for SAID. The upper horizontal axis (blue bars) describe the association of the data set with a given pathway (-log (<i>p</i> value)). The cutoff threshold value (defined as <i>p</i> = 0.05) is shown by the vertical red line. The ratio of the number of genes from the data set that map to a given pathway divided by the total number of molecules that comprise the pathway is shown on the horizontal axis (green triangles). Previously reported viral-related signaling pathways are in bold print.</p

IFN-regulated genes that were significantly differentially expressed between SAID probands and their unaffected twins^*

<p>*Listed are the 41 IFN-regulated genes of 107 viral infection-related genes identified by the INTERFEROME program that were statistically significant (<i>q</i> < 0.05) and had a fold change >1.5 between disease discordant twins. Fold change values indicate increase (positive) or decreased (negative) levels of gene expression in probands related to unaffected twins.</p><p>IFN-regulated genes that were significantly differentially expressed between SAID probands and their unaffected twins^{<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0142486#t002fn001" target="_blank">*</a>}</p

Flow chart of the experimental design and data analyses.

<p>Flow chart of the experimental design and data analyses.</p

Graphical summary of variation in a presumptive regulatory VNTR containing region.

<p>Panel <b>A</b> shows an overview of approximately 2 Mb locus on 21q, around four genes encoding functionally related receptors; <i>IFNAR2</i>, <i>IL10Rβ</i>, <i>IFNAR1</i> and <i>IFNGR2</i>. Panel <b>B</b> is zooming on the position of the hypervariable region (HVR, red box), which is located approximately 4 kb upstream from the transcription start site of the <i>IFNAR1</i> gene and is flanked by CpG-islands (green boxes). The last three and the first three exons of <i>IL10Rβ</i>, and <i>IFNAR1</i>, respectively, are shown as grey boxes. Panel <b>C</b> is showing the size and position of HVR according to the most common allele (HVR1098, see below panel D) in relation to the CpG island. Positions of PCR and sequencing primers used in the analysis of the locus are also displayed. Yellow boxes indicate the position of the non-repetitive anchor 1 (A1) and anchor 2 (A2) sequences, that are immediately flanking the repeated segments and were used for alignments of sequence reads. Panel <b>D</b> shows a summary of eight HVR-alleles from the studied samples, which were identified based on Sanger sequencing results of PCR fragments sub-cloned in plasmids. The displayed alleles are ordered from longest to shortest according to size from anchor 1 (A1) to anchor 2 (A2) sequences. Summary of sizes for all 14 different HVR-alleles is shown in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0067752#pone-0067752-t001" target="_blank">Table 1</a>. Sizes of fragments (in base pairs) are given between non-repetitive A1 and A2 sequences and between primers p1 and p8, which were used for PCR amplification from genomic DNA. Asterisk (*) indicates the most frequent allele (HVR1098), which is in agreement with the reference sequence according to NCBI sequence build 36.3. The allele frequency shown here is taking into account only the nine alleles, where the entire sequence could be unequivocally determined using Sanger sequencing. The most common variation encompasses the variable number of 32 base pair segments; i.e. indel 2, indel 3, indel 4, and indel 5. The latter indel 5 is composed of 6 repeated 32 base pair segments (HVR1066). However, there are also indels containing shorther segments; e.g. indel 1, indel 6 and indel 7. Panel <b>E</b> illustrates the positions of two of the four probes from Illumina beadchips, which are aligned onto the NCBI reference sequence for this locus (top sequence with an asterisk, representing HVR1098). The two probes shown here are from Illumina 610 SNP array; cnvi0010761 (green) and cnvi0010759 (blue). All four Illumina probes from <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0067752#pone-0067752-g001" target="_blank">Figure 1</a>, which were used for initial identification of variation in this region are located within hypervariable region. As shown here for two of these four probes, the probeA sequences (as called by Illumina and used for capturing of genomic DNA on beadchips) are shifted only by two bases. The core 32 bp repeat motif is shown in brackets.</p

Size and distribution of 14 HVR-alleles identified by sequencing and gel electrophoresis.

<p>The •/× indicate that the size of the allele was determined by both agarose gel images and sequencing, whereas filled circles (•) denote that the allele size was estimated from agarose gel images.</p><p>The two most common alleles (HVR1098 and HVR1700) are highlighted in bold and underlined text.</p><p>The sample indicated by a single asterisk (*) are from breast cancer patients. BL, PT and UM indicate peripheral blood DNA, primary breast tumor and healthy morphologically normal breast tissue from a patient affected with breast cancer, respectively.</p><p>The samples indicated by two asterisks (**) are monozygotic twin pairs.</p

Indications of post-zygotic copy number variation in a region between <i>IL10Rβ</i> and <i>IFNAR1</i>.

<p>Results from eight Illumina genotyping experiments are shown using blood DNA from two pairs of monozygotic twins (panel A for twin 012_01 versus co-twin 012_02, and in panel B for twin 159201 versus co-twin 159202) and two unrelated individuals, where two different tissues were analyzed from each subject (Panels C and D; subjects ML36 and SK58, respectively). Abbreviations BL, PT and UM indicate peripheral blood DNA, primary breast tumor and healthy morphologically normal breast tissue from a patient affected with breast cancer, respectively. Illumina 610 or 660W SNP arrays were used, which also contain so called “intensity only probes” (often with cvni-prefix), only useful for copy number analyses. Therefore, only Log R Ratio (LRR) windows of Illumina experiments are shown here, since the B Allele Frequency (BAF) values are not informative for this type of probes. LRR values below and above zero suggest a deletion or a gain, respectively. The four array probes showing variation between the studied samples are labeled as red dots in yellow fields.</p