Study design

Nine patients who presented to the Charité outpatient clinic because of suspected post-COVID syndrome and two patients that were hospitalized at the Charité were included between June 2020 and November 2021 based on the following criteria:

(1)

Age >18 years

(2)

PCR-proven SARS-CoV-2 infection

(3)

Persistent muscular fatigue and post exertional malaise (PEM) first manifesting after infection with SARS-CoV-2 and lasting for at least 6 months

(4)

Exclusion of other causes explaining the symptoms listed under (3)

(5)

Approval for and absence of contraindication for vastus lateralis muscle biopsy

All patients signed informed consent before study inclusion and the study was approved by the Ethics Committee of the Charité—Universitätsmedizin Berlin (EA2/066/20 and EA2/163/17) in accordance with the 1964 Declaration of Helsinki and its later amendments.

Nine of the eleven patients were part of a larger prospective observational study [14]. These patients were seen at least once in the outpatient clinic, when a detailed clinical evaluation and neurological examination (muscle strength testing of major muscle groups, handgrip strength test, reflexes and sensory testing, 6-min-walk-test) was performed and serum samples were obtained. On at least one other occasion, study participants responded to online questionnaires (Bell and Chalder fatigue questionnaires, EQ-5D-5L) [87, 88] hosted in a secure REDCap database as previously described [14]. Proximal lower extremity MRI was performed in these patients on the same day or close to the day of the biopsy. Two other patients (PCS-7 and -10) were included out of an inpatient setting based on the above-mentioned inclusion criteria. After inclusion, one of them (PCS-7) was diagnosed with a rheumatoid arthritis and primary biliary cholangitis. Clinical records were consulted for age, sex, preexisting medical conditions, onset and nature of acute and chronic clinical symptoms, laboratory results, therapeutic measures, and complications. All included patients were contacted by telephone beginning of December 2022 for a final follow-up evaluation, assessing the current state of perceived symptoms.

Vastus lateralis muscle biopsy

Biopsies were taken according to standard procedures as previously described [89]. In short, after informed consent for the procedure had been granted, open biopsy of the vastus lateralis muscle was performed under local anesthesia (lidocaine 2%). After circumscribed incision of the skin with removal of a 2 × 5 mm section of the cutis, the muscle was carefully removed. A 15 × 15 mm part of the skeletal muscle, and a 5 × 10 mm part of the muscle fascia were acquired for histopathological assessment. After biopsy procedure at the Department of Neurosurgery at Charité – Universitätsmedizin Berlin, muscle tissue, and fascia were processed immediately at the Department of Neuropathology, Charité – Universitätsmedizin Berlin.

Control cohorts

Cryopreserved skeletal muscle specimens obtained from the vastus lateralis muscle were selected based on the following inclusion criteria:

(1)

Age between 18 and 65 years

(2)

Vastus lateralis muscle biopsy prior to December 2019

(3)

Absence of known inflammatory disease, cancer or mitochondriopathy

(4)

Absence of increased creatinine kinase levels, pathological EMG, corticosteroid or other immunosuppressive therapy

(5)

For the healthy disease control (HDC) cohort: absence of any histopathological or immunohistochemical abnormality in the excised muscle tissue

(6)

For the type-2b atrophy control cohort (2BA): presence of a selective atrophy of type-2b-fibers but other than that absence of any histopathological or immunohistochemical abnormality in the excised muscle tissue

These biopsies had been performed for routine diagnostic reasons, and patients had consented to further processing of their samples for scientific purposes. We chose the term “healthy disease control (HDC)” as these patients were clinically diseased (symptoms reported in Additional file 1: Table S1) justifying a muscle biopsy, but did not show any histological abnormalities. Due to the high prevalence of women in our PCS cohort, we preferentially selected samples of women fulfilling the above-mentioned criteria. However, due to restricted numbers of available samples from women that fulfilled all our inclusion criteria, the HDC group was composed of five men and three women. The type-2b-fiber atrophy control cohort on the other hand consisted of women only.

The mean age of the HDC group was 42.6 years (SD 10.9; median 45 years), the one of the 2BA group was 44.6 years (SD 14.0; median 44.5 years).

Magnetic resonance imaging (MRI)

