ミオシン重鎖を超えたヒト骨格筋線維の異質性

nature.com にアクセスいただきありがとうございます。ご利用のブラウザバージョンでは、CSS のサポートが制限されています。最適なエクスペリエンスを得るには、最新バージョンのブラウザをご利用いただくか、Internet Explorer の互換モードを無効にしてください。また、継続的なサポートを確保するため、このサイトではスタイルシートと JavaScript は使用しません。
骨格筋は、主に筋原線維からなる不均一な組織であり、ヒトでは通常、1 つの「遅筋」（タイプ 1）と 2 つの「速筋」（タイプ 2A および 2X）の 3 種類に分類されます。しかし、従来の筋原線維タイプ間およびタイプ内の不均一性については、まだ十分に解明されていません。私たちは、ヒト外側広筋由来の 1,050 個と 1,038 個の筋原線維に対して、それぞれトランスクリプトームおよびプロテオームの手法を適用しました。プロテオーム研究には男性が含まれ、トランスクリプトーム研究には男性 10 名と女性 2 名が含まれていました。ミオシン重鎖アイソフォームに加えて、代謝タンパク質、リボソームタンパク質、および細胞接合タンパク質が、多次元的な筋原線維間の変動の原因であることがわかりました。さらに、遅筋と速筋のクラスターが特定されたにもかかわらず、私たちのデータは、タイプ 2X 線維が他の速筋線維と表現型的に区別できないことを示唆しています。さらに、ミオシン重鎖に基づく分類は、ネマリンミオパチーにおける筋線維の表現型を記述するには不十分である。全体として、我々のデータは、ミオシン重鎖アイソフォームを超えた変異源を伴う、多次元的な筋線維の不均一性を示唆している。
細胞の異質性はすべての生物系に内在する特徴であり、細胞が組織や細胞の異なるニーズを満たすために特殊化することを可能にしてきました。1 骨格筋線維の異質性に関する従来の見解は、運動ニューロンが運動単位内の線維タイプを定義し、線維タイプ（ヒトではタイプ1、タイプ2A、タイプ2X）はミオシン重鎖（MYH）アイソフォームの特性によって決定されるというものでした。2 これは当初、pH ATPaseの不安定性3,4 に基づき、後にMYHの分子発現5 に基づきました。しかし、複数のMYHをさまざまな割合で共発現する「混合」線維が特定され、その後受け入れられるようになったため、骨格筋線維は、明確な線維タイプではなく、連続体として見られることが多くなりました。6 それにもかかわらず、この分野では、筋線維分類の主な分類基準としてMYHに大きく依存しており、これは、MYH発現プロファイルと線維タイプの範囲が異なる初期のげっ歯類研究の限界と大きなバイアスの影響を受けている可能性があります。 2 人間の骨格筋はそれぞれ多様な繊維タイプを示すため、状況はさらに複雑になっています。7 外側広筋は、中間的な（したがって代表的な）MYH 発現プロファイルを持つ混合筋です。7 さらに、サンプル採取が容易なため、人間の筋肉の中で最も研究されている筋肉となっています。
したがって、強力な「オミクス」ツールを用いて骨格筋線維の多様性を偏りなく研究することは重要であると同時に、骨格筋線維の多核性という性質もあって困難を伴う。しかしながら、近年の様々な技術進歩により、トランスクリプトミクス8,9 およびプロテオミクス10 技術の感度は飛躍的に向上し、単一線維レベルの解像度で骨格筋を解析することが可能となった。その結果、単一線維の多様性、そして萎縮刺激や老化に対する反応の特性評価において大きな進歩が遂げられた11,12,13,14,15,16,17,18。重要なのは、これらの技術進歩が臨床応用され、疾患に関連する調節異常をより詳細かつ正確に評価できるようになることである。たとえば、最も一般的な遺伝性筋疾患の 1 つであるネマリンミオパチー (MIM 605355 および MIM 161800) の病態生理は複雑でわかりにくいものです。19,20 そのため、骨格筋線維の調節不全をより適切に特徴付けることで、この疾患に対する理解が大きく進む可能性があります。
我々は、ヒト生検標本から手作業で単離した単一骨格筋線維のトランスクリプトームおよびプロテオーム解析法を開発し、数千の線維に適用することで、ヒト骨格筋線維の細胞内異質性を調査することができました。この研究を通して、筋線維のトランスクリプトームおよびプロテオームによる表現型解析の有効性を実証し、代謝タンパク質、リボソームタンパク質、および細胞接合タンパク質が線維間変異の重要な原因であることを明らかにしました。さらに、このプロテオーム解析ワークフローを用いて、単一骨格筋線維における線虫性ミオパチーの臨床的意義を明らかにし、MYHに基づく線維の種類に関わらず、非酸化性線維への協調的な移行を明らかにしました。
ヒト骨格筋線維の不均一性を調べるために、我々は単一の骨格筋線維のトランスクリプトームおよびプロテオーム解析を可能にする 2 つのワークフローを開発した (図 1A および補足図 1A)。サンプルの保管と RNA およびタンパク質の完全性の保持から各アプローチのスループットの最適化まで、いくつかの方法論的ステップを開発し、最適化した。トランスクリプトーム解析では、逆転写の最初のステップでサンプル固有の分子バーコードを挿入することでこれを実現し、96 の線維をプールして効率的な下流処理を可能にした。従来の単一細胞アプローチと比較してより深いシーケンス (線維あたり±100 万リード) により、トランスクリプトームデータがさらに強化された。21 プロテオミクスでは、短いクロマトグラフィーグラジエント (21 分) と timsTOF 質量分析計での DIA-PASEF データ取得を組み合わせて、高スループットを維持しながらプロテオームの深さを最適化した。 22,23 健康な骨格筋線維の不均一性を調べるために、14 人の健康な成人ドナーから採取した 1,050 個の線維のトランスクリプトームと、5 人の健康な成人ドナーから採取した 1,038 個の線維のプロテオームを特徴付けました (補足表 1)。この論文では、これらのデータセットをそれぞれ 1,000 線維トランスクリプトームとプロテオームと呼びます。私たちのアプローチでは、1,000 線維のトランスクリプトーム解析とプロテオーム解析で合計 27,237 個のトランスクリプトと 2,983 個のタンパク質が検出されました (図 1A、補足データセット 1 ～ 2)。検出された遺伝子が 1,000 個を超え、線維あたり 50% の有効値を持つトランスクリプトームとプロテオームのデータセットをフィルタリングした後、トランスクリプトームとプロテオームのそれぞれ 925 個と 974 個の線維に対してバイオインフォマティクス解析を実行しました。フィルタリング後、線維あたり平均4257 ± 1557個の遺伝子と2015 ± 234個のタンパク質（平均値±標準偏差）が検出され、個人間変動は限定的でした（補足図1B～C、補足データセット3～4）。しかし、被験者内変動は参加者間でより顕著であり、これはおそらく、異なる長さや断面積を持つ線維間のRNA/タンパク質収量の違いによるものと考えられます。ほとんどのタンパク質（2000個以上）について、変動係数は20%未満でした（補足図1D）。どちらの方法も、筋収縮に重要な高発現シグネチャー（例：ACTA1、MYH2、MYH7、TNNT1、TNNT3）を持つ転写産物とタンパク質の幅広いダイナミックレンジを捕捉することができました（補足図1E～F）。特定された特徴のほとんどは、トランスクリプトミクス データセットとプロテオーム データセット間で共通しており (補足図 1G)、これらの特徴の平均 UMI/LFQ 強度はかなりよく相関していました (r = 0.52) (補足図 1H)。
トランスクリプトミクスとプロテオミクスのワークフロー（BioRender.comで作成）。BD MYH7、MYH2、MYH1のダイナミックレンジ曲線と、繊維タイプ割り当てのための計算された閾値。E、F トランスクリプトミクスとプロテオミクスのデータセットにおける繊維全体にわたるMYH発現の分布。G、H MYHベースの繊維タイプ別に色分けされたトランスクリプトミクスとプロテオミクスのUniform Diversity approximation and Projection (UMAP)プロット。I、J トランスクリプトミクスとプロテオミクスのデータセットにおけるMYH7、MYH2、MYH1発現を示す特徴プロット。
我々は当初、オミクスデータセットにおけるMYH発現の高感度およびダイナミックレンジを活用する最適化されたアプローチを用いて、各線維にMYHベースの線維タイプを割り当てることに着手した。これまでの研究では、異なるMYHの発現の固定割合に基づいて、任意の閾値を用いて線維を純粋なタイプ1、タイプ2A、タイプ2X、または混合として分類していた11,14,24。我々は、線維を分類するために使用したMYHによって各線維の発現をランク付けするという異なるアプローチを用いた。MYH7、MYH2、およびMYH1は、それぞれタイプ1、タイプ2A、およびタイプ2X線維に対応する。次に、得られた各曲線の下側の変曲点を数学的に算出し、それを閾値として使用して、各MYHについて線維を陽性（閾値以上）または陰性（閾値以下）として割り当てた（図1B～D）。これらのデータは、MYH7（図1B）とMYH2（図1C）が、タンパク質レベルと比較してRNAレベルでより明確なオン/オフ発現プロファイルを有することを示しています。実際、タンパク質レベルでは、MYH7を発現しない線維はごくわずかで、MYH2を100%発現する線維はありませんでした。次に、事前に決定した発現閾値を使用して、各データセットのすべての線維にMYHベースの線維タイプを割り当てました。たとえば、MYH7+/MYH2-/MYH1-線維はタイプ1に割り当てられ、MYH7-/MYH2+/MYH1+線維は混合タイプ2A/2Xに割り当てられました（詳細は補足表2を参照）。すべての線維をプールすると、RNAレベル（図1E）とタンパク質レベル（図1F）の両方でMYHベースの線維タイプの分布が著しく類似していることがわかりましたが、MYHベースの線維タイプの相対的な構成は予想どおり個人間で異なっていました（補足図2A）。ほとんどの線維は純粋なタイプ1（34～35%）またはタイプ2A（36～38%）に分類されましたが、タイプ2A/2X混合線維も相当数検出されました（16～19%）。顕著な違いは、純粋なタイプ2X線維はRNAレベルでのみ検出され、タンパク質レベルでは検出されなかったことです。これは、MYHの急速な発現が少なくとも部分的に転写後に制御されていることを示唆しています。
抗体を用いたドットブロッティングを用いて、プロテオミクスに基づくMYH線維タイピング法を検証しました。両法とも、純粋なタイプ1線維とタイプ2A線維の同定において100%の一致を達成しました（補足図2B参照）。しかし、プロテオミクスに基づくアプローチは、より感度が高く、混合線維の同定と各線維における各MYH遺伝子の割合の定量においてより効率的でした。これらのデータは、客観的かつ高感度なオミクスに基づくアプローチを用いて骨格筋線維のタイプを特徴づけることの有効性を実証しています。
次に、トランスクリプトミクスとプロテオミクスから得られた情報を組み合わせて、完全なトランスクリプトームまたはプロテオームに基づいて筋線維を客観的に分類しました。均一多様体近似および投影法（UMAP）を使用して次元を6つの主成分に削減し（補足図3A〜B）、トランスクリプトーム（図1G）とプロテオーム（図1H）の筋線維の変動を視覚化できました。注目すべきは、トランスクリプトミクスデータセットでもプロテオミクスデータセットでも筋線維が参加者（補足図3C〜D）または試験日（補足図3E）ごとにグループ化されていなかったことです。これは、骨格筋線維の被験者内変動が被験者間変動よりも大きいことを示唆しています。UMAPプロットでは、「速い」筋線維と「遅い」筋線維を表す2つの明確なクラスターが現れました（図1G〜H）。 MYH7+（遅筋）筋線維はUMAP1の正極にクラスターを形成し、MYH2+およびMYH1+（速筋）筋線維はUMAP1の負極にクラスターを形成しました（図1I～J）。しかし、MYHの発現に基づく速筋線維のタイプ（すなわち、タイプ2A、タイプ2X、または2A/2X混合）の区別は認められませんでした。これは、MYH1（図1I～J）やACTN3、MYLK2（補足図4A～B）などの他の古典的な2X筋線維マーカーの発現は、トランスクリプトーム全体またはプロテオーム全体を考慮すると、異なる筋線維タイプを区別しないことを示しています。さらに、MYH2およびMYH7と比較して、MYH1と正の相関を示す転写産物またはタンパク質は少なかった（補足図4C～H）。これは、MYH1の存在量が筋線維のトランスクリプトーム／プロテオームを完全に反映していないことを示唆している。3つのMYHアイソフォームの混合発現をUMAPレベルで評価した場合も同様の結論に達した（補足図4I～J）。したがって、2X線維はMYHの定量のみに基づいて転写産物レベルで同定できるものの、トランスクリプトームまたはプロテオーム全体を考慮すると、MYH1+線維は他の速筋線維と区別できない。
MYH以外の遅線維の異質性に関する初期調査として、TPM3、TNNT1、MYL3、およびATP2A22という4つの確立された遅線維タイプ特異的タンパク質を評価しました。遅線維サブタイプは、トランスクリプトミクス（補足図5A）およびプロテオミクス（補足図5B）の両方で、MYH7と高いものの、完全ではないピアソン相関を示しました。トランスクリプトミクス（補足図5C）およびプロテオミクス（補足図5D）において、それぞれ約25％および33％の遅線維がすべての遺伝子/タンパク質サブタイプによって純粋な遅線維として分類されませんでした。したがって、複数の遺伝子/タンパク質サブタイプに基づく遅線維分類は、線維タイプ特異的であることがわかっているタンパク質の場合でも、さらなる複雑さをもたらします。これは、単一の遺伝子/タンパク質ファミリーのアイソフォームに基づく線維分類が、骨格筋線維の真の異質性を適切に反映していない可能性があることを示唆しています。
ヒト骨格筋線維の表現型の多様性をオミックスモデル全体の規模でさらに調査するために、主成分分析（PCA）を用いてデータの不偏次元削減を行った（図2A）。UMAPプロットと同様に、参加者も試験日もPCAレベルでは線維のクラスタリングに影響を与えなかった（補足図6A〜C）。どちらのデータセットでも、MYHベースの線維タイプはPC2によって説明され、PC2は遅筋タイプ1線維のクラスターと、速筋タイプ2A、タイプ2X、および混合2A/2X線維を含む2番目のクラスターを示した（図2A）。どちらのデータセットでも、これら2つのクラスターは少数の混合タイプ1/2A線維によって接続されていた。予想通り、主要なPCドライバーの過剰表現分析により、PC2は収縮および代謝シグネチャによって駆動されていることが確認された（図2Bおよび補足図6D〜E、補足データセット5〜6）。全体的に、MYH ベースの繊維タイプは、高速クラスター内のトランスクリプトーム全体に分布するいわゆる 2X 繊維を除いて、PC2 に沿った連続的な変化を説明するのに十分であることがわかりました。
A. MYHに基づいて繊維タイプ別に色分けされたトランスクリプトームおよびプロテオームデータセットの主成分分析（PCA）プロット。B. PC2およびPC1におけるトランスクリプトおよびタンパク質ドライバーのエンリッチメント分析。統計解析は、clusterProfilerパッケージとBenjamini-Hochberg補正p値を用いて実施しました。C、D. トランスクリプトーム内の細胞間接着遺伝子オントロジー（GO）用語とプロテオーム内のコスタメアGO用語に基づいて色分けされたPCAプロット。矢印は、トランスクリプトおよびタンパク質ドライバーとその方向を示しています。E、F. 臨床的に関連する特徴の均一多様体近似および投影（UMAP）特徴プロット。これは、遅い/速い繊維タイプとは無関係に発現勾配を示しています。G、H. トランスクリプトームおよびプロテオームにおけるPC2およびPC1ドライバー間の相関関係。
予想外にも、MYHに基づく筋線維型は、2番目に高い変動性（PC2）のみを説明した。これは、MYHに基づく筋線維型（PC1）とは無関係な他の生物学的因子が、骨格筋線維の不均一性を制御する上で重要な役割を果たしていることを示唆している。PC1における主要な因子の過剰発現解析により、PC1の変動性は主に細胞間接着とトランスクリプトーム中のリボソーム含量、およびプロテオーム中のコスタメアとリボソームタンパク質によって決定されていることが明らかになった（図2Bおよび補足図6D～E、補足データセット7）。骨格筋において、コスタメアはZ板と筋鞘を連結し、力の伝達とシグナル伝達に関与している。細胞間接着（トランスクリプトーム、図 2C）およびコスタメア（プロテオーム、図 2D）の特徴を使用した 25 注釈付き PCA プロットでは、PC1 の強い左シフトが明らかになり、これらの特徴が特定の繊維に豊富に含まれていることが示されました。
UMAPレベルでの筋線維のクラスタリングをより詳細に調べた結果、ほとんどの特徴は、筋線維サブクラスター特異的ではなく、筋線維タイプに依存しないMYHベースの発現勾配を示していることが明らかになりました。この連続性は、CHCHD10（神経筋疾患）、SLIT3（筋萎縮）、CTDNEP1（筋疾患）など、病態に関連するいくつかの遺伝子で観察されました（図2E）。この連続性は、神経疾患（UGDH）、インスリンシグナル伝達（PHIP）、転写（HIST1H2AB）に関連するタンパク質を含むプロテオーム全体にわたっても観察されました（図2F）。これらのデータを総合すると、異なる筋線維間で、線維タイプに依存しない遅筋/速筋の異質性に連続性があることが示唆されます。
興味深いことに、PC2のドライバー遺伝子は良好なトランスクリプトーム-プロテオーム相関（r = 0.663）を示した（図2G）。これは、遅筋線維と速筋線維のタイプ、特に骨格筋線維の収縮特性と代謝特性が転写制御されていることを示唆している。しかし、PC1のドライバー遺伝子はトランスクリプトーム-プロテオーム相関を示さなかった（r = -0.027）（図2H）。これは、遅筋線維/速筋線維のタイプに関連しない変異は、主に転写後に制御されていることを示唆している。PC1の変異は主にリボソーム遺伝子オントロジー用語によって説明されており、リボソームはタンパク質翻訳に積極的に関与し、影響を与えることで細胞内で極めて重要かつ特殊な役割を果たしていることから31、次に、この予期せぬリボソームの不均一性について調査することにした。
まず、GOCC用語「細胞質リボソーム」内のタンパク質の相対的存在量に応じて、プロテオミクス主成分分析プロットを色分けしました（図3A）。この用語はPC1の正側に濃縮されているため、勾配は小さいですが、リボソームタンパク質はPC1の両方向への分配を促進します（図3A）。PC1の負側に濃縮されているリボソームタンパク質にはRPL18、RPS18、RPS13（図3B）が含まれ、RPL31、RPL35、RPL38（図3C）はPC1の正側の主な推進力でした。興味深いことに、RPL38とRPS13は、他の組織と比較して骨格筋で高く発現していました（補足図7A）。PC1のこれらの特徴的なリボソームシグネチャーはトランスクリプトームでは観察されず（補足図7B）、転写後調節を示しています。
A. プロテオーム全体にわたる細胞質リボソーム遺伝子オントロジー（GO）用語に基づいて色分けされた主成分分析（PCA）プロット。矢印はPCAプロットにおけるタンパク質を介した変動の方向を示す。線の長さは、特定のタンパク質の主成分スコアに対応する。B、C. RPS13およびRPL38のPCA特徴プロット。D. 細胞質リボソームタンパク質の教師なし階層的クラスタリング解析。E. 骨格筋線維における異なる存在量を持つリボソームタンパク質を強調表示した80Sリボソーム（PDB: 4V6X）の構造モデル。F. mRNA出口チャネル付近に局在する、異なる化学量論を持つリボソームタンパク質。
リボソームの不均一性と特殊化の概念は以前にも提唱されており、異なるリボソームサブポピュレーション（リボソームの不均一性）の存在が、特定の mRNA 転写プールの選択的翻訳 34 （リボソームの特殊化）を介して、異なる組織 32 および細胞 33 でのタンパク質翻訳に直接影響を与える可能性があるとされています。骨格筋線維で共発現しているリボソームタンパク質のサブポピュレーションを特定するために、プロテオーム中のリボソームタンパク質の教師なし階層的クラスタリング解析を実施しました（図 3D、補足データセット 8）。予想どおり、リボソームタンパク質は MYH に基づく線維タイプごとにクラスター化されませんでした。ただし、リボソームタンパク質の 3 つの異なるクラスターを特定しました。最初のクラスター（ribosomal_cluster_1）は RPL38 と共制御されているため、PC1 プロファイルが陽性の線維で発現が増加しています。 2番目のクラスター（ribosomal_cluster_2）はRPS13と共発現しており、PC1プロファイルが陰性の線維で増加しています。3番目のクラスター（ribosomal_cluster_3）は骨格筋線維において協調的な差次的発現を示さず、「コア」骨格筋リボソームタンパク質と考えることができます。リボソームクラスター1と2はどちらも、代替翻訳を制御すること（例：RPL10A、RPL38、RPS19、RPS25）や、発達に機能的に影響を与えること（例：RPL10A、RPL38）が以前に示されたリボソームタンパク質を含んでいます。34,35,36,37,38 PCAの結果と一致して、線維全体にわたるこれらのリボソームタンパク質の不均一な表現も連続性を示しました（補足図7C）。
リボソーム内の異種リボソームタンパク質の位置を可視化するために、ヒト80Sリボソームの構造モデル（タンパク質データバンク：4V6X）を使用しました（図3E）。異なるリボソームクラスターに属するリボソームタンパク質を単離したところ、それらの位置は密接に整列していませんでした。これは、私たちのアプローチではリボソームの特定の領域/画分の濃縮に失敗したことを示唆しています。しかし興味深いことに、クラスター2の大サブユニットタンパク質の割合は、クラスター1および3よりも低かったです（補足図7D）。骨格筋線維において化学量論が変化したタンパク質は、主にリボソーム表面に局在していることがわかりました（図3E）。これは、異なるmRNA集団内の内部リボソーム進入部位（IRES）エレメントと相互作用して選択的翻訳を調整する能力と一致しています。 40, 41 さらに、骨格筋線維中の化学量論が変化した多くのタンパク質は、mRNA 出口トンネル (図 3F) などの機能領域の近くに位置しており、これは特定のペプチドの翻訳伸長と停止を選択的に制御します。42 まとめると、私たちのデータは、骨格筋リボソームタンパク質の化学量論が不均一性を示し、その結果、骨格筋線維間に違いが生じることを示唆しています。
次に、速筋線維と遅筋線維のシグネチャーを特定し、その転写制御機構を解明することに着手しました。2つのデータセット（図1G～Hおよび4A～B）においてUMAPによって定義された速筋線維と遅筋線維のクラスターを比較したところ、トランスクリプトーム解析とプロテオーム解析により、それぞれ1,366個と804個の異なる量の遺伝子が同定されました（図4A～B、補足データセット9～12）。サルコメア（トロポミオシン、トロポニンなど）、興奮収縮連関（SERCAアイソフォーム）、エネルギー代謝（ALDOA、CKBなど）に関連するシグネチャーに、予想通りの違いが見られました。さらに、タンパク質のユビキチン化を制御する転写産物とタンパク質は、速筋線維と遅筋線維で異なる発現を示しました（USP54、SH3RF2、USP28、USP48など）（図4A～B）。さらに、微生物タンパク質遺伝子RP11-451G4.2（DWORF）は、子羊の筋線維タイプ間で発現が異なること43 が示されており、心筋のSERCA活性を増強することが44 示されていますが、骨格筋遅線維で有意にアップレギュレーションされていました（図4A）。同様に、個々の線維レベルでは、代謝関連乳酸脱水素酵素アイソフォーム（LDHAおよびLDHB、図4Cおよび補足図8A）45,46 などの既知のシグネチャー、およびこれまで知られていなかった線維タイプ特異的シグネチャー（IRX3、USP54、USP28、およびDPYSL3など）で有意な差が観察されました（図4C）。トランスクリプトームデータセットとプロテオームデータセット間では、発現が異なる特徴の有意な重複が見られ（補足図8B）、サルコメア特徴のより顕著な発現差によって主に駆動される倍率変化の相関も見られました（補足図8C）。特に、一部のシグネチャ（USP28、USP48、GOLGA4、AKAP13など）は、プロテオームレベルでのみ強力な転写後調節を示し、遅筋/速筋線維タイプ特異的な発現プロファイルを持っていました（補足図8C）。
図1G～Hの均一多様体近似および投影（UMAP）プロットによって識別された低速クラスターと高速クラスターを比較したAおよびBのボルケーノプロット。色のついた点はFDR < 0.05で有意に異なる転写産物またはタンパク質を表し、濃い点はlog変化> 1で有意に異なる転写産物またはタンパク質を表します。双方向統計分析は、DESeq2 Wald検定とBenjamini–Hochberg調整済みp値（トランスクリプトミクス）またはLimma線形モデル法と経験的ベイズ分析に続いて多重比較のためのBenjamini–Hochberg調整（プロテオミクス）を使用して実行されました。C 低速繊維と高速繊維の間で選択された差次的発現遺伝子またはタンパク質のシグネチャプロット。D 有意に差次的発現した転写産物およびタンパク質のエンリッチメント解析。重複する値は両方のデータセットでエンリッチされ、トランスクリプトーム値はトランスクリプトームでのみエンリッチされ、プロテオーム値はプロテオームでのみエンリッチされます。統計解析は、clusterProfilerパッケージを用いて、Benjamini-Hochberg補正p値を用いて実施しました。E. SCENICが導出した調節因子特異性スコアと繊維タイプ間のmRNA発現差に基づいて、SCENICによって同定された繊維タイプ特異的転写因子。F. 遅繊維と速繊維間で発現差のある特定の転写因子のプロファイリング。
次に、発現レベルが異なる遺伝子およびタンパク質の過剰発現解析を実施しました（図4D、補足データセット13）。2つのデータセット間で異なる特徴に関するパスウェイエンリッチメントを実施した結果、脂肪酸β酸化およびケトン体代謝プロセス（遅線維）、ミオフィラメント／筋収縮（それぞれ速線維および遅線維）、炭水化物分解プロセス（速線維）といった、予想通りの差異が明らかになりました。速線維では、セリン／スレオニンタンパク質ホスファターゼ活性も上昇しており、これはグリコーゲン代謝を制御することが知られている調節性ホスファターゼサブユニット（PPP3CB、PPP1R3D、およびPPP1R3A）などの特徴によって引き起こされました（47）（補足図8D～E）。速線維に多く含まれる他の経路としては、プロテオーム（補足図8F）中のプロセシング（P-）体（YTHDF3、TRIM21、LSM2）（転写後制御に関与する可能性がある（48）、およびトランスクリプトーム（補足図8G）中の転写因子活性（SREBF1、RXRG、RORA）が挙げられる。遅線維には、酸化還元酵素活性（BDH1、DCXR、TXN2）（補足図8H）、アミド結合活性（CPTP、PFDN2、CRYAB）（補足図8I）、細胞外マトリックス活性（CTSD、ADAMTSL4、LAMC1）（補足図8J）、および受容体-リガンド活性（FNDC5、SPX、NENF）（補足図8K）が多く含まれる。
遅筋/速筋線維タイプ特性の基礎となる転写調節についてのさらなる知見を得るために、SCENIC49（補足データセット 14）を用いて転写因子エンリッチメント解析を行った。多くの転写因子が速筋線維と遅筋線維の間で有意にエンリッチメントしていた（図 4E）。これには、以前に速筋線維の発達に関連付けられている MAFA50 などの転写因子や、これまで筋線維タイプ特異的遺伝子プログラムとは関連付けられていなかったいくつかの転写因子が含まれていた。これらのうち、PITX1、EGR1、および MYF6 は速筋線維で最もエンリッチメントされていた転写因子であった（図 4E）。対照的に、ZSCAN30 および EPAS1（別名 HIF2A）は遅筋線維で最もエンリッチメントされていた転写因子であった（図 4E）。これに一致して、MAFA は速筋線維に対応する UMAP 領域で高レベルで発現していたのに対し、EPAS1 は逆の発現パターンを示した（図 4F）。
既知のタンパク質コード遺伝子に加えて、ヒトの発達と疾患の制御に関与している可能性のある非コードRNAバイオタイプが多数存在します。51, 52 トランスクリプトームデータセットでは、いくつかの非コードRNAが繊維タイプ特異性を示します（図5Aおよび補足データセット15）。これには、遅筋に非常に特異的で、ミトコンドリアミオパチー患者の筋肉で減少していると報告されているLINC01405が含まれます。53 対照的に、lnc-ERCC5-5遺伝子（https://lncipedia.org/db/transcript/lnc-ERCC5-5:2）54に対応するRP11-255P5.3は、速筋タイプ特異性を示します。 LINC01405 (https://tinyurl.com/x5k9wj3h) と RP11-255P5.3 (https://tinyurl.com/29jmzder) はどちらも骨格筋特異性を示し（補足図 9A–B）、1 Mb のゲノム近傍には既知の収縮遺伝子が存在しないことから、これらの遺伝子は隣接する収縮遺伝子の制御ではなく、線維タイプの制御に特化した役割を果たしていることが示唆されます。LINC01405 と RP11-255P5.3 のそれぞれにおける遅筋線維／速筋線維タイプ特異的な発現プロファイルは、RNAscope を用いて確認されました（図 5B–C）。
A. 非コードRNA転写産物は、遅筋線維と速筋線維において顕著に制御されている。B. LINC01405とRP11-255P5.3の遅筋線維と速筋線維の種類特異性を示す代表的なRNAscope画像。スケールバー = 50 μm。C. RNAscopeによって測定された筋線維の種類特異的な非コードRNA発現の定量（独立した個人から3つの生検試料を取得し、各個人内の速筋線維と遅筋線維を比較）。統計解析は、両側スチューデントt検定を用いて実施した。箱ひげ図は中央値と第1四分位値、第3四分位値を示し、ヒゲは最小値と最大値を示している。D. 微生物タンパク質のde novo同定ワークフロー（BioRender.comで作成）。 E. 微生物タンパク質LINC01405_ORF408:17441:17358は、骨格筋遅筋線維に特異的に発現している（独立した参加者から採取した5つの生検検体で、各参加者の速筋線維と遅筋線維を比較）。統計解析は、リム線型モデル法と経験的ベイズ法を組み合わせ、さらにベンジャミニ・ホッホバーグ法を用いてp値調整を行った多重比較を行った。箱ひげ図は中央値、第1四分位値、第3四分位値を示し、ヒゲは最大値／最小値を示している。
最近の研究では、多くの推定上の非コード転写産物が転写された微生物タンパク質をコードしており、そのいくつかは筋肉機能を調節していることが示されています。44, 55 繊維タイプ特異性の可能性のある微生物タンパク質を同定するために、1000繊維トランスクリプトームデータセットで見つかった非コード転写産物 (n = 305) の配列を含むカスタム FASTA ファイルを使用して、1000繊維プロテオームデータセットを検索しました (図 5D)。22 の異なる転写産物から 197 の微生物タンパク質を同定し、そのうち 71 は骨格筋の遅筋と速筋の間で異なる調節を受けていました (補足図 9C および補足データセット 16)。LINC01405 については、3 つの微生物タンパク質産物が同定され、そのうちの 1 つはその転写産物と同様の遅筋特異性を示しました (図 5E および補足図 9D)。
我々は、個々の筋線維の大規模なプロテオーム特性評価のための包括的なワークフローを開発し、健常状態における線維の不均一性の調節因子を特定しました。このワークフローを適用して、ネマリンミオパチーが骨格筋線維の不均一性にどのように影響するかを理解しました。ネマリンミオパチーは、筋力低下を引き起こす遺伝性の筋疾患で、罹患した小児では、呼吸困難、脊柱側弯症、四肢運動制限などのさまざまな合併症を呈します。19,20 典型的には、ネマリンミオパチーでは、アクチンアルファ1（ACTA1）などの遺伝子の病原性変異によって遅筋線維の筋線維構成が優位になりますが、この影響は不均一です。注目すべき例外の1つは、速筋線維が優位であるトロポニンT1ネマリンミオパチー（TNNT1）です。したがって、ネマリンミオパチーで観察される骨格筋線維の調節不全の根底にある異質性をより深く理解することは、これらの疾患と筋線維のタイプの複雑な関係を解明するのに役立つ可能性があります。
健康な対照群（1群あたりn=3）と比較して、ACTA1およびTNNT1遺伝子に変異を持つネマリンミオパチー患者から単離された筋線維は、顕著な筋線維萎縮またはジストロフィーを示しました（図6A、補足表3）。これは、利用可能な材料の量が限られているため、プロテオーム解析にとって大きな技術的課題を提示しました。それにもかかわらず、272本の骨格筋線維で2485個のタンパク質を検出することができました。線維あたり少なくとも1000個の定量化されたタンパク質をフィルタリングした後、250本の線維をその後のバイオインフォマティクス解析にかけました。フィルタリング後、線維あたり平均1573 ± 359個のタンパク質が定量されました（補足図10A、補足データセット17～18）。注目すべきことに、線維サイズが大幅に減少したにもかかわらず、ネマリンミオパチー患者サンプルのプロテオーム深度はわずかに減少しただけでした。さらに、これらのデータを独自のFASTAファイル（非コード転写産物を含む）を用いて処理することで、ネマリンミオパチー患者の骨格筋線維中に5種類の微生物タンパク質を同定することができました（補足データセット19）。プロテオームのダイナミックレンジは有意に広く、対照群の総タンパク質量は、以前の1000線維プロテオーム解析の結果とよく相関していました（補足図10B–C）。
A. ACTA1およびTNNT1型ネマリンミオパチー（NM）における線維萎縮またはジストロフィー、およびMYHに基づく異なる線維タイプの優位性を示す顕微鏡画像。スケールバー = 100 μm。ACTA1およびTNNT1患者における染色の再現性を確保するため、3つの患者生検を2～3回染色（症例ごとに4切片）し、代表画像を選択した。B. MYHに基づく参加者における線維タイプの割合。C. ネマリンミオパチー患者と対照群の骨格筋線維の主成分分析（PCA）プロット。 D. 図2で分析した1000本の線維から決定したPCAプロットに投影された、ネマリンミオパチー患者とコントロールの骨格筋線維。たとえば、ACTA1およびTNNT1ネマリンミオパチーの参加者とコントロールの間、およびACTA1およびTNNT1ネマリンミオパチーの参加者間の差異を比較するボルケーノプロット。色付きの円はπ < 0.05で有意に異なるタンパク質を示し、暗い点はFDR < 0.05で有意に異なるタンパク質を示します。統計分析は、Limma線形モデル法と経験的ベイズ法を使用して実行され、その後、Benjamini-Hochberg法を使用して多重比較のp値調整が行われました。H. プロテオーム全体およびタイプ1と2Aの線維における有意に異なる発現タンパク質のエンリッチメント分析。統計分析は、clusterProfilerパッケージとBenjamini-Hochberg調整p値を使用して実行されました。 I、J. 細胞外マトリックスとミトコンドリア遺伝子オントロジー (GO) 用語で色分けされた主成分分析 (PCA) プロット。
ネマリンミオパチーは骨格筋におけるMYH発現筋線維タイプの割合に影響を及ぼす可能性があるため19,20、まずネマリンミオパチー患者と対照群におけるMYH発現筋線維タイプを調べた。筋線維タイプは、1000筋線維アッセイで以前に報告したバイアスのない方法を用いて決定したが（補足図10D～E）、純粋な2X筋線維を同定することはできなかった（図6B）。ネマリンミオパチーが筋線維タイプに及ぼす不均一な影響を観察した。ACTA1変異を有する患者2名ではタイプ1筋線維の割合が増加したのに対し、TNNT1ネマリンミオパチーを有する患者2名ではタイプ1筋線維の割合が減少した（図6B）。実際、ACTA1-ネマリンミオパチーではMYH2および速筋トロポニンアイソフォーム（TNNC2、TNNI2、TNNT3）の発現が低下していたのに対し、TNNT1-ネマリンミオパチーではMYH7の発現が低下していました（補足図11A）。これは、ネマリンミオパチーにおける不均一な筋線維タイプのスイッチングに関する過去の報告19,20と一致しています。これらの結果を免疫組織化学染色で確認したところ、ACTA1-ネマリンミオパチー患者ではタイプ1筋線維が優位であったのに対し、TNNT1-ネマリンミオパチー患者では逆のパターンが認められました（図6A）。
単一線維プロテオームレベルでは、ACTA1およびTNNT1ネマリンミオパチー患者の骨格筋線維は、対照線維の大部分とクラスターを形成し、TNNT1ネマリンミオパチー線維は概して最も重篤な影響を受けている（図6C）。これは、各患者の疑似膨張線維の主成分分析（PCA）プロットにおいて特に顕著であり、TNNT1ネマリンミオパチー患者2および3は対照サンプルから最も離れた位置に現れた（補足図11B、補足データセット20）。ミオパチー患者の線維と健常線維の比較をより深く理解するために、健常成人参加者の1,000本の線維のプロテオーム解析から得られた詳細な情報を用いた。我々は、ミオパチーデータセット（ACTA1およびTNNT1ネマリンミオパチー患者および対照）の繊維を、1000繊維プロテオーム解析から得られたPCAプロットに投影した（図6D）。対照繊維におけるPC2に沿ったMYH繊維タイプの分布は、1000繊維プロテオーム解析から得られた繊維分布と類似していた。しかし、ネマリンミオパチー患者のほとんどの繊維は、本来のMYH繊維タイプに関わらず、PC2を下方にシフトし、健康な速筋繊維と重なっていた。つまり、ACTA1ネマリンミオパチー患者はMYHベースの方法を用いて定量化した場合、タイプ1繊維へのシフトを示したが、ACTA1ネマリンミオパチーとTNNT1ネマリンミオパチーの両方において、骨格筋繊維プロテオームは速筋繊維の方へシフトしていた。
次に、各患者群を健常者と直接比較し、ACTA1およびTNNT1ネマリンミオパチーにおいてそれぞれ256および552の発現が異なるタンパク質を同定しました（図6E～Gおよび補足図11C、補足データセット21）。遺伝子エンリッチメント解析により、ミトコンドリアタンパク質の協調的な減少が明らかになりました（図6H～I、補足データセット22）。驚くべきことに、ACTA1およびTNNT1ネマリンミオパチーにおける繊維タイプの優位性の違いにもかかわらず、この減少はMYHに基づく繊維タイプとはまったく無関係でした（図6Hおよび補足図11D～I、補足データセット23）。ACTA1またはTNNT1ネマリンミオパチーでは、3つの微生物タンパク質も制御されていました。これらのマイクロタンパク質のうちの2つ、ENSG00000215483_TR14_ORF67（LINC00598またはLnc-FOXO1としても知られる）とENSG00000229425_TR25_ORF40（lnc-NRIP1-2）は、タイプ1筋線維でのみ差別的な存在量を示しました。ENSG00000215483_TR14_ORF67は、以前に細胞周期の調節に役割を果たすことが報告されています。56 一方、ENSG00000232046_TR1_ORF437（LINC01798に相当）は、健常対照群と比較して、ACTA1-ネマリンミオパチーのタイプ1およびタイプ2A筋線維の両方で増加していました（補足図12A、補足データセット24）。対照的に、リボソームタンパク質はネマリンミオパチーによる影響をほとんど受けませんでしたが、RPS17はACTA1ネマリンミオパチーでダウンレギュレーションされました（図6E）。
エンリッチメント解析により、ACTA1およびTNNT1ネマリンミオパチーにおける免疫系プロセスのアップレギュレーションも明らかになり、TNNT1ネマリンミオパチーでは細胞接着も増加しました（図6H）。これらの細胞外因子のエンリッチメントは、細胞外マトリックスタンパク質がPC1およびPC2のPCAを負の方向（最も影響を受けた線維の方向）にシフトしたことに反映されていました（図6J）。両患者群とも、アネキシン（ANXA1、ANXA2、ANXA5）57,58やそれらの相互作用タンパク質S100A1159など、免疫応答や筋形質修復機構に関与する細胞外タンパク質の発現増加を示しました（補足図12B–C）。このプロセスは以前にも筋ジストロフィーで増強されることが報告されていますが60、われわれの知る限り、ネマリンミオパチーと関連付けられたことはありませんでした。この分子機構の正常な機能は、損傷後の筋線維鞘の修復と、新たに形成された心筋細胞と筋線維の融合に必要である58,61。したがって、両患者群におけるこのプロセスの活性上昇は、筋線維の不安定性によって引き起こされた損傷に対する修復反応を示唆している。
各ネマリンミオパチーの影響はよく相関しており（r = 0.736）、適度な重複が見られました（補足図11A–B）。これは、ACTA1とTNNT1ネマリンミオパチーがプロテオームに同様の影響を及ぼすことを示しています。しかし、一部のタンパク質はACTA1またはTNNT1ネマリンミオパチーでのみ制御されていました（補足図11AおよびC）。線維化促進タンパク質MFAP4は、TNNT1ネマリンミオパチーで最もアップレギュレーションされたタンパク質の1つでしたが、ACTA1ネマリンミオパチーでは変化がありませんでした。HOX遺伝子転写の制御を担うPAF1C複合体の構成要素であるSKIC8は、TNNT1ネマリンミオパチーでダウンレギュレーションされましたが、ACTA1ネマリンミオパチーでは影響を受けませんでした（補足図11A）。 ACTA1とTNNT1ネマリンミオパチーを直接比較したところ、TNNT1ネマリンミオパチーではミトコンドリアタンパク質の減少が大きく、免疫系タンパク質の増加が顕著であることが明らかになりました（図6G～Hおよび補足図11C、11H～I）。これらのデータは、TNNT1ネマリンミオパチーではTNNT1ネマリンミオパチーと比較して萎縮／ジストロフィーが顕著であること（図6A）と一致しており、TNNT1ネマリンミオパチーがより重篤な疾患であることを示唆しています。
観察されたネマリンミオパチーの効果が全筋肉レベルで持続するかどうかを評価するために、TNNT1ネマリンミオパチー患者の同じコホートから採取した筋肉生検のバルクプロテオーム解析を行い、対照群（1群あたりn=3）と比較しました（補足図13A、補足データセット25）。予想通り、対照群は主成分解析で密接に関連していたのに対し、TNNT1ネマリンミオパチー患者は、単一線維解析で見られるものと同様の高いサンプル間変動を示しました（補足図13B）。バルク解析では、個々の線維の比較によって強調された差次的発現タンパク質（補足図13C、補足データセット26）と生物学的プロセス（補足図13D、補足データセット27）が再現されましたが、異なる線維タイプを区別することができず、線維間の不均一な疾患効果を説明できませんでした。
これらのデータを総合すると、単一筋線維プロテオミクスは、免疫ブロット法などの標的法では検出できない臨床生物学的特徴を解明できることが示唆されます。さらに、これらのデータは、アクチン線維タイピング（MYH）のみを用いて表現型の適応を記述することの限界を浮き彫りにしています。実際、アクチンネマリンミオパチーとトロポニンネマリンミオパチーでは線維タイプのスイッチングは異なりますが、どちらのネマリンミオパチーも、MYH線維タイピングと骨格筋線維代謝を切り離し、より速く、より酸化の少ない筋プロテオームへと誘導します。
細胞の異質性は、組織が多様な要求を満たす上で極めて重要です。骨格筋では、これはしばしば、力発揮や疲労性の程度が異なる線維タイプとして説明されます。しかし、これでは骨格筋線維の多様性のほんの一部しか説明できないことは明らかです。骨格筋線維の多様性は、これまで考えられていたよりもはるかに多様で、複雑かつ多面的です。技術の進歩により、骨格筋線維を制御する因子が明らかになりつつあります。実際、私たちのデータは、タイプ2X線維が明確な骨格筋線維のサブタイプではない可能性を示唆しています。さらに、代謝タンパク質、リボソームタンパク質、細胞関連タンパク質が、骨格筋線維の異質性の主な決定因子であることを特定しました。さらに、このプロテオーム解析ワークフローを線虫性ミオパチーの患者サンプルに適用することで、MYHに基づく線維タイピングは、特にシステムが乱れた場合には、骨格筋の異質性を完全には反映しないことを実証しました。実際、MYH ベースの繊維タイプに関係なく、線虫性ミオパシーはより速く、より酸化の少ない繊維への移行をもたらします。
骨格筋線維は19世紀から分類されてきました。近年のオミクス解析により、様々なMYH線維タイプの発現プロファイルと、それらの刺激に対する反応について理解が深まりつつあります。本稿で述べたように、オミクス手法は、従来の抗体を用いた方法よりも線維タイプマーカーの定量感度が高く、単一（または少数）のマーカーの定量に依存せずに骨格筋線維タイプを定義できるという利点もあります。私たちは相補的なトランスクリプトミクスおよびプロテオミクスワークフローを用い、その結果を統合することで、ヒト骨格筋線維における線維異質性の転写および転写後制御を解析しました。このワークフローの結果、健康な若年男性コホートの外側広筋において、純粋な2X型線維をタンパク質レベルで同定することはできませんでした。これは、健康な外側広筋において純粋な2X型線維が1%未満であることが示された過去の単一線維研究と一致していますが、将来的には他の筋肉でも同様の結果が得られるはずです。 mRNAレベルではほぼ純粋な2X線維が検出されたのに、タンパク質レベルでは2A/2X線維の混合のみが検出されたという矛盾は不可解である。MYHアイソフォームmRNAの発現は概日リズムではないため67、RNAレベルでは一見純粋な2X線維におけるMYH2発現シグナルを「見逃した」可能性は低いと考えられる。考えられる説明の1つは、純粋に仮説的ではあるが、MYHアイソフォーム間のタンパク質および/またはmRNAの安定性の違いである可能性がある。確かに、どのMYHアイソフォームでも速筋線維は100%純粋ではなく、MYH1 mRNAの発現レベルが70～90%の範囲にある場合、タンパク質レベルでMYH1とMYH2の存在量が等しくなるかどうかは不明である。しかし、トランスクリプトームまたはプロテオーム全体を考慮すると、クラスター分析によって、正確なMYH構成に関係なく、骨格筋の速筋線維と遅筋線維を表す2つの明確なクラスターのみを確実に識別できる。これは、典型的には 2 つの異なる筋核クラスターのみを同定する単核トランスクリプトミクス手法を用いた解析と一致しています。68, 69, 70 さらに、これまでのプロテオーム研究ではタイプ 2X 繊維が同定されていますが、これらの繊維は他の速筋繊維とは別個にクラスター化されておらず、MYH に基づく他の繊維タイプと比較して、少量の差次的存在量のタンパク質しか示されていません。14 これらの結果は、20 世紀初頭の筋繊維分類の見解に戻るべきであることを示唆しています。この見解では、ヒトの骨格筋繊維を MYH に基づく 3 つの異なるクラスではなく、代謝特性と収縮特性に基づいて 2 つのクラスターに分類していました。63
さらに重要なのは、筋線維の不均一性は多次元的に考察されるべきであるということです。これまでの「オミクス」研究はこの方向性を示しており、骨格筋線維は個別のクラスターを形成するのではなく、連続体に沿って配列していることを示唆しています。11, 13, 14, 64, 71 本研究では、骨格筋の収縮特性や代謝特性の違いに加えて、細胞間相互作用や翻訳機構に関連する特徴によって筋線維を区別できることを示します。実際、骨格筋線維において、遅筋線維と速筋線維の種類とは無関係に、リボソームの不均一性が不均一性に寄与していることを発見しました。遅筋線維と速筋線維の種類に関わらず、この著しい筋線維の不均一性の根本原因は未だ解明されていないが、特定の力や負荷に最適に反応する筋束内の特殊な空間構成72、筋微小環境における他の細胞種との特殊な細胞特異的または臓器特異的なコミュニケーション73,74,75、あるいは個々の筋線維内のリボソーム活性の違いなどが原因として考えられる。実際、リボソームヘテロプラスミーは、RPL3およびRPL3Lの相同置換を介するか、rRNAの2'O-メチル化レベルのいずれかで、骨格筋肥大と関連することが示されている76,77。マルチオミクスおよび空間アプリケーションと個々の筋線維の機能的特徴付けを組み合わせることで、マルチオミクスレベルでの筋生物学の理解がさらに深まるであろう78。
ネマリンミオパチー患者の単一筋線維のプロテオームを解析することで、単一筋線維プロテオミクスの有用性、有効性、そして骨格筋の臨床病態生理学の解明への適用可能性も実証しました。さらに、私たちのワークフローをグローバルプロテオーム解析と比較することで、単一筋線維プロテオミクスはグローバル組織プロテオミクスと同等の詳細な情報を提供し、線維間の不均一性と筋線維の種類を考慮することで、より詳細な情報が得られることを実証しました。ACTA1およびTNNT1ネマリンミオパチーでは、健常対照群と比較して線維タイプ比に予想通り（ただし変動はあるものの）差が認められましたが、19 MYHを介した線維タイプのスイッチングとは独立して、酸化的リモデリングと細胞外リモデリングも観察されました。 TNNT1 ネマリンミオパチーでは線維化が以前にも報告されています。19 しかし、私たちの分析はこの発見に基づいており、ACTA1 および TNNT1 ネマリンミオパチー患者の筋線維において、筋線維鞘修復機構に関与するアネキシンなどの細胞外分泌ストレス関連タンパク質のレベルが上昇していることも明らかにしました。57,58,59 結論として、ネマリンミオパチー患者の筋線維におけるアネキシン レベルの上昇は、重度に萎縮した筋線維を修復するための細胞反応を表している可能性があります。
本研究は、これまでで最大規模のヒト単一線維全筋オミクス解析ではあるが、限界がないわけではない。比較的小規模で均質な参加者サンプルから単一の筋肉（外側広筋）から骨格筋線維を単離した。そのため、筋肉の種類や筋生理学的に極端な状態において、特定の線維集団が存在する可能性を排除することは不可能である。例えば、高度に訓練された短距離走者や筋力アスリート79、あるいは筋肉が活動していない期間66,80に、超高速線維のサブセット（純粋な2X線維など）が出現する可能性を排除することはできない。さらに、線維タイプの比率は男性と女性で異なることが知られているため、参加者のサンプルサイズが限られていたため、線維の不均一性における性差を調査することができなかった。さらに、同じ筋線維または同じ参加者のサンプルに対してトランスクリプトーム解析やプロテオーム解析を実施できなかった。私たちと他の研究者が、オミックス解析を使用して単一細胞および単一筋線維解析を最適化し、超低サンプル入力を実現し続けるにつれて（ここではミトコンドリアミオパチー患者の線維解析で実証されているように）、単一筋線維内でマルチオミックス（および機能的）アプローチを組み合わせる機会が明らかになります。
全体として、我々のデータは、骨格筋の異質性を引き起こす転写および転写後因子を特定し、説明するものである。具体的には、古典的なMYHに基づく線維タイプの定義に関連する、骨格筋生理学における長年の定説に疑問を投げかけるデータを提示する。我々は議論を新たにし、最終的には骨格筋線維の分類と異質性に関する理解を再考することを期待している。
白人参加者14名（男性12名、女性2名）が、自発的に本研究への参加に同意しました。本研究は、ゲント大学病院倫理委員会（BC-10237）の承認を受け、2013年ヘルシンキ宣言に準拠し、ClinicalTrials.gov（NCT05131555）に登録されました。参加者の一般的な特徴は補足表1に示されています。口頭および書面によるインフォームドコンセントを得た後、参加者は最終的に研究に参加する前に健康診断を受けました。参加者は若年（22～42歳）、健康（持病なし、喫煙歴なし）、中程度の身体活動量でした。最大酸素摂取量は、前述のように、身体能力評価用のステップエルゴメーターを用いて測定しました。81
筋生検サンプルは、安静時および絶食時に14日間間隔をあけて3回採取されました。これらのサンプルはより大規模な研究の一環として採取されたため、参加者は生検の40分前にプラセボ（乳糖）、H1受容体拮抗薬（フェキソフェナジン540mg）、またはH2受容体拮抗薬（ファモチジン40mg）を摂取しました。これらのヒスタミン受容体拮抗薬は安静時の骨格筋フィットネスに影響を与えないことを我々は以前に実証しており81、品質管理プロットにおいて状態関連のクラスタリングは観察されませんでした（補足図3および6）。各実験日の48時間前から標準食（体重1kgあたり41.4 kcal、炭水化物5.1 g/kg、タンパク質1.4 g/kg、脂肪1.6 g/kg）を摂取し、実験日の朝には標準朝食（体重1kgあたり1.5 gの炭水化物）を摂取した。局所麻酔（エピネフリンを含まない1%リドカイン0.5 ml）下で、経皮的ベルグストローム吸引法を用いて外側広筋から筋生検を実施した。82 筋サンプルは直ちにRNAlaterに包埋し、用手線維解剖（最大3日間）まで4℃で保存した。
新鮮に単離した筋線維束を培養皿内の新鮮なRNAlater培地に移した。次に、実体顕微鏡と細いピンセットを用いて、個々の筋線維を手作業で解剖した。各生検組織から25本の線維を解剖し、生検組織の異なる部位から線維を選択するように特に注意を払った。解剖後、各線維をプロテイナーゼKおよびDNase酵素を含む溶解バッファー（SingleShot Cell Lysis Kit、Bio-Rad）3 μlに静かに浸し、不要なタンパク質とDNAを除去した。次に、軽くボルテックスミキサーで撹拌し、液体をマイクロ遠心分離機で遠心分離し、室温で10分間インキュベートすることで、細胞溶解およびタンパク質/DNAの除去を開始した。次に、溶解液をサーマルサイクラー（T100、Bio-Rad）で37°Cで5分間、75°Cで5分間インキュベートし、次の処理まで直ちに-80°Cで保存した。
QuantSeq-Pool 3′ mRNA-Seq Library Prep Kit（Lexogen社）を用いて、2 µlの筋線維ライセートからイルミナ対応ポリアデニル化RNAライブラリを調製しました。詳細な手順については、製造元のマニュアルをご覧ください。このプロセスは、逆転写によるファーストストランドcDNA合成から始まります。この過程で、ユニーク分子識別子（UMI）とサンプル固有のi1バーコードが導入され、サンプルのプールが確保され、下流処理における技術的なばらつきが低減されます。次に、96本の筋線維からcDNAをプールし、磁気ビーズで精製します。その後、RNAを除去し、ランダムプライマーを用いてセカンドストランド合成を行います。ライブラリは磁気ビーズで精製され、プール固有のi5/i7タグが付加され、PCR増幅されます。最後の精製ステップでイルミナ対応ライブラリが生成されます。各ライブラリ プールの品質は、高感度小断片 DNA 分析キット (Agilent Technologies、DNF-477-0500) を使用して評価されました。
Qubit定量に基づき、プールはさらに等モル濃度（2 nM）でプールされました。得られたプールは、NovaSeq 6000システムの標準モードで、NovaSeq S2試薬キット（1 × 100ヌクレオチド）を用いて2 nMローディング（4% PhiX）でシーケンスされました。
このパイプラインは、LexogenのQuantSeq Poolデータ分析パイプライン（https://github.com/Lexogen-Tools/quantseqpool_analysis）に基づいています。データはまず、i7/i5インデックスに基づいてbcl2fastq2（v2.20.0）でデマルチプレックスされました。次に、リード2はi1サンプルバーコードに基づいてidemux（v0.1.6）でデマルチプレックスされ、umi_tools（v1.0.1）でUMI配列が抽出されました。その後、cutadapt（v3.4）を使用して複数回のトリミングを行い、短いリード（長さ20未満）またはアダプター配列のみで構成されるリードを削除しました。その後、STAR（v2.6.0c）を使用してリードをヒトゲノムにアライメントし、SAMtools（v1.11）を使用してBAMファイルをインデックス化しました。umi_tools（v1.0.1）を使用して重複リードを削除しました。最後に、Subread (v2.0.3) の featureCounts を使用してアライメントカウントを実行しました。パイプラインの複数の中間段階では、FastQC (v0.11.9) を使用して品質管理を実行しました。
バイオインフォマティクス処理と可視化は、主にSeurat（v4.4.0）ワークフローを用いてR（v4.2.3）で行われた。83 そのため、個々のUMI値とメタデータマトリックスはSeuratオブジェクトに変換された。全繊維の30%未満で発現する遺伝子は除去された。低品質サンプルは、UMI値1000個と検出された遺伝子1000個という最小閾値に基づいて除去された。最終的に、925本の繊維がすべての品質管理フィルタリングステップを通過した。UMI値は、検出された7418個の特徴すべてを含むSeurat SCTransform v2メソッド84を用いて正規化され、参加者間の差異は回帰分析によって除去された。関連するすべてのメタデータは、補足データセット28にある。