Transcriptomic analysis of human ALS skeletal muscle reveals a disease-specific pattern of dysregulated circRNAs

Human sample collection

Spinal cord and frontal cortex tissue samples were harvested post-mortem from patients enrolled in an ALS tissue collection program (directed by Peter H. King) approved by the University of Alabama at Birmingham (UAB) Institutional Review Board (IRB-100908007 and IRB-091222037). Muscle biopsy samples were obtained from an archive of tissue remnants in the UAB Neuromuscular Division. Neuropathy and myopathy samples were chosen based on histological evidence of denervation or myopathic changes, as determined by a neuromuscular pathologist. All samples were stored at −80°C until use. All patients had a definite diagnosis of ALS based on revised El Escorial criteria.

Mouse sample collection

All animal protocols were approved by the Institutional Animal Care and Use Committee at UAB in compliance with the National Research Council Guide for the Care and Use of Laboratory Animals. B6.Cg-Tg (SOD1^G93A) 1 Gur/J mice (the Jackson Laboratory, ME) were mated with C57BL/6J to generate hemizygous SOD1^G93A offspring. Age-matched wild-type (WT) littermates were used as controls. Clinical progression was evaluated as described previously [18]. Tissue samples were collected from SOD1^G93A and WT littermate controls at post-natal days 60, 125, and 150 as previously described [18].

Motor neuron differentiation

Human induced pluripotent stem cell (iPSC) lines were maintained at 70% confluency in supplemented media purchased from STEMCELL Technologies [37]. To begin iPSCs differentiation, cells were dissociated into single cells with ReLeSR and then seeded into Matrigel-coated 6-well plates. The cells were then treated with Iscove’s modified Dulbecco’s medium (IMDM) supplemented with B27–vitamin A (2%), N2 (1%), SB431542, CHIR99021, and LDN193189 for 6 days to force them towards a neuroectodermal fate. Fifty percent of the culture medium was replenished every other day. To push cells toward neural progenitors, the cells were gently dissociated into single cells by accutase treatment on day 6 and grown in neural differentiation media supplemented with caudalizing factors all-trans retinoic acid (Sigma R2625) and smoothened agonist (SAG; Cayman 11914) for the next 6 days. Starting at day 12, caudalized neural progenitor cells (NPCs) were maintained in IMDM containing neurogenic small molecules brain-derived neurotrophic factor (BDNF; PeproTech) and glial cell-derived neurotrophic factor (GDNF; Peprotech) for the next 20 days which differentiated NPCs into motor neurons. At day 32, cells were utilized in downstream analyses.

RNA isolation

For muscle tissues, samples were homogenized in TRIzol reagent (Thermo Fisher Scientific). Chloroform was added to the lysate, and then thoroughly mixed by vortexing. After centrifugation of muscle sample lysates, the aqueous phase was transferred to an RNAspin Mini Column (GE Healthcare) to isolate total RNA following the manufacturer’s protocol. The cultured iPSC-derived motor neurons were harvested and lysed in TRIzol followed by mixing with chloroform, then a Qiagen RNeasy mini kit was used for the following steps. Briefly, after centrifugation, the aqueous phase containing RNA was transferred to a new microfuge tube, and the RNA was precipitated with ethanol. Precipitated RNA was dissolved, passed through the RNeasy column, and treated with DNase to remove genomic DNA contamination. The RNA was eluted in RNase/DNase-free water and quantified using a Nanodrop instrument.

circRNA enrichment, library preparation, and sequencing

RNA quality and quantity were assessed using the Agilent RNA 6000 nano kit in 2100 Bioanalyzer System (Agilent). High-quality RNA (500 ng) was used to prepare the sequencing library with an Illumina TruSeq Stranded Total RNA Library prep kit following the manufacturer’s protocol. Briefly, after rRNA depletion and cDNA generation, the cDNAs were subjected to 3’ end adenylation, adapter ligation, and purification with AMPure beads (Beckman Coulter). The products were size-selected with SPRIselect beads (Beckman Coulter), and the selected cDNAs were enriched by PCR and repurified with SPRIselect beads to generate the final libraries. The quality and quantity of sequencing libraries were checked using Agilent DNA 1000 Screen Tape on the Agilent TapeStation. Single-read sequencing was performed for 124 cycles with an Illumina NovaSeq 6000 sequencer.