MRI scans were acquired on a 3 Tesla scanner (MAGNETOM PRISMA®, Siemens, Erlangen, Germany). The subjects were examined in the supine position and feet first using a 28-channel sensitivity encoding torso array coil placed anteriorly. Total scan duration was approximately 35 min and included qualitative imaging by axial and coronal T2-weighted turbo spin echo (TSE), axial T1-weighted TSE, axial T1-weighted 3D volumetric interpolated breath-hold examination (VIBE) with Dixon fat suppression and reconstruction of in-/opposed phase, water- and fat-based images, axial 2D Spin-Echo (SE) T2 mapping, details are contained in Additional file 1: Fig. S1. The field of view (FOV) was at mid-thigh level and anatomic T1w/T2w imaging, T2 maps, and DTI were acquired using the same FOV and geometry. The mean diffusivity (MD), T1 and T2 relaxation times, and muscle quantitative fat fraction (MFF [%]) were evaluated using Visage Imaging Client (Software Release v7.1, Visage Imaging). Manual seeding of regions-of-interest (ROI) at mid-thigh level avoiding areas of fatty infiltration or vascular structures in the biceps femoris (BF), semitendinosus (ST), semimembranosus (SM), and vastus lateralis (VL) muscle was conducted by a radiologist with more than 8 years’ experience in musculoskeletal MRI, blinded to the patients’ clinical data. The MFF was calculated using axial 3D gradient echo-modified two-point Dixon-based MRI with a chemical shift-encoded reconstruction of the water and fat signal as i) SIFAT / (SIFAT + SIWATER) × 100 and reported as mean value of all pixels within the ROI.

One patient (PCS-3) interrupted the examination before the MD and T2 relaxation times were acquired.

Virological analysis

Unfixed, cryopreserved muscle samples were used for detection and quantification of SARS-CoV-2 RNA by quantitative reverse transcription–polymerase chain reaction (RT-qPCR). Only samples with at least two positive results were considered positive. Oligonucleotides targeting the leader transcriptional regulatory sequence and a region within the single-guide RNA encoding the SARS-CoV-2 E gene were used to detect single-guide RNA as described previously [90, 91].

Anti-SARS-CoV-2 IgG enzyme-linked immunosorbent assays with S1 and NCP domain substrate were performed in available serum samples according to the manufacturer’s instructions (Euroimmun AG®, Lübeck, Germany). In addition, electrochemiluminescence immunoassay (ECLIA) antigen tests were performed to detect S- and N-antigens according the manufacturer’s instructions (Elecsys®, Roche, Basel, Switzerland).

Autoimmunity assays

Antinuclear antibody assays (HEp2 –IFT), myositis-specific autoantibodies (anti–nuclear matrix protein-2 [anti-NXP2], anti–transcriptional intermediary factor 1γ [anti-TIF1γ], anti–melanoma differentiation-associated gene 5 [anti-MDA5], anti–signal recognition particle [anti-SRP], anti-Mi2, anti-isoleucyl-transfer RNA [tRNA] synthetase [anti-OJ], anti–glycyl-tRNA synthetase [anti-EJ], anti–threonyl-tRNA synthetase [anti-PL7], anti–alanyl-tRNA synthetase [anti-PL12], anti-histidyl-tRNA synthetase [anti-Jo1], and anti–small ubiquitin-like modifier-1 activating enzyme [anti-SAE]), and myositis-associated autoantibodies (anti-Ku, anti-PM75, anti-PM100, and anti-Ro52) were performed in available serum samples according to the manufacturer’s instructions (ANA-Mosaik 1A EUROPattern and EUROLINE ANA-Profil 3 (IgG), EUROIMMUN Medizinische Labordiagnostika AG, Lübeck, Germany).

Histology & immunohistochemistry

Unfixed biopsy specimens were snap-frozen in a container with isopentane in liquid nitrogen and stored at − 80 °C until further workup. Stainings on cryopreserved samples were performed on 7 μm thick cryomicrotome sections. Routine histological and enzymological staining (hematoxylin–eosin, Gömöri trichrome, periodic acid–Schiff, ATPases) were carried out according to standard procedures. Immunohistochemical staining was performed on a Benchmark XT autostainer (Ventana Medical Systems), as described previously [92].

