RNA-Seq for enrichment and analysis of IRF5 transcript expression in SLE

Abstract

Polymorphisms in the interferon regulatory factor 5 (IRF5) gene have been consistently replicated and shown to confer risk for or protection from the development of systemic lupus erythematosus (SLE). IRF5 expression is significantly upregulated in SLE patients and upregulation associates with IRF5-SLE risk haplotypes. IRF5 alternative splicing has also been shown to be elevated in SLE patients. Given that human IRF5 exists as multiple alternatively spliced transcripts with distinct function(s), it is important to determine whether the IRF5 transcript profile expressed in healthy donor immune cells is different from that expressed in SLE patients. Moreover, it is not currently known whether an IRF5-SLE risk haplotype defines the profile of IRF5 transcripts expressed. Using standard molecular cloning techniques, we identified and isolated 14 new differentially spliced IRF5 transcript variants from purified monocytes of healthy donors and SLE patients to generate an IRF5 variant transcriptome. Next-generation sequencing was then used to perform in-depth and quantitative analysis of full-length IRF5 transcript expression in primary immune cells of SLE patients and healthy donors by next-generation sequencing. Evidence for additional alternatively spliced transcripts was obtained from de novo junction discovery. Data from these studies support the overall complexity of IRF5 alternative splicing in SLE. Results from next-generation sequencing correlated with cloning and gave similar abundance rankings in SLE patients thus supporting the use of this new technology for in-depth single gene transcript profiling. Results from this study provide the first proof that 1) SLE patients express an IRF5 transcript signature that is distinct from healthy donors, 2) an IRF5-SLE risk haplotype defines the top four most abundant IRF5 transcripts expressed in SLE patients, and 3) an IRF5 transcript signature enables clustering of SLE patients with the H2 risk haplotype.

Figure 1. Identification of IRF5 transcript variants expressed in immune cells of SLE patients and healthy donors.

A, A profile of IRF5 transcripts expressed in purified Mo of healthy donors (n = 3) and SLE patients (n = 6) was generated by molecular cloning and sequencing. PCR amplification was from non-coding Ex1a through Exon 9. At least 20 clones from each donor were obtained. Data are represented as a pooled % of total clones after calculating % of clones per variant per donor. Statistical analysis was performed using the two-way Anova with Bonferroni post-test. Known transcript variants V1, V4, V5 and V6 are depicted with novel variants (NV1-NV14). B, Scheme illustrating all IRF5 transcripts cloned from purified Mo of healthy donors and SLE patients in a display generated using IGV software. Transcripts are comprised of exons (blue rectangles) spliced from the full IRF5 genomic sequence (horizontal lines with arrowheads). Transcript NV12 is not shown since IGV software displays each transcript using the IRF5 genomic sequence as a reference and is thus unable to represent duplications; NV12 contains a 12-bp internal duplication of a portion of exon 8 but is otherwise identical to transcript NV1. C, Diagram of exon 6 divided into portions 6A-6D as found in alternatively spliced IRF5 transcripts. Ex6a and/or Ex6c (white rectangles) are deleted in several of the transcripts. The intronic insertion of transcript NV10 is shown as well.

Figure 2. Pile-up view of reads from deep sequencing of PCR-amplified IRF5 fragments.

A, Pile-up view of IRF5 sequencing reads mapped to the IRF5 genomic sequence was generated with IGV. Peak heights represent the relative frequency of mappings to the indicated regions - hills reflect a high mapping frequency of reads to a particular region and valleys indicate low mapping frequency or deletions (example shown with star); arrows point to areas of reads aligning to previously unannotated exons. Results shown are from 4 independent healthy donors and 5 independent SLE patients. B, Reads from healthy donors and SLE patients that mapped to exon 7, shown by the star in A above, and reflect true internal deletions of exon 7 are shown. Junction reads are aligned by position and labeled as to whether it was detected in a healthy or SLE sample; the exact number of reads per sample that contained the junction is shown in parentheses.

Figure 3. Workflow for enrichment and analysis of single gene profiling by RNA-Seq.