CircRNA was enriched as previously described [38] with minor modifications. Briefly, 2 μg of RNA were treated with 20 U of RNase R (RNR07250, Epicentre) for 1 h at 37°C, after which RNA was isolated using Direct-zol RNA kit (Zymo Research) following the manufacturer’s instructions; the isolated RNA was polyadenylated using a poly(A) tailing kit (Invitrogen) and then selected using Oligo(dT)₂₅ Dynabeads (Invitrogen) following the manufacturer’s instructions. From the supernatant, which contained the enriched circRNA pool, RNA was isolated using Direct-zol RNA kit (Zymo Research) following the manufacturer’s instructions; libraries were generated using the entire material, omitting the rRNA depletion steps, and quantified as described above. Paired-end sequencing was performed for 105 cycles with an Illumina HiSeq 2500 sequencer.

Reverse transcription (RT) followed by quantitative (q)PCR analysis

For reverse transcription (RT) followed by quantitative PCR (qPCR) analysis, 500 ng of total RNA was used. For qPCR analysis, 0.1 μl cDNA was employed with 250 nM of primers (Supplementary File 1) and KAPA SYBR^® FAST qPCR Kits (KAPA Biosystems) as described [39]. Divergent primers spanning the circRNA junctions of interest were designed using CircInteractome [40]. RT-qPCR analysis was carried out on a QuantStudio 5 Real-Time PCR System (Thermo Fisher Scientific) with a cycle setup of 3 min at 95°C and 40 cycles of 5 s at 95°C, and 20 s at 60°C. For the mRNAs examined, primers were designed in exons different from those predicted to constitute the circRNA body. Relative RNA levels were calculated after normalizing to the levels of RPS9 mRNA (human) or Rps9 mRNA (mouse) using the 2^−ΔΔCt method. The levels of housekeeping TBP mRNA (human) or Tbp mRNA (mouse) were measured as additional controls. The same procedure was followed in both human and mouse samples.

Bioinformatic analysis of sequencing data

BCL files were de-multiplexed and converted to FASTQ files using bcl2fastq program (v2.20.0.422). FASTQ files were trimmed for adapter sequences using Cutadapt version v1.18 and aligned to human genome hg19 Ensembl v82 using STAR software v2.4.0j [41]; featureCounts (v1.6.4) [42] were used to create gene counts from the samples for linear RNA analysis. The chimeric junction file obtained from STAR software was parsed for fusion junctions and analyzed using CIRCexplorer v1.1.10 [43] to obtain the circularizing junction counts for circRNA analysis as well as for circRNA annotation. The RNA sequencing data are deposited in GSE215424.

For both the linear and the circRNA transcripts, read counts were normalized and differential abundance between healthy and ALS individuals was assessed using the DESeq2 package version 1.30.0 pipeline [44] in R (version 4.0.3). Briefly, the dispersion of samples was explored by using a principal component analysis (PCA) scattergram of regularized logarithm (rlog) transformed counts. Statistical testing was performed using the Wald test. Linear RNAs were defined as differentially regulated with an absolute log2 fold change > 1 and Benjamini–Hochberg adjusted p-value < 0.05. Given the low counts and highly variable abundance of circRNAs, we focused on identifying only abundant circRNAs by requiring one or more counts in at least 40% of the ALS muscle samples. Using this stringent cutoff, 250 circRNA transcripts were selected and subjected to further analysis based on their fold change; significance was established based on an unadjusted p-value < 0.05.

Immunostaining