For quantification of immune cell populations (CD68-, CD169-, CD206-, CD45- and CD8-positive cells, respectively) and semiquantitative scoring (degree of MHC class I & II upregulation and of type-2b-fiber atrophy), ten random fields of vision were examined independently by two experienced morphologists (W.S. and T.A.) at × 400 magnification with an Olympus BX50 microscope (Ocular WH10X-H/22). Semiquantitative scoring of type 2b fiber atrophy: 0 = no atrophic 2b fibers; 1 = atrophy of < 15% of 2b fibers; 2 = atrophy of 15–60% of 2b fibers; 3 = atrophy of > 60% of 2b fibers + many fibers < 20 µm in diameter. Positively stained immune cells were counted manually in 10 high-power fields of vision. Positive staining results with MHC class I and MHC class II were defined as a clear upregulation at the sarcolemma with capillaries and arterioles serving as internal positive controls. The following antibodies were used: MHC class I (DAKO; clone W6/32, 1:100), MHC class II (DAKO; M0775, 1:100) CD45 (DAKO; clone UCHL1, 1:100), CD68 (DAKO; clone EBM11, 1:100), CD8 (DAKO; clone C8/144B, 1:100), NKp46 (R&D Systems; clone MAB1850, 1:100), Siglec-1/CD169 (Novus Biologicals; clone HSn 7D2, 1:200), C5b-9 (DAKO/M777; clone aE11 1:100), CD206 (Abnova clone 5C11; 1:50).

Transmission electron microscopy and capillary morphometry

Electron microscopy was performed as described previously [91]. In short, after fixation in 2.5% glutaraldehyde in 0.1 M sodium cacodylate buffer muscle samples were incubated with 1% osmium tetroxide in 0.05 M sodium cacodylate and embedded in Renlam resin after dehydration by a graded acetone series. Semithin sections of 500 nm were cut with an ultramicrotome (Ultracut E, Reichert-Jung) and a Histo Jumbo diamond knife (Diatome) and stained with toluidine blue at 80 °C. Ultrathin sections of 70 nm were cut using the same ultramicrotome and an Ultra 35° diamond knife (Diatome) and stained with uranyl citrate. Standard transmission EM was performed using a Zeiss 906 microscope in conjunction with a 2 k CCD camera (TRS).

For each patient, between 20 and 30 randomly selected capillaries from at least two different Renlam resin blocks were photographed at a final magnification of 7000 × . Blurry images with unclear basement membrane borders as well as images of abnormally large microvessels or capillaries with very high pericyte coverage were excluded from the analysis. Capillary basement membrane (CBM) thickness was measured at six distinct sites per capillary with Image J v1.53c (NIH, Bethesda, MD, United States), omitting areas close to profiles of pericyte processes where the CBM is irregular and usually thicker at these sites. For the evaluation of compartmental organization of capillaries, micrographs showing capillary profiles with an aspect ratio (ratio of the smallest to largest diameter) of more than 1.2 were considered too obliquely sectioned and were excluded from morphometric evaluation, as previously recommended [93].

Tablet-based image analysis (TBIA) was performed for capillary morphometry, as previously described [94]. On electron micrographs of capillary profiles, lines were drawn with a digital pen using ImageJ around the capillary lumen (blood:EC transition), along the abluminal EC surface (EC:BM transition), along the basement membrane (BM):endomysium transition, and around the PC surface to obtain values for areas and circumferences. Absolute arithmetic values for the lumen radius, and the EC and BM thicknesses were calculated using formulae previously described [95].

For subjective scoring of capillary alterations, images were interpreted in a blinded fashion, side-by-side by two neuropathologists (W.S. and H.H.G.) and independently by one physiologist (O.B.) with long-standing experience in ultrastructural analysis of skeletal muscle capillaries. The following scores were used per capillary: 0 =  < 4 pericyte processes; no reduplication; no ensheathment; maximum 2–3 endothelial cells; vesicles and mitochondria within normal range; 1 = single prominent pericytes; single or no ensheathment; slightly increased number and size of vesicles and mitochondria; 2 = prominent pericytes with irregular structure; ensheathment, prominent vesicles and mitochondria; 3 = same as 2 but more pronounced; 4 = loss; necrosis of endothelial cells; capillary remnants; cellular debris.

Morphometry of light micrographs for the evaluation of capillarity

Semithin sections were prepared as described above and per patient at least 10 images at a magnification of 400 × were taken with a Keyence BZ-X810. After anonymization of the obtained images, each of these light micrographs was overlaid with a digital counting grid consisting of 10 × 10 test lines, as previously reported [96]. For calculation of the capillary-to-fiber (C/F) ratio, the number of capillary profiles and that of muscle fibers were counted within the counting grid, taking into account the so-called forbidden line rule. For an estimate of the mean cross-sectional fiber area (MCSFA), the number of test points falling on fiber profiles was divided by the total number of test points and multiplied by the total area of the counting grid (= total area of muscle fiber profiles) divided by the number of muscle fiber profiles that fell within the confines of the grid. The capillary density was determined by dividing the number of capillary profiles that fell within the boundaries of the grid by the total area of muscle fiber profiles.