Following PCR amplification of RNA of a single gene, cDNA is sequenced by next-generation sequencing and reads are mapped to sequences of alternatively spliced transcripts to obtain expression estimates and to compare transcript expression profiles across samples.

Figure 4. Expression estimates of IRF5 alternatively spliced transcripts in healthy donors and SLE patients.

A, Gibbs expression estimates were obtained with MMSeq and converted to raw read counts, then normalized to total read counts. Normalized read counts for each transcript in 5 SLE patients and a representative healthy donor (#1) are plotted; read counts were compared between the two sets using a student's t-test. *p<0.05; **p<0.01. B, Scatter plot showing the correlation between MMSeq-derived expression estimates and true (assigned) expression levels of the 18 transcript variants. R² = 0.6532 at alpha = 0.05 and p<0.0001 indicating a highly significant correlation.

Figure 5. Correlation between Gibbs expression estimates and clone count.

Gibbs expression estimates and Monte Carlo standard errors were generated by MMSeq. Bar graphs show Gibbs expression measures from next-generation sequencing versus clone counts from traditional cloning and sequencing. Data from n = 2 healthy donors and n = 5 SLE patients is shown; at least 20 IRF5-positive clones from each donor were used for this analysis.

Figure 6. Top four most abundant transcripts are shared by SLE patients with the H2 full risk haplotype.

A, Relative abundance ranking of transcripts from 4 healthy donors, 1 SLE patient with the H3 full protective haplotype, and 4 SLE patients with the H2 full risk haplotype is shown. Transcripts were ranked from most to least abundant using MMSeq-generated Gibbs expression estimates; individual transcripts are represented by distinct colors. B, Hierarchial clustering was performed on the nine samples using expression profiles of the 18 IRF5 transcript variants. Height represents the difference between samples. Red values are AU p-values and green values are BP values. Clusters with AU larger than 95% are statistically significant with a p<0.05.

Figure 7. Profiling IRF5 alternative splicing using junction counts and de novo junction discovery.

A, Counts were obtained for the 6 junctions unique to single cloned transcripts by aligning reads from healthy donor (n = 4) and SLE patient samples (n = 5) to 94-bp junction sequences. Proportionally more transcripts containing each junction were found in either healthy donors (open circles) or SLE patients (closed circles); plotted are the p values obtained by the Pearson’s χ²-test. B, De novo junction discovery was performed using TopHat to identify novel alternative splicing events. Representative junctions are displayed using IGV with the IRF5 genomic sequence as a reference. Each pair of red blocks represents one junction, and lines joining two blocks denote spliced-out portions of IRF5. Junctions 1-3 show evidence of splicing to unannotated intronic regions of IRF5; junctions 4–9 have novel inter- and intra-exonic deletions. C, Schematic illustration of junctions in B confirmed by qPCR.

Figure 8. Distinct transactivation potential and stability for IRF5 isoforms associated with the IRF5-SLE H2 risk haplotype.

A, Flag-tagged IRF5 V1, V2/V6, V3/V4, V5, and empty vector (EV) were co-transfected with the indicated luciferase reporters - 5xISRE, IFNα4, IFNβ, IL6, TNFα, or IL12p40, and promoter transactivation determined by dual luciferase assay. Independent experiments were repeated in duplicate wells at least three times; mean ± S.D. are plotted. Statistical analysis was performed by one-way ANOVA with Bonferroni post-test; *p<0.05, **p<0.01, ***p<0.001. Promoter induction by each IRF5 isoform was significantly greater than empty vector (EV) controls in all experiments (astericks not shown). B, Relative half-lives of IRF5 isoforms V1, V2/V6, V3/V4, and V5 were determined by cycloheximide chase experiments (see Materials & Methods). Similar to A, GFP-tagged IRF5 variants were transfected to Hek293T cells; V8 was included as a comparative control lacking the PEST-rich exons 5–7. Representative data from 2–3 independent experiments is shown as mean fluorescence intensity (MFI) of GFP-IRF5 over the indicated time period. Analysis of covariance (ANCOVA) was used to determine if regression lines between isoforms differed significantly.