Deciphering the largest disease-associated transcript isoforms in the human neural retina with advanced long-read sequencing approaches

Tissue collection for PacBio Iso-Seq long-read mRNA sequencing

Human donor eyes, deemed unsuitable for corneal transplantation, were used in this study. Written consent was obtained from the donors’ next of kin in accordance with the guidelines of the Italian National Transplant Centre (Centro Nazionale Trapianti, Rome, Italy). The procurement and processing of these tissues followed Italian law and adhered to the principles of the Declaration of Helsinki and the guidelines of the European Eye Bank Association. Three human neural retinal samples were obtained from nonvisually impaired individuals through the Fondazione Banca degli Occhi del Veneto (Venice, Italy). Before this study, we performed WGS on DNA isolated from these retinal samples to screen for and exclude any known genetic variants associated with inherited retinal disorders. The eyes were enucleated within 2–12 h postmortem, and the retinal extraction followed established protocols (Niyadurupola et al. 2011; Osborne et al. 2016). After cornea removal, the eyeball was dissected at the ora serrata; the iris, lens, and vitreous body were removed before carefully detaching the retina from the sclera and retinal pigment epithelium (RPE) by cutting at the optic nerve head. Subsequently, retinal samples were snap-frozen in cryovials using liquid nitrogen and shipped on dry ice. Detailed information about the donors is presented in Supplemental Table S5.

RNA isolation and PacBio Iso-Seq library preparation following the standard workflow

To preserve the integrity of long mRNA molecules, all samples were handled with care during RNA isolation and library preparation, avoiding vortexing or vigorous pipetting to prevent shearing of mRNA molecules. Total RNA was extracted from the human neural retina samples by adding 500 µL of TRIzol reagent to each sample, which was then homogenized in a 2 mL tube containing a sterile glass bead using a TissueLyser (Qiagen, Aarhus, Denmark) in two cycles at 30 Hz for 30 sec. Following a 5 min incubation at room temperature (RT), 100 µL chloroform was added to the samples which were then mixed and incubated for another 3 min at RT before being centrifuged at 12,000g for 15 min. The resulting aqueous phase was collected and mixed with 1 µL glycogen (5 µg/µL) and an equal volume of isopropanol. This mixture was incubated at 20°C for 75 min followed by centrifugation at 12,000g for 30 min at 4°C. Subsequently, the supernatant was removed and the remaining RNA pellet was dissolved in MQ water and further purified and DNase treated using the Nucleospin RNA Clean-up Kit (Macherey-Nagel, Düren, Germany) according to the manufacturer's instructions. The total isolated RNA quantity was measured with a Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). Additionally, RNA integrity number (RIN) values were determined using a 2100 Bioanalyzer (Aligent Technologies, Santa Clara, CA, USA). The RIN values of all three samples exceeded 7.0 (Supplemental Table S5), and 300 ng of RNA input was used to generate Iso-Seq SMRTbell libraries following the Iso-Seq-Express-Template-Preparation protocol PN 101-763-800 version 2.0 (Pacific Biosciences, CA, USA). These libraries were prepared using the standard workflow, suitable for samples with a predominant transcript size of ∼2 kb, and the SMRTbell library binding kit 2.1 (Pacific Biosciences). The samples were not labeled with a barcode, and 500 ng of cDNA was used for the subsequent steps in the procedure. The quantification and assessment of the SMRTbell library for each sample were conducted using Qubit (Thermo Fisher Scientific) and an Agilent Bioanalyzer 2100 employing HS RNA screentape. The on-plate loading concentration of the final Iso-Seq SMRTbell libraries was set at 80 pM, and sequencing was carried out on a Sequel IIe system (Pacific Biosciences) with a movie time of 24 h.

PacBio Iso-Seq library preparation following adjusted workflow optimized for large transcripts