Muscle sample preparation for proteomics

Muscle specimens were lysed in 200 µl of 50 mM TEAB (pH 8.5) buffer, 5% SDS, and complete ULTRA protease inhibitor (Roche) using the Bioruptor® (Diagenode) for 10 min (30 s on, 30 s off, 10 cycles) at 4 °C. To ensure complete lysis we conducted an additional sonication step using an ultra-sonic probe (30 s, 1 s/1 s, amplitude 40%) followed by centrifugation at 4 °C and 20,000 g for 15 min. Protein concentration of the supernatant was determined by BCA assay according to the manufacturer’s protocol. Disulfide bonds were reduced by addition of 10 mM TCEP at 37 °C for 30 min, and free sulfhydryl bonds were alkylated with 15 mM IAA at room temperature (RT) in the dark for 30 min. 100 µg protein of each sample was used for proteolysis using the S-Trap protocol (Protifi) and using a protein to trypsin ratio of 20:1. The incubation time for trypsin was changed to 2 h at 47 °C.

All proteolytic digests were checked for complete digestion after desalting by using monolithic column separation (PepSwift monolithic PS-DVB PL-CAP200-PM, Dionex) on an inert Ultimate 3000 HPLC (Dionex, Germering, Germany) by direct injection of 1 μg sample. A binary gradient (solvent A: 0.1% TFA, solvent B: 0.08% TFA, 84% ACN) ranging from 5 to 12% B in 5 min and then from 12 to 50% B in 15 min at a flow rate of 2.2 μL/min and at 60 °C, was applied. UV traces were acquired at 214 nm [97].

Proteomic analysis of muscle samples

An UltiMate 3000 RSLC nano UHPLC coupled to a QExactive HF mass spectrometer was used for the analysis of all samples, and the total amount of peptides used was always 1 µg. Samples were first transferred to a 75 µm × 2 cm, 100 Å, C18 precolumn at a flow rate of 10 µl/min for 20 min followed by separation on the 75 µm × 50 cm, 100 Å, C18 main column with a flow rate of 250 nl/min and a linear gradient composed of solution A (99.9% water, 0.1% formic acid) and solution B (84% acetonitrile, 15.9% water, 0.1% formic acid), with a pure gradient length of 120 min (3–45% solution B). The gradient was applied as follows: 3% B for 20 min, 3–35% for 120 min, followed by 3 wash steps, each reaching 95% buffer B for 3 min. After the last wash step, the instrument was equilibrated for 20 min. MS data were acquired in data independent acquisition (DIA) mode using an in-house generated spectral library for the corresponding tissue or body fluid. Each sample was mixed with an appropriate amount of iRT standard (Biognosys). Full MS scans were acquired from 300 to 1100 m/z at a resolution of 60,000 (Orbitrap) using the polysiloxane ion at 445.12002 m/z as the lock mass. The automatic gain control (AGC) was set to 3E6 and the maximum injection time was set to 20 ms. The full MS scans were followed by 23 DIA windows, each covering a range of 28 m/z with an overlap of 1 m/z, starting at 400 m/z, acquired at a resolution of 30,000 (Orbitrap) with an AGC of 3E6 and an nCE of 27 (CID).

For the analysis of samples acquired by nano-LC–MS/MS in DIA mode, the data were entered into Spectronaut software (Biognosys) and analyzed using a library-based search. The library utilized were the in-house created spectral libraries, depending on the tissue or body fluid type. The search and extraction settings were kept as default (BGS Factory settings). Human proteome data from UniProt (www.uniprot.org) with 20,374 entries were selected as the proteome background. For reliable label-free quantification, only proteins with ≥ 2 unique peptides were considered for further analysis. The average normalized abundances were determined using Spectronaut for each protein and used to determine the ratio between the patients' samples and the corresponding controls.

Quantitative reverse transcription PCR (qRT-PCR)