At day 32, iPSCs differentiated in 96-well plates were fixed in 4% (v/v) paraformaldehyde (PFA) and then rinsed with phosphate-buffered saline (PBS). Fixed cultures were blocked in 5% (v/v) donkey serum with 0.03% (v/v) Triton X-100 and 2% fetal bovine serum (FBS). Blocked cell preparations were incubated with primary antibodies anti-SMI32 (BioLegend 801701, dilution 1:250) raised in mouse and anti-Isl1 (Abcam ab109517, dilution 1:500) raised in rabbit. After overnight incubation with the primary antibodies at 4°C, fixed cells were rinsed with PBS and incubated in species-specific secondary antibodies conjugated with Alexa Fluor 488 or Alexa Fluor 594 (dilution 1:1000) at room temperature for 2 h. Nuclei were counterstained with DAPI followed by rinsing with PBS. Cells were visualized using a fluorescent microscope (EVOS M5000) at 10× magnification.

Statistical analysis

Quantitative data are represented as the means ± SD of the number of samples indicated in each case; statistical significance was established using unpaired Student’s t-test in GraphPad Prism (9.0). A p-value of < 0.05 was considered statistically significant and was indicated in the figures as ^*p < 0.05, ^**p < 0.01, ^***p < 0.001. Graphs were generated using GraphPad Prism (9.0).

Transcriptomic analysis of normal and ALS skeletal muscle

To begin to identify differentially expressed RNAs in ALS, we used muscle biopsy samples from ALS patients followed in our clinic with a definite diagnosis of ALS (Table 1). Histologically normal muscle biopsies were used as a control. For later validation of our findings, we also included additional normal and ALS muscle biopsies, as well as muscle samples from neuropathy and myopathy disease controls (Table 1 and Supplementary Table 1). We isolated total RNA from muscle biopsies (5 healthy and 5 ALS individuals) and performed total RNA sequencing (RNA-seq) analysis (Figure 1A). The RNA sequencing data are deposited in GSE215424. Principal component analysis (PCA) using the linear transcripts revealed a distinct separation of the ALS cohort compared to the normal cohort (Figure 1B).

To gain a better insight into the ALS cohort, we assessed the differentially expressed linear coding and noncoding transcripts (Figure 1C and Supplementary File 2). We found that several mRNAs previously associated with ALS were significantly elevated (adjusted p-value < 0.05), including collagen type XIX, alpha 1 (COL19A1) mRNA, previously associated with higher mortality and fast progression [45], and SMAD9 (also known as SMAD8) mRNA, encoding a muscle biomarker that is broadly increased in ALS muscle and tracks disease progression [18]. Other transcripts, including inter-alpha-trypsin inhibitor heavy chain family member 6 (ITIH6) mRNA and myosin light chain kinase (MYLK2) mRNA, which have roles in muscular function, were significantly (adjusted p-value < 0.05) downregulated [46, 47]. Taken together, our ALS cohort shared common transcriptomic features with previous reports, supporting its suitability for further assessment of circRNAs.

Differentially expressed circRNAs in normal and ALS skeletal muscle

To identify circRNAs in the muscle biopsies, we enriched our RNA samples by first digesting them with RNase R and then polyadenylating any remaining 3’-OH ends in order to eliminate any residual linear RNAs, as explained ([38] and Materials and Methods). RNA-seq analysis of the circRNA-enriched samples initially identified 4,656 unique junctions; however, after adopting a stringent cutoff criterion requiring at least one junction count in at least 40% of the donors (Material and Methods), we narrowed the results to only 250 circRNAs (Supplementary Figure 1 and Supplementary File 3). Among them, only 4 and 6 circRNAs were significantly (unadjusted p-value < 0.05) up- or downregulated, respectively (Figure 1D). Considering the notorious difficulties in detecting circRNAs by RNA-seq analysis and the intrinsic heterogeneity of ALS cohorts, we decided to test the differential expression of at least 50 circRNAs, 25 increased (Supplementary Figure 2A) and 25 decreased (Supplementary Figure 2B), that showed more pronounced fold changes in the RNA-seq analysis using a broader sampling of circRNAs in ALS muscle biopsies. Using the CircInteractome [40] and circBase [48] databases, we identified the circRNA aliases, as well as the counterpart mRNAs arising from the shared linear precursor transcripts, and performed RT-qPCR analysis to assess both forms, the circRNA and the corresponding mRNA. Among the tested circRNAs, 8 were significantly upregulated whereas 10 were significantly downregulated in the ALS cohort (Figure 2); the additionally tested circRNAs did not significantly change or changed in a direction that was opposite to that suggested by the RNA-seq analysis (Supplementary Figure 3). Notably, among the differentially expressed circRNAs based on RT-qPCR analysis, 3 and 4 were predicted to be significantly upregulated or downregulated, respectively, based on the RNA-seq analysis.

