We define matrix q-analogues of the higher-order Mersenne numbers via holomorphic functional calculus. For \(q\in (0,1),\) \(k\in \mathbb {N},\) and a complex square matrix A, let \(B_{k,q}(A):=[A^k]_q=(\varvec{I}-q^{A^k})/(1-q)\) and set \(M^{(k)}_{n,q}(A)=\sum _{j=0}^{n-1}B_{k,q}(A)^j.\) All single-matrix identities are polynomial or holomorphic identities in \(B_{k,q}(A)\) and hence hold for arbitrary A, while commutativity is only required in genuinely multi-parameter matrix settings. We derive Binet-type formulas, matrix recurrences, a factorized Jackson q-exponential generating function, and a -embedding. A transfer-matrix method yields matrix q-Cassini, q-Catalan, and q-d’Ocagne identities, together with a determinantal Abel–Cassini identity for the nonautonomous system. We also prove spectral mapping results for \(\sigma \!\big (M^{(k)}_{n,q}(A)\big )\) and include explicit \(2\times 2\) examples, including a Jordan-block case. As \(q\rightarrow 1^-,\) the scalar specialization \(r=1\) and \(A=2\) recovers the higher-order Mersenne formulas of Prasad et al. [9].