Total RNA was extracted from muscle specimens using a trizol-chloroform method as described previously [92] and cDNA was synthesized using the High-Capacity cDNA Archive Kit (Applied Biosystems, Foster City, CA, USA). For qPCRs, 10 ng of cDNA was used for subsequent analysis, using an Applied Biosystems™ QuantStudio™ 6 Flex Real-Time PCR System (ThermoFischer, Waltham, MA, USA) with the following running conditions: 95 °C for 20 s, 95 °C for 1 s, 60 °C for 20 s, for 45 cycles (values above 40 cycles were defined as not expressed). All targeted transcripts were run as triplicates. For each of these runs, the reference gene GAPDH has been included as an internal control to normalize the relative expression of the targeted transcripts. The qPCR assay identification numbers, TaqMan® Gene Exp Assay from Life Technologies/ ThermoFisher are listed as follows: GAPDH Hs02786624_g1, COL4A1 Hs00266237_m1, COL4A2 Hs05006309_m1, NID1 Hs00915875_m1, NID2 Hs00201233_m1, HSPG2 Hs01078536_m1. The ΔCT of HDCs was subtracted from the ΔCT of COVID patients muscles to determine the differences (ΔΔCT) and fold change (2^-ΔΔCT) of gene expression.

Bulk RNA-sequencing

RNA was isolated using Trizol (ThermoFisher) and the DirectZol kit (Zymo) according to manufacturer’s instruction. Poly-A RNA sequencing libraries were then prepared using the NEBNext Ultra II Directional RNA Library Prep Kit (NEB), and sequenced to a depth of 20–30 million paired-end reads on a NovaSeq 6000 device (Illumina) with 2 × 109 sequencing length. Since the samples were sequenced together with samples containing high amounts of SARS-CoV-2 RNA on the same Novaseq sequencing lane, there were some reads aligning to the SARS-CoV-2 genome due to index hopping on the flowcell, i.e. misassignement of samples with high SARS-CoV-2 amount to the muscle samples. Due to this effect, quantification of SARS-CoV-2 RNA from sequencing is less reliable than from RT-qPCR.

RNA sequencing analysis

Due to data protection restrictions, raw sequencing reads cannot be made readily available, however a read count table is provided as online Additional file 1: data. Raw sequencing reads were aligned to version hg19 of the human genome using hisat2 [98] with standard parameters. Read counts were quantified based on the hg19 RefSeq annotation using quasR [99]. Differential expression values were calculated using edgeR [100]. Further analysis was done and plots were generated using the R packages PCAtools [PCAtools: Blighe, K., Lun, A. PCAtools: everything principal components analysis. https://github.com/kevinblighe/PCAtools. R package version 2.10.0, https://github.com/kevinblighe/PCAtools], ComplexHeatmap [ComplexHeatmap: Gu Z (2022). “Complex Heatmap Visualization.” iMeta. https://doi.org/10.1002/imt2.43.], as well as packages from the Tidyverse [https://doi.org/10.21105/joss.01686]. All code is available on https://github.com/landthalerlab.

Pathway enrichement analysis of proteomic and transcriptomic datasets

For pathway analysis of the proteomic datasets, significantly regulated proteins (p-value ≤ 0.05) with either positive or negative regulation (calculated fold change of ≥ 1.2 or ≤ 0.8) were included. For pathway analysis of the transcriptomic datasets, significantly regulated genes (p-value ≤ 0.01) with either positive or negative regulation (calculated log fold change of ≥ 0.8 or ≤ -0.8) were included using Metascape (https://metascape.org) [101].

Statistical analysis

Statistical analyses were done and graphs created with GraphPad Prism 9 (GraphPad Software) and R Software (R Core Team (2022) R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2022. Available from: https://www.R-project.org/).

Normality testing was performed (D’Agostino & Pearson, Anderson–Darling, Shapiro–Wilk and Kolmogorov–Smirnov tests) and a Gaussian distribution was confirmed for the following parameters: immune cell quantification (CD68, CD169, CD206, CD45, CD8), capillary-to-fiber-ratio, MCFSA, capillary density, CBM thickness. One-way Anova with Tukey’s multiple comparison test was used to compare differences in these parameters between the three cohorts. Data are presented as counts, percentages or means (SDs). Values were considered significant at p ≤ 0.05. The significance of the abundance changes in the proteomics and transcriptomic datasets were calculated with DESeq [102]. Abundance changes were considered significant at p ≤ 0.05 (proteomics) and at p ≤ 0.001 (transcriptomics). Heat map and dot plot visualizations were created with R Software and tables with Excel 2016 (Microsoft).