For the upregulated group, 7 out of 8 circRNAs, hsa_circ_0000119, hsa_circ_0000567, hsa_circ_0007778, hsa_circ_0000099, hsa_circ_0005171, hsa_circ_0000231, and hsa_circ_0001451, were altered without significant changes in the linear counterpart transcripts (MAN1A2, SETD3, SLC25A16, AMY2B, TMEM38B, ARHGAP12, and FBXW7 mRNAs, respectively) suggesting that these circRNAs may have a role in ALS, independent of the associated mRNAs; only hsa_circ_0007099 and its counterpart (ABHD2 mRNA) showed joint increases in abundance (Figure 2A). Interestingly, hsa_circ_0000567 was previously reported as a blood biomarker in ALS [49]. To evaluate if these circRNAs were specifically elevated in ALS, we assessed their expression levels in muscle biopsies from control myopathy and neuropathy samples (Table 1 and Figure 2B). Notably, the levels of hsa_circ_0000119, hsa_circ_0000567, hsa_circ_0007778, hsa_circ_0005171, and hsa_circ_0007099 did not change significantly (or changed in the opposite direction) in these disease controls, supporting the possibility that their elevation was specific to ALS. Other circRNAs and/or mRNA counterparts changed in myopathy or neuropathy, suggesting that perhaps their altered levels in ALS were linked to common pathways of neuromuscular dysfunction.

For the downregulated circRNAs, 6 out of 10 circRNAs in ALS muscle biopsies, hsa_circ_0047886, hsa_circ_0141401, hsa_circ_0117010, hsa_circ_0000033, hsa_circ_0000247, and hsa_circ_0009027, were selectively reduced while the linear counterparts (ZCCHC2, ALPK2, NPHP1, CEP85, MCU, and ZNF362 mRNAs, respectively) were not. In addition, for hsa_circ_0056856, hsa_circ_0001948, and hsa_circ_0000944, both the circRNA and the linear counterparts (ITGB6, CD99L2, and CCDC9 mRNAs, respectively) decreased, whereas hsa_circ_0006633 also decreased while the linear counterpart (FGGY mRNA) increased (Figure 2C). Interestingly, many of the downregulated circRNAs were also lower in myopathy and neuropathy disease controls, including hsa_circ_0056856, hsa_circ_0047886, hsa_circ_0141401, hsa_circ_0117010, hsa_circ_0000944, and hsa_circ_0009027 (Figure 2D). Taken together, we have identified 7 selectively increased and 6 decreased circRNAs in ALS muscle biopsies, without changes in their counterpart mRNAs; neuropathy and myopathy controls had distinct patterns of expression, sometimes mirroring the differential expression of the circRNAs, sometimes not.

ALS muscle circRNA patterns are distinct from ALS central nervous system (CNS) circRNA patterns

