RNA-seq analysis of impact of PNN on gene expression and alternative splicing in corneal epithelial cells

Cell lines

These studies employed an SV40-immortalized human corneal epithelial cell line (HCET), with properties similar to normal corneal epithelial cells, such as cornea-specific keratins and abundant aldehyde dehydrogenase activity [1]. Cell lines were authenticated using short tandem repeat (STR) analysis as described in 2012 in ANSI Standard (ASN-0002) Authentication of Human Cell Lines: Standardization of STR Profiling by the ATCC Standards Development Organization (SDO) and by Capes-Davis et al. [27]. Seventeen STR loci plus the gender-determining locus, amelogenin, were amplified using the commercially available PowerPlex® 18D Kit from Promega, Madison, WI. The cell line sample was processed using the ABI Prism® 3500×l Genetic Analyzer. Data were analyzed using GeneMapper® ID-X v1.2 software (Applied Biosystems, Inc., Grand Island, NY). Appropriate positive and negative controls were run and confirmed for each sample submitted. The HCET cell line used in this study (ATCC FTA bar code STRA0571) was found to be human and a match to two other HCET cell lines, RCB1384 and RCB2280, in the DSMZ STR database (Appendix 1).

For doxycycline-inducible PNN-knockdown HCET cells, human TRIPZ shRNAmir clones (V2THS_170187) were purchased from Open Biosystems (Lafayette, CO). HCET cells transfected with either shRNAmir clone were then selected with 10 μg/ml of puromycin (#61-385-RA, Cellgro, Manassas, VA). To induce shRNA expression, the cells were treated with doxycycline (Doxy) at a concentration of 1 μg/ml with a daily change of fresh doxycycline medium. Using western blots, we detected the typical knockdown at 48 h post-Doxy, D3447, Sigma Aldrich, St. Louis, MO in the range of 46–75% reduction.

RNA-seq library preparation

RNA-seq libraries were prepared beginning with total RNA (RIN 8.60–10) extracted from control and PNN knockdown HCET cells, using RNeasy Plus Mini Kit (Qiagen, Valencia, CA). Before the RNA-seq library was constructed, rRNA was removed from the total RNA samples using an Illumina Ribo-Zero kit. Following rRNA depletion, RNA-seq library construction used ScriptSeq Complete v2 (Illumina, San Diego, CA). The ScriptSeq protocol used sequence-tagged random hexamer priming to produce RNA-seq libraries that retained the identity of the original transcribed strand.

DNA sequence determination

Each library was barcoded and then quantified. The libraries were then pooled into a single pool, and the pool was run on four independent lanes [28] of the Illumina NextSeq 500. Quality control was performed using a combination of FASTQC and Lab Scripts that count the number of duplicate reads, the number of reads with an adaptor, and the number of reads with homopolymer sequences. The results indicated the data were of good quality and all of the libraries were comparable. The native Illumina bcl DNA sequence output format was converted to the fastq format and demultiplexed by index reads using bcl2fastq version 2.15 (Illumina). The average number of reads per sample was about 40 million (average 39,336,132.75 reads/sample), and the distribution of read numbers is given in Appendix 2. All of the reads have been deposited in the SRA (accession # SRP063860).

Read alignments for viewing in UCSC Genome Browser

Paired-end reads acquired as described were processed using Trimmomatic version 0.32 [29] to remove oligonucleotide adaptors used during RNA-seq library construction and to trim low-quality sequences. DNA sequence reads that had a minimum length of 50 nt were retained for downstream processing. After filtering, an average of 35 million paired reads were retained per sample. Paired reads were mapped to the GRCh37 human reference genome (Genome) using STAR aligner version 2.0.4c [30] with default mapping parameters. Mapped paired reads were filtered to remove duplicates using the Genome Analysis Toolkit [31]. SAMtools [32] was used to sort mapped reads by chromosomal coordinate, create, and index human chromosome-specific map files, and convert the file format from SAM to BAM. The bam2wig.py Python script from RSeQC [33] was used to create chromosome-specific wig-format files, which were uploaded for viewing in the University of California, Santa Cruz (UCSC) Genome Browser.

The stranded RNA-seq data were aligned to the GRCh37/hg19 version (release 73) of the complete human genome using the BWA-MEM algorithm [34]. Genomic coverage was summarized using AceView annotations [35]. There were 678,664 exons in this annotation, and the genome positions of approximately 68% of these exons overlapped with at least one other exon (e.g., alternative donor/acceptor sites). Overlapping exons in the same direction were combined, and the minimum start and maximum end positions were used to estimate expression in 322,096 unique exonic regions. For ease, these exonic regions are referred to as exons throughout the text. Exons were determined to be “expressed” if they had coverage in at least three samples per group.

Annotating alternative splicing events

We created a complete, annotated catalog of splicing events from genomic feature databases (GFF files) using gene, transcript, and exon-level references for RNA-seq alignments and coverage calculations. We used the latest (2010) release of AceView gene models for the hg19 human genome build [35,36]. A list of possible exon–exon junctions was generated from all logical combinations of exon pairs within a gene in a 5′-to-3′ manner.

For AS event expression analysis, RNA-seq data were aligned to the sequences of all exon junctions and retained intron events using Bowtie [37]. AS event coverage was summarized using the annotations defined above. AS events were deemed “expressed” if they had more than ten aligned reads in at least three samples per group.

We fit the model Y_ij=μ + t_j + ϵ_ij, where Y_ij is the log₂-transformed average depth per nucleotide (APN) for exon or AS event i, sample j, t_j is the treatment status of sample j, and ϵ_ij is the random error. A false discovery rate (FDR) was used to correct results for multiple testing [38], and a level of 0.05 was considered to be significant.

Gene ontology analysis

Gene Ontology was analyzed and demonstrated using GOrilla a tool for identifying and visualizing enriched GO terms in ranked lists of genes. GOrilla searches for enriched GO terms that appear densely at the top of a ranked list of genes. The system utilizes 13,033 genes that are associated with a GO term GO database and other sources. Graphical representation of the GO results was accomplished with the Graphviz layout programs, which take descriptions of graphs in a simple text language and create diagrams in useful formats.

Droplet digital PCR

Total RNA was isolated from cultured HCET cells, and semiquantitative reverse transcription (RT)-PCR was performed as described previously [39,40], with RNeasy Plus Mini Kit (#74134, Qiagen, Valencia, CA) and treated with RNase-free DNase I (#79254, Qiagen). One microgram of total RNA was reverse transcribed with the SuperScript III First-Strand Synthesis kit (#18080051, Invitrogen, Carlsbad, CA) using oligo-dT primers. The droplet digital PCR (ddPCR) assays were performed according to the Bio-Rad User Guide. Briefly, each of the 20 μl reactions contained 10 μl ddPCR Supermix (Bio-Rad, Hercules, CA), 250 nM gene-specific primers, and 0.5 μl of the cDNA sample (about 5 ng). Primer sequences for PNN, MMP1 (Gene ID:4312; OMIM 120353), MMP9 (Gene ID:4318; OMIM 120361), MMP13 (Gene ID:4322; OMIM 600108), VEGF (Gene ID:7422; OMIM 192240), and GAPDH are listed in Appendix 3. Each reaction was mixed with 20 μl of Droplet Generation Oil (Bio-Rad), partitioned into about 20,000 droplets in a QX200 Droplet Generator (Bio-Rad), transferred to 96-well plates (Eppendorf, Hauppauge, NY), and sealed. PCR was performed in a PTC-100 Thermal Cycler (Bio-Rad) with the following cycling conditions: 1 × (95 °C for 5 min); 40 cycles: (95 °C for 30 s, 60 °C for 60 s); 1 × (4 °C for 5 min; 90 °C for 5 min); 4 °C (hold). The fluorescence intensity of individual droplets was measured with the QX200 Droplet Reader (Bio-Rad). The data analysis was performed with QuantaSoft droplet reader software (Bio-Rad). Positive and negative droplet populations were detected automatically. The target mRNA concentrations were calculated using the Poisson statistics and background-corrected based on the no template control data. The absolute transcript levels were initially computed in [copies/μl PCR] and in [copies/ng cDNA]. In the latter case, the correction for the input cDNA amount was applied to each sample. The data are presented as expression levels normalized to GAPDH expression, p value, and log₂ fold change, for MMP1, MMP9, MMP13, VEGF-A, and PNN.

Immunofluorescence

Doxycycline-inducible PNN-knockdown HCET cells were plated on coverslips in six-well dishes at 2×10⁵ cells per well in a Dulbecco's Modification Eagle's Medium/Ham's F-12 50/50 mix, (#10-090-CV, Cellgro, Manassas, VA); 5% fetal bovine serum [FBS, #35-011-CV, Cellgro], 5 μg/ml insulin, 0.1 μg/ml cholera toxin (C8052, Sigma Aldrich), 10 ng/ml human epidermal growth factor [hEGF, PHG0311, Gibco, Frederick, MD], 0.5% dimethyl sulfoxide [DMSO, D2650, Sigma Aldrich]). After 48 h treatment in media with and without doxycycline (1 μg/ml), cells were fixed in freshly prepared 2% paraformaldehyde (D2650, Sigma Aldrich) in PBS (1X; 137 mM NaCl, 2.7 mM KCl, 8 mM Na₂HPO₄, 1.5 mM KH₂PO₄, pH 7.4) for 1 h, quenched in 50 mM ammonium chloride (A-661, Fisher Scientific), 0.1%Tween-20 (BP-337, Fisher Scientific) in PBS for 10 min, rinsed and incubated for 2 h at room temperature in 1:200 rabbit anti-MMP9 antibody (#13667, Cell Signaling, Danvers, MA). After rinsing with PBS, the cells were incubated in 1:200 fluorescein isothiocyanate (FITC)-labeled goat anti-rabbit antibody (#F6005, Sigma-Aldrich, St. Louis, MO) for 1 h. The coverslips were mounted in VectaShield with 4',6-diamidino-2-phenylindole (DAPI; #H1200, Vector Labs, Burlingame, CA) and examined using a Leica inverted fluorescent microscope (Model DM-IRBE, Leica Microsystems, Buffalo Grove, IL) and OpenLab (Improvision, Coventry, England) software.

RNA-seq and GO analysis

To evaluate the effect of PNN on mRNA expression and RNA processing on a global scale, we performed mRNA-seq experiments in HCET cells in which the levels of PNN were reduced with specific siRNA. RNA-seq-based transcriptional profiling was employed to compare Wt-HCET cells with doxycycline-inducible PNN-knockdown HCET cells (50% knockdown). Data were deposited in GEO under accession number GSE73060. The RNA-seq profiles of 36,275 transcripts were obtained, out of which 90 are significant at P-adj (p<0.01). Of these, 61 genes exhibited >1.0 log₂ fold change, and 28 genes showed <–1.0 log₂ fold change. The P-adj is an adjusted p value, taking into account the FDR, which is the proportion of the total rejections that are false positives [41]. The FDR has been widely used for multiple test corrections in various applications, such as when measuring thousands of variables (e.g., gene expression levels) from a small sample set [42]. Table 1 and Table 2 depict the top 20 significant genes with the largest fold changes (supplementary data provided in Appendix 4 and Appendix 5). GO analyses were performed with Gorilla, a tool for identifying and visualizing enriched GO terms in ranked lists of genes [43,44]. Table 3 and Table 4 show the GO analyses for upregulated and downregulated genes, respectively. The upregulated genes included a large proportion of genes that are associated with enhanced cell migration and ECM remodeling processes, including MMPs, ADAMs, HAS2, LAMA3, and UNC5C (Table 3, Figure 1 and supplementary data Appendix 6). In addition, the upregulated genes included members of the CXCR chemokine receptor binding family, metalloendopeptidase activities, and extracellular matrix organization and structure. GO analysis revealed that a significant subset of genes, which impact cellular differentiation, were included in those that exhibited downregulation, such as IGFBP5, FGD3, FGFR2, PAX6, RARG, and SOX10 (Table 4). Taken together, these downregulated genes may contribute to the loss in corneal epithelial differentiation seen in downregulation and knockouts of PNN [3-6,39].

The upregulation of MMPs and VEGF mRNA, observed via RNA-seq subsequent to PNN knockdown, was confirmed with quantitative PCR, using droplet digital technology [45]. The downregulation of PNN (log₂ fold change, −1.93) was noted at 48 h post-doxycycline induction of PNN-shRNA, as was significant upregulation of MMP1 (log₂ fold change, 1.69), MMP9 (log₂ fold change, 0.78), MMP13 (log₂ fold change, 2.0) and VEGF-A (log₂ fold change, 0.49) at 48 h post-Doxy (Table 5). Immunofluorescence microscopy confirmed the upregulation of MMP9 protein in HCETs 48 h after Doxy treatment (Figure 2). Interestingly, when downregulation of PNN was less (log₂ fold change, −0.82), there was no observed significant upregulation of MMP1 (log₂ fold change, −0.02) or MMP9 (log₂ fold change, 0.23), but significant changes were still observed for MMP13 (log₂ fold change, 1.73) and VEGF-A (log₂ fold change, 0.43) at 48 h post-Doxy (supplementary data in Appendix 7). Although outside the scope of this manuscript, we suggest that these data may reflect a dosage effect of PNN knockdown. Future experiments will explore the dosage response to PNN knockdown, as well as the differential sensitivity seen across the MMPs.

Annotating alternative splicing events

The evaluation of the RNA-seq analyses of the PNN-knockdown cells revealed the remarkable dynamics of RNA processing of the expressed genes coding and non-coding RNAs. The results are broadly consistent with previous observations from hGlue3_0 transcriptome array analyses cGH arrays [39,40]. One of the main advantages of RNA-seq is the ability to identify and quantify AS events like exon–exon junctions and retained introns. We created a complete, annotated catalog of all possible splicing events using genomic feature databases (GFF files) used when building gene, transcript, and exon-level references used for RNA-seq alignments and coverage calculations. Annotations for each event were classified based on four main alternative splicing event types: alternative donor, alternative acceptor, skipped exon, and retained intron.

For this study, we used the latest (2010) release of AceView gene models for the hg19 human genome build, due to the noted higher accuracy of AceView gene models as well as the higher abundance of gene models per gene compared with Ensembl and RefSeq [35,36]. Our annotation method described here is not specific to the AceView annotations. The scripts to build the AS catalog are available at GitHub and can be used for any existing database of genomic feature annotations.

A list of possible exon–exon junctions were generated from all logical combinations of exon pairs within gene in a 5′-to-3′ manner. Junctions were not generated if they were not biologically plausible (e.g., an exon donor site to an earlier exon acceptor site or a junction between two overlapping exons), nor were junctions created between genes. To facilitate the classification of AS events, overlapping exons were grouped together, and the longest exon in each group was arbitrarily assigned as a reference exon and used to define alternative donor and acceptor sites within each exon group. There were several occurrences where two non-overlapping exons overlapped with another exon. In this case, the longest exon per group was used as the reference, and the junction between the two internal exons was created. We classified junctions as exon skipping if there was at least one exonic region between the exons two exons forming a junction. For each exon-skipping junction, we also generated a list of possible skipped exons. Skipped exons were included on this list if they were unambiguously skipped (i.e., exons overlapping the donor or acceptor exon were not included). We also annotated junctions that were part of at least one known transcript. Intron retention (IR) events were generated by using the donor site of an exon and extending the sequence into the neighboring intron. For groups of overlapping exons, we generated IR events only from the 3′-most exon to ensure that only a uniquely intronic sequence was used for the intronic segment of the event. In total, there were 6,390,703 possible AS junctions and 232,249 IR events annotated in this database: 5,837,064 junctions were annotated as an exon-skipping junction, 1,240,520 junctions were classified as having an alternative donor, 1,257,096 junctions were classified as having an alternative acceptor, and 903,513 junctions were classified as having an alternative donor and an alternative acceptor. In addition, 292,753 junctions were annotated to at least one known transcript. The database is available on Github.

Figure 3 summarizes the impact of PNN knockdown on alternative splicing in HCET cells. We plotted the mean expression (as log₂ APN) of AS events detected in PNN knockdown and control cells (58,817 of 80,813 total detected splicing events) and noted that the expression of AS events was generally higher in the PNN-knockdown cells compared to controls (Figure 3). In total, we observed that 5,063 AS events were differentially expressed, with 4,683 AS events (92.5%) significantly higher after PNN knockdown and the remaining 380 significant AS events higher (7.5%) in the control cells, indicating that reduction of PNN expression results in upregulation of AS events (p<0.05; Figure 3). For alternative splicing events detected in both conditions and differentially expressed, all types of AS events were more likely to be upregulated in PNN knockdown HCET cells: 92.5% of all AS events, 91.8% of exon-skipping events (635 AS events were unregulated of the 692 that were significantly different), 94.1% of alternative donor sites (531 AS events were unregulated of the 564 that were significantly different), 93.0% of alternative receptor sites (490 AS events were unregulated of the 527 that were significantly different), 93.7% of AS events involving alternative donors and acceptors (554 AS events were unregulated of the 591 that were significantly different), and 62.9% of intron retention events (95 IR events were unregulated of the 151 that were significantly different). The detailed list of differentially expressed events is presented in the supplementary data (Table 6 and supplementary data in Appendix 8). The knockdown of PNN resulted in a dramatic increase in the number of splicing events, suggesting that altering levels of PNN has a broad influence on splicing. This widespread impact has been observed after the knockdown of other members of the exon junction complex [25]. Interestingly, genes that were upregulated by PNN knockdown displayed a much greater increase in altered splicing events than genes not changed or downregulated. These data speak to the coupling of splicing to transcription.

Transcript assembly and quantification with RNA-seq revealed bone fide isoform switching and the generation of unannotated transcripts subsequent to PNN knockdown. Table 7 (more complete list, supplementary data (Appendix 9) reveals the observed splicing events (p<0.05) associated with splicing variants after PNN knockdown. The upregulated variant and log fold change are listed on the left, while the downregulated variant and fold change are listed on the right in the table. Interestingly, in addition to variants in the ADAM family (membrane disintegrin and metalloproteinase), many splicing variations occur in the abundant solute carrier (SLC) group of membrane transport proteins after PNN knockdown. The observed high correlation with SLC group members may reflect the abundance of SLC mRNA and proteins in mammals.

RNA splicing occurs cotranscriptionally, and we hypothesize that PNN may serve as a functional bridge between transcription and RNA splicing. Therefore, we examined the correlation of change in gene expression and the occurrence of altered RNA splicing. Table 8 and Table 9 and supplementary data (Appendix 10) show the splicing events (p<0.05) that are correlated with upregulated and downregulated genes, respectively. Most interestingly, a large proportion of the upregulated genes, which GO analyses indicate are involved in enhanced cell migration such as MMPs, ADAMs, SERPINs, ECM, and matrix receptors, also demonstrated altered RNA splicing. Although the list of downregulated genes with correlated alteration RNA events is smaller, of significant note is the change in splicing events of two unexpected cadherins in PNN knockdown, endothelial cadherin 5 and minor cadherin 11. Although the roles of these cadherins in the corneal epithelial context is uncertain, these data are similar to that we reported previously for the coregulation of transcription and splicing of E-cadherin by PNN [46] and may indicate a more global role of PNN in RNA processing on cadherin genes.