In total, 300 ng of RNA sample 2 was used to generate the Iso-Seq SMRTbell library using the “Iso-Seq-Express-Template-Preparation” protocol PN 101-763-800 version 2.0 (Pacific Biosciences). The library was prepared using the PacBio long transcript workflow suitable for samples with a transcript size larger than 3 kb using the SMRTbell library binding kit 2.1. To further enhance the enrichment for large transcript sizes, an additional size selection step was included using diluted AMPure PB beads. This was performed after the standard purification of amplified cDNA for long transcripts (pages 6–7 of protocol PN 101-763-800 version 2) and after performing an additional five PCR cycles as outlined in Appendix 1 of PN 101-763-800 (page 11). Next, instead of the recommended ProNex bead purification (page 12), the additional size selection was carried out using a 3.3× ratio of diluted AMPure PB beads according to the PacBio Procedure “Using AMPure PB Beads for Size Selection” protocol PN 101-854-900 version 2.0 (Pacific Biosciences). The concentration and size distribution of the resulting SMRTbell library was evaluated using Qubit (Thermo Fisher Scientific) and an Agilent Bioanalyzer 2100 with HS RNA screentape. The on-plate loading concentration of the final Iso-Seq SMRTbell libraries was set at 120 pM, and sequencing was carried out on a Sequel IIe system (Pacific Biosciences) with a movie time of 30 h.

PacBio long-read Iso-Seq data analysis