Given the extensive connectivity between the peripheral neuromuscular system and the CNS, we assessed the levels of the validated circRNAs (Figure 2) in post-mortem spinal cord (SC) and frontal cortex (FC) tissues from normal and ALS individuals (Table 2 and Figure 3). From the group of upregulated muscle circRNAs, none was increased in either SC or FC (Figure 3A and Supplementary Figure 4A). In fact, several circRNAs, including hsa_circ_0000119, hsa_circ_0000567, hsa_circ_0005171, and hsa_circ_0000231, showed modest but significant reductions in SC; hsa_circ_0000119, hsa_circ_0000567, and hsa_circ_0000231 also decreased in FC. One circRNA, hsa_circ_0000099, was unchanged in SC but decreased in FC; none of the counterpart mRNAs in the entire group was altered in SC or FC (Figure 3A and Supplementary Figure 4A).

For the downregulated muscle circRNAs, on the other hand, 8 circRNAs were similarly reduced in SC, including hsa_circ_0047886, hsa_circ_0141401, hsa_circ_0001948, hsa_circ_0000033, hsa_circ_0009027, hsa_circ_0006633, hsa_circ_0000944, and hsa_circ_0009027, whereas hsa_circ_0056856 and hsa_circ_0117010 were unchanged (Figure 3B and Supplementary Figure 4B). Surprisingly, in FC, hsa_circ_0056856, hsa_circ_0141401, and hsa_circ_0000944, as well as their corresponding mRNAs, were significantly increased (Figure 3B and Supplementary Figure 4B). In summary, there was an interesting disparity in the circRNAs elevated in ALS muscle (Figure 2A), as their expression was lower or unchanged in ALS SC and/or FC (Figure 3A); in contrast, circRNAs reduced in ALS muscle (Figure 2C) were generally unchanged or reduced in ALS CNS (Figure 3B).

Shared circRNA profiles between ALS muscle and ALS iPSC-derived motor neurons

We extended our analysis to motor neurons derived from iPSCs of ALS patients carrying hexanucleotide repeat expansions in C9ORF72 (C9-ALS) (Table 3 and Figure 4). iPSC-derived motor neurons from healthy subjects served as control cells. We performed immunocytochemistry on day-32 neurons and found that they expressed two markers of mature motor neurons, ISL1 and SMI32, confirming their motor neuron phenotype (Figure 4A). Among those circRNAs upregulated in ALS skeletal muscle, C9-ALS motor neurons showed similar increases in hsa_circ_0000119, hsa_circ_0000567, hsa_circ_0007778, hsa_circ_0000099, and hsa_circ_0005171 relative to control motor neurons, while the mRNA counterparts for these circRNAs were not altered; the greatest fold change was seen for hsa_circ_0000119 and hsa_circ_0000099 (Figure 4B). The other circRNAs that were upregulated in skeletal muscle in the ALS cohort did not change in the ALS motor neurons. Regarding the circRNAs that were reduced in ALS muscle, hsa_circ_0056856, hsa_circ_0117010, and hsa_circ_0006633 showed a similar decrease in motor neurons; interestingly, the mRNA counterparts of hsa_circ_0056856 and hsa_circ_0141401 (ITGB3 and ALPK2 mRNAs, respectively) were modestly but significantly elevated in the ALS motor neurons, in contrast to what was observed in skeletal muscle (Figure 4C). Most other circRNAs reduced in ALS muscle did not significantly decrease in ALS iPSCs (Figure 4C). Taken together, the pattern of increased circRNAs in iPSC-derived motor neurons of ALS patients shared patterns similar to those observed in ALS skeletal muscle.

Muscle circRNAs can track disease progression in the SOD1^G93A mouse

The SOD1^G93A transgenic mouse recapitulates many features of ALS including progressive weakness, motor neuron loss, and muscle denervation compared to age-matched wild-type (WT) littermate controls, thus providing an opportunity to assess temporal patterns of biomarkers with disease progression [50–52]. We assessed the expression of altered circRNA candidates in mouse skeletal muscle at three ages which reflect different disease stages based on rotarod performance and weight measurements [18]: 60 days old [presymptomatic (PS)], 125 days old [early symptomatic (ES)], and 150 days old [late symptomatic (LS)] (Figure 5A, skeletal muscle). We first identified mouse circRNA orthologs using the flanking exons of the human circRNA candidates to find the conserved exons in the respective mouse linear ortholog. We then designed divergent primers spanning the predicted junction in mice. A summary of the aliases of the host genes and exons predicted to constitute the circRNA bodies in each species, as well as the predicted sequence overlaps are shown in Supplementary Figure 5.

At the PS stage, four mouse circRNA orthologs, circ_0007099, circ_0005171, circ_0056856, and circ_0000033, showed differential expression in gastrocnemius muscle compared to age-matched WT controls (Figure 5B). Circ_0007099 and circ_0005171 were upregulated in the SOD1^G93A mouse, similar to the human circRNA ortholog, and gradually increased with disease progression through the LS stage, whereas the respective mRNA counterparts were unchanged. Likewise, circ_0056856 and circ_0000033 were downregulated in the SOD1^G93A mouse, again similar to what was seen for the human orthologs, but the suppression was relatively constant through the LS stage. With the exception of the PS stage for circ_0056856, the corresponding mRNAs were unchanged. In the ES stage, circ_0000119 and circ_0047886 started to show differential expression in the same direction as in human ALS muscle, with the former showing upregulation and the latter showing downregulation; notably, the expression changes did not progress in the LS stage (Figure 5C). Similar to the human ortholog, the linear counterpart of circ_0000119 (Man1a2 mRNA) was unchanged compared to WT, whereas the linear counterpart of circ_0047886 (Zcchc2 mRNA) was significantly lower (Figure 5C). Circ_0000567, previously reported as a biomarker in ALS [49], was the only one differentially expressed at the LS stage without a change in the mRNA counterpart (Figure 5D); the circRNA levels increased in ways similar to those observed in human ALS muscle. The other candidates, predicted to either increase or decrease in human ALS muscle, showed no differential expression at any stage (Supplementary Figure 6).

We then assessed whether the spinal cord tissues reflected the differential expression of the mouse muscle circRNAs that showed disease stage-associated patterns (Figure 5E, spinal cord). Four circRNAs, circ_0000119, circ_0000567, circ_00047886, and circ_0000033, showed differential expression at the ES stage. Two of the upregulated muscle circRNAs, circ_0000119 and circ_0000567, were lower in the spinal cord, similar to their human counterparts at the ES stage; the reduced levels persisted but did not further decline through the LS stage. The levels of the linear counterpart of circ_0000119 (Man1a2 mRNA) did not change, while the levels of the linear counterpart of circ_0000567 (Setd3 mRNA) decreased in parallel through the LS stage (Figure 5F). The downregulated muscle circRNAs circ_0047886 and circ_0000033 decreased in the spinal cord similar to what was seen for the human orthologs. The linear counterpart for circ_0000033 (Cep85 mRNA) also decreased, albeit modestly, in contrast to what we saw in human ALS spinal cord (Figure 5F). Lastly, circ_0005171 was downregulated in the LS stage (Figure 5G), whereas circ_0007099 and circ_0056856 did not change at any disease stage (Figure 5H). Among other prominent muscle circRNA candidates, circ_0141041 showed differential expression at the PS stage and progressively decreased through the LS stage (Supplementary Figure 7A), whereas circ_0000231, circ_0000099, circ_0001948, and circ_0000944 decreased at the LS stage in spinal cord (Supplementary Figure 7B). Other candidate muscle circRNAs had variable patterns or no significant change (Supplementary Figure 7C).

A summary of the changes for circRNAs based on RT-qPCR analysis in human and mouse tissues is shown in Table 4. Several circRNAs showed consistent changes in human and mouse skeletal muscle (5 upregulated and 4 downregulated circRNAs). Several circRNAs also showed consistent changes (all decreased) in human and mouse spinal cord tissues, including hsa_circ_0000119, hsa_circ_0000567, hsa_circ_0005171, hsa_circ_0047886, hsa_circ_0001948, and hsa_circ_0000033. We note an interesting difference in direction of abundance, with those circRNAs elevated in muscle (top half of Table 4) showing reduced abundance in spinal cord.