Four PacBio long-read RNA Iso-Seq samples (three samples prepared following the PacBio standard workflow, 1 sample prepared following the optimized PacBio long transcript workflow) were analyzed with IsoQuant (Prjibelski et al. 2023). Three separate data sets were created (Fig. 1): the first data set included the IsoQuant analysis of the three standard workflow samples, like we previously performed in Riepe et al. (2024). The second data set combined the analysis of PacBio standard workflow retina samples 1, 2, and 3 with the retina sample prepared following the optimized long transcript workflow. These two data set used CCS3 reads. The third data set consists only of the long transcript workflow sample, analyzed with less strict filtering settings and using CCS0 reads under the assumption that the longest transcripts might not reach CCS3 as frequently as an average-sized transcript. The sequencing data resulting from the PacBio Iso-Seq approaches that were used for IsoQuant analyses are available via the European Genome-phenome Archive (EGA; https://ega-archive.org) with the identifier EGAD50000000720. Detailed properties of the IsoQuant analyses can be found in the Supplemental Code document.

Data analysis and transcript visualization

RStudio software (v4.3.2) (R Core Team 2021) was used to obtain read and transcript counts from the IsoQuant output files. ggtranscript package (Gustavsson et al. 2022) was used to visualize transcripts from the GTF output files. ORFs were predicted with SnapGene. Details and R scripts are posted on GitHub (https://github.com/erikdevrieze/USH_retina), and can be found in the Supplemental Code document.

Targeted Iso-Seq for USH2A and ADGRV1 transcripts using the Samplix Xdrop Sort

The Samplix Xdrop Sort (Samplix ApS, Birkerød, Denmark) was employed for a targeted Iso-Seq analysis, enriched for USH2A and ADGRV1 transcripts. Of each retina sample, 300 ng of RNA was used for cDNA synthesis using the NEBNext Single Cell/Low Input cDNA Synthesis kit (New England Biolabs). For this purpose, the “Preparing Iso-Seq libraries using SMRTbell prep kit 3.0” protocol 102-396-000 version 2.0 (Pacific Biosciences) was followed until the cDNA amplification step. The synthesized cDNA of all three samples was pooled and the final DNA concentration of the resulting pool was determined using the Qubit with the single-stranded DNA kit (Thermo Fisher Scientific) according to the standard workflow. To perform the indirect target enrichment, single cDNA molecules were packaged in double emulsion (DE) droplets together with detection sequence primers using the Samplix Xdrop Sort. The detection sequence primers were designed to target a 100–150 bp region located in either the 5′, middle or 3′ region of full-length USH2A (ENST00000307340.8) and ADGRV1 (ENST00000405460.9) transcripts, respectively, as to increase the likelihood of capturing full-length transcripts as well as known and novel shorter isoforms. Primer sequences and targets are provided in Supplemental Table S6.

The encapsulation process was performed separately for each of the six dPCR primer sets to ensure the targeting of a single region during enrichment. Samples were run in technical duplicates, and two positive controls for amplification and sorting provided by Samplix were included on each Xdrop DE20 Cartridge (Samplix ApS). The input concentration was 2.8 ng of pooled cDNA, to which 20 µL of DE PCR mix (Samplix ApS), 0.8 µL of the respective forward and reverse dPCR primers, and a volume of nuclease-free water was added, resulting in a sample mix with a total volume of 40 µL per sample. For the positive controls, the Samplix positive control DNA and dPCR primers were employed. To ensure proper DE droplet production, the sample and additional reagents were loaded onto the Xdrop DE20 Cartridge in the following order: Firstly, 300 µL of DE PCR buffer (Samplix ApS) (diluted in a 1:1 ratio using nuclease-free water) was loaded into well #A of the cartridge. Next, 40 µL of diluted DE PCR buffer was loaded on the shelf in well #D. Subsequently, 40 µL of the sample mix was pipetted into well #C. Lastly, 100 µL of DE Droplet Oil (Samplix ApS) was loaded into well #B. The cartridge was sealed using a gasket and placed into the Xdrop Sort, after which the DE20 droplet production program was run. Afterward, the DE20 droplets were collected from well #D and dispensed into four equal aliquots per sample in PCR-tubes, which were then placed in a thermal cycler in which the following program was run: 30°C for 5 sec, 94°C for 3 min, 40 cycles of 94°C for 3 sec followed by 65°C for 30 sec, ending in a hold at 4°C. Amplification takes place exclusively in droplets containing cDNA molecules in which the detection sequence targeted by the dPCR primers is present. The encapsulated DNA can be stained using an intercalating dye, resulting in a stronger fluorescent signal in the DE20 droplets containing a target transcript, which enables the sorting of these droplets using the Xdrop Sort. To do so, the four aliquots of each sample were pooled again and the aqueous phase was removed before adding 1 mL of DE staining buffer (Samplix ApS). Samples were incubated for a duration of 15 min at room temperature in the absence of light. An Xdrop DE20 Sort Cartridge (Samplix ApS) was sealed with the accompanying sorting foil, ensuring a tight seal around all wells. To prevent residual fluorescence from neighboring lanes, the sorting process was performed in two separate steps: first sorting of the uneven lanes took place, followed by a second sorting run for the even lanes, taking along one of the previously mentioned positive controls in each step. Using the provided lane opener, the lanes to be used in the sorting run were punched open and the cartridge was loaded with the sample and additional reagents in the following order: 5 µL of Xdrop Blank Oil Droplets (Samplix ApS) were loaded in well #Out, 600 µL of DE Sorting Buffer (Samplix ApS) (diluted in a 1:1 ratio using nuclease-free water) was loaded in well #B1, 300 µL of diluted DE Sorting Buffer was loaded in well #B2 and lastly 300 µL of sample in staining buffer was pipetted in well #In. The cartridge was sealed using a gasket and left to settle in the Xdrop Sort for 5 min, after which the Sorting program was run for the selected lanes. Thresholds for detection and sorting were determined based on the displayed signal of droplets and background fluorescence. As these values differed per sample, the guidelines for picking a threshold as established in the Xdrop Sort Manual (Samplix ApS) were consulted. Directly after sorting, droplets were removed from the #Out well of the cartridge. After incubation for 5 min at room temperature, the aqueous phase was removed in order to start the process of breaking the droplets to prepare the DNA for the subsequent MDA. The droplets were washed twice using 200 µL of Droplet Sorting Wash Buffer (Samplix ApS), after which all the wash buffer was carefully removed. Subsequently, 20 µL of Droplet Break Solution (Samplix ApS) was added as well as 1 µL of Droplet Break Color (Samplix ApS). After vortexing, the clear break solution phase could be removed, leaving the colored aqueous phase containing the enriched DNA.

Since the DNA amount after sorting was insufficient for direct sequencing, the captured transcripts were amplified using MDA, employing non-sequence-specific primers. Additionally, the Xdrop Sort was used to package single transcripts in single emulsion (SE) droplets to avoid bias toward smaller transcripts. For this purpose, 10 µL of enriched DNA from each sample was transferred to separate PCR-tubes. Furthermore, 1 pg of unenriched cDNA dissolved in 10 µL nuclease-free water was taken along as a positive control, as well as a nontemplate control. To each of the samples, 1 µL of SE MDA enzyme (Samplix ApS), 4 µL of SE MDA mix (5×) (Samplix ApS), and 5 µL of nuclease-free water were added on ice. Subsequently, 20 µL of each sample was loaded onto the Xdrop SE85 Cartridge (Samplix ApS) according to the manufacturer's instructions, after which 75 µL of SE Droplet Oil (Samplix ApS) was administered to each inlet well of the cartridge. The cartridge was then placed into the Xdrop Sort and the SE droplet production program was run, after which droplets were collected in PCR-tubes. All but 1–2 mm of droplet oil was removed from the bottom of the PCR-tubes, which were then placed in a thermal cycler to be incubated at 30°C for 16 h, 65°C for 10 min, and 4°C until the release of DNA from the droplets. To break the SE droplets, 20 µL of Droplet Break Solution is added to each of the tubes, together with 1 µL of Droplet Break Color. After vortexing, the clear break solution phase could be removed, leaving the colored aqueous phase containing the enriched and amplified DNA.

Enrichment validation was performed using qPCR to compare enriched samples with the original, unenriched cDNA pool. The amplified 1 pg of unenriched cDNA and nontemplate control from the MDA served as positive and negative controls, respectively. qPCR primers, designed to target 100–125 bp regions close to the targeted regions of the dPCR primers, were used in technical replicates to negate bias from unintended target capture. Primer compositions and targeted exons are detailed in Supplemental Table S6. Each qPCR reaction contained 1 µL of template DNA, 10 µL of GoTaq 2x Master Mix (Promega Corporation), 1 µL of forward and reverse qPCR primers, and 7 µL of nuclease-free water. Reactions were run on a QuantStudio 3 Real-Time PCR System (Thermo Fisher Scientific) with the following program: 50°C for 1 sec, 95°C for 10 min, 40 cycles of 95°C for 15 sec, and 60°C for 30 sec, then 95°C for 15 sec, 60°C for 1 min, and finally 95°C for 15 sec.

PacBio library preparation and sequencing—adjusted workflow for Samplix-enriched transcripts

To remove the branched structures present in the amplified transcripts generated through the MDA process, an enzymatic digestion was performed. To do so, the enriched samples were diluted with nuclease-free water to a total volume of 25.5 µL. To each sample, 3 µL NEBuffer 2 (New England Biolabs) and 1.5 µL of T7 Endonuclease (New England Biolabs) were added before incubation at 37°C for 15 min in a thermal cycler. Afterward, 20 µL of TE buffer (pH8) was added to each of the samples and a custom bead cleanup was performed. To do so, 40 µL of AMPure PB beads were resuspended and pelleted on a magnet. After the supernatant had been completely removed the beads were washed with nuclease-free water twice and resuspended in an equal volume of the custom bead buffer consisting of the following components: 20 µL of 1 M Tris-HCL (Sigma-Aldrich), 4 µL 0.5 M EDTA pH8 (Thermo Fisher Scientific), 640 µL 5 M NaCl (Thermo Fisher Scientific), 550 µL 40% PEG8000 (Sigma-Aldrich), and 778 µL nuclease-free water. Of this custom bead suspension, 35 µL was added to each sample and incubated at room temperature for 20 min. Subsequently, the samples were pelleted on a magnet and the supernatant was removed. Samples were then washed with 70% ethanol and dried for 30 sec. Afterward, samples were removed from the magnetic rack and resuspended in nuclease-free water before being incubated for 1 min at 50°C and 5 min at room temperature. After spinning the samples down, the beads were pelleted on the magnetic rack until the eluate was clear and the purified DNA could be collected. A total of 500 ng of MDA-amplified and purified DNA was used to create multiplexed SMRTbell amplicon libraries, following the manufacturer's instructions for the SMRTbell prep kit 3.0 102-359-000 (Pacific Biosciences). If the sample had <500 ng of DNA, linearized plasmids not containing the region of interest were added to spike the samples. The samples were then prepared for sequencing using the Sequel IIe Binding Kit 3.2 (Pacific Biosciences), following the recommended protocols from SMRT Link 11.0.0.146107. Finally, 115 µL of the final mix was loaded per well, and long-read sequencing was performed using the Sequel IIe system (Pacific Biosciences). Following sequencing, the reads were processed following a cDNA analysis: subreads were demultiplexed using lima V.2.5.0, and combined to create a consensus sequence using CCS V.6.3.0. These reads were filtered RQ 0.99 to obtain HiFi reads, and HiFi reads were then mapped along the GRCh38 reference genome using pbmm2 V1.8.0 with the –preset Iso-Seq mode. Instead of performing an IsoQuant analysis, we visualized the data in IGV, because despite the enzymatic digestion following the MDA procedure, debranched reads can still contain multiple fragments of the captured transcript and are incompatible with algorithms such as IsoQuant.

Tissue collection for ONT sequencing

For the independent ONT validation data set, the rest material from human donors (n = 3) without any known or clinical evidence of retinal disease was collected from the tissue banks of either Ghent University Hospital or Antwerp University Hospital. This collection adhered to the ethical standards of the Declaration of Helsinki and received approval from the Ethics Committee of Ghent University Hospital (IRB approval B670201837286). Details about the donors can be found in Supplemental Table S7. The eyes were transported in CO₂ Independent Medium (Gibco) before dissection. To preserve RNA integrity and minimize the effects of autolysis, retinas were only harvested from eyes with a total postmortem interval of <20 h. Following a visual inspection to ensure no contamination from the RPE, the neural retinas were either processed immediately for total RNA isolation or snap-frozen and stored at −80°C for later use.

RNA isolation and ONT sequencing

Total RNA was extracted from the postmortem adult human neural retina samples using the RNeasy Mini kit (Qiagen), following the manufacturer's instructions. The extracted RNA then underwent DNase treatment (ArcticZymes, Tromsø, Norway) followed by poly(A) capture. Samples of poly(A) mRNA with sufficient quality (RNA Integrity Value, RIN > 8.0) were used for direct-cDNA library preparation using the SQK-DCS109 kit (ONT), with minor modifications to the supplier's protocol. Each prepared library was then loaded onto a FLO-PRO002 flow cell (ONT) and sequenced on an ONT PromethION device for 72 h. Details regarding the number of reads can be found in Supplemental Table S7.

ONT data analysis of selected genes

MinKNOW version 5.1.0 was used to produce FAST5 files, which were subsequently base-called with Guppy version 6.1.5. The reads were then aligned to the Human Reference Genome GRCh38/hg38 using minimap2 (Li 2021) version 2.24 with the -ax splice flags for spliced alignments. Alignment files were converted to BAM format, sorted, and indexed using SAMtools version 1.15 (Li et al. 2009) to produce the data set investigated here.

Relative expression of MYO7A and WHRN isoforms

Quantitative PCR analysis was conducted to determine the relative expression levels of the different MYO7A- and WHRN transcript isoforms identified in the PacBio Iso-Seq data sets across three human neural retina samples. The same RNA used for preparing the PacBio libraries was employed as the template, with 250 ng of RNA being used for cDNA synthesis using the SuperScript IV Reverse Transcriptase kit (Thermo Fisher Scientific 18090200). Quantitative PCR analysis was performed using GoTaq qPCR Master Mix (Promega), following the manufacturer's protocol. For MYO7A, transcript-specific primers were designed and validated to target the 5′ canonical start site, the 5′ alternative start site, and the 3′ end of the identified MYO7A isoforms, as well as the reference gene GUSB. Amplifications were carried out using the QuantStudio 3 Real-Time PCR system (Applied Biosystems, Waltham, MA, USA), with PCR reactions performed in triplicate. Relative gene expression levels compared to the reference gene GUSB were determined using the 2^−ΔCt method.

For WHRN, we designed and validated transcript-specific primers to target isoforms with intron 4 retention and those containing exon 7B. In the absence of a universally present WHRN transcript region, we included a primer pair targeting WHRN exons 8–9 to capture a segment present in nearly all isoforms, designated as “total” WHRN. Expression levels of WHRN transcripts were determined using the 2^−ΔCt method, and plotted relative to the expression of the exons 8–9 target representing “total WHRN transcripts” set at 1. The primers used are listed in Supplemental Table S8.

Analysis of protein isoforms

To predict the 2D protein domain architecture of encoded proteins in silico analyses were performed using the SMART online tool (Letunic et al. 2021). 3D protein structures were modeled using AlphaFold (Jumper et al. 2021) using the Google Colab notebook (v1.5.3) with standard settings. The 3D structures were visualized with YASARA (Krieger and Vriend 2014), and structural alignments were conducted using the MUSTANG algorithm (Konagurthu et al. 2006).