Update: the extension is correct since for $j = 1, \dots, k$
$$
\begin{align}\sigma_j(M_1^\mathrm{T} A M_2) &= \sigma_j\left(\begin{bmatrix}M_1& O_{n, n-k}\end{bmatrix}^\mathrm{T} A \begin{bmatrix}M_2& O_{n, n-k}\end{bmatrix}\right)\\
& \leq \sigma_j(A) \cdot \left\|\begin{bmatrix}M_1& O_{n, n-k}\end{bmatrix}\right\| \cdot \left\|\begin{bmatrix}M_2& O_{n, n-k}\end{bmatrix}\right\| \\
&= \sigma_j(A).
\end{align}
$$

```
Below is an earlier proof.
```

By considering the SVD of $A$, it is w.l.o.g. to assume that $A = D =\mathrm{diag}(d_1, \dots, \dots, d_n)$ with $d_1 \geq \dots \geq d_n \geq 0$. Suppose we have the SVD $M_1^\mathrm{T} D M_2 = U \hat{D} V^\mathrm{T}$, where $\hat{D} = \mathrm{diag}(\hat{d}_1, \dots, \hat{d}_k)$, $U$ and $V$ are orthogonal matrices. Let us define $n \times k$ matrices $Q = M_1 U, S = M_2V$. Then we have $Q^\mathrm{T}Q = S^\mathrm{T}S = I_k$ and
$\hat{D} = Q^\mathrm{T} D S$, which implies that
$$
\hat{d}_j = \sum_{i = 1}^n q_{ij} s_{i j} d_i, \quad j = 1, \dots, k,
$$
where $q_{ij}$ and $s_{i j}$ are the $(i, j)$-th entries of $Q$ and $S$, respectively.

As a result,
$$
\begin{align}
\sum_{j = 1}^k \left[\sigma_j\left(M_1^\mathrm{T}DM_2\right)\right]^p
= \sum_{j = 1}^k (\hat{d}_j)^p
&= \sum_{j = 1}^k \left(\sum_{i = 1}^n q_{ij} s_{i j} d_i\right)^p\\
&\leq \sum_{j = 1}^k \left(\sum_{i = 1}^n |q_{ij} s_{i j}| d_i\right)^p\\
&= \sum_{j = 1}^k \left(\sum_{i = 1}^n \frac{|q_{ij} s_{i j}|}{\sum_{i = 1}^n |q_{ij} s_{i j}|} d_i\right)^p \cdot \left(\sum_{i = 1}^n |q_{ij} s_{i j}|\right)^p\\
&\leq \sum_{j = 1}^k \sum_{i = 1}^n \frac{|q_{ij} s_{i j}|}{\sum_{i = 1}^n |q_{ij} s_{i j}|} d_i^p \cdot \left(\sum_{i = 1}^n |q_{ij} s_{i j}|\right)^p\\
&=\sum_{j = 1}^k \sum_{i = 1}^n |q_{ij} s_{i j}| d_i^p \cdot \left(\sum_{i = 1}^n |q_{ij} s_{i j}|\right)^{p-1}\\
&\leq \sum_{j = 1}^k \sum_{i = 1}^n \left(\frac{q_{ij}^2 + s_{ij}^2}{2}\cdot d_i^p \cdot \left(\sum_{i = 1}^n \frac{q^2_{ij}+ s^2_{i j}}{2}\right)^{p-1}\right)\\
&= \sum_{i = 1}^n c_i\cdot d_i^p,
\end{align}
$$
where we have defined
$\displaystyle c_i = \sum_{j = 1}^k \frac{q_{ij}^2 + s_{ij}^2}{2}, i= 1, \dots, n$. To obtain the second inequality, we have used the convexity of $x \mapsto x^p$. The last equality follows from the fact $$ \sum_{i = 1}^n q_{ij}^2 = \sum_{i = 1}^n s_{ij}^2 = 1, \quad\text{for } j = 1, \dots, k.$$
Since $0 \leq c_i \leq 1$ for each $i = 1, \dots, n$, and $\displaystyle \sum_{i = 1}^n c_i = \|Q\|_{\mathrm{F}}^2 + \|S\|_{\mathrm{F}}^2 = k$, we have
$$\sum_{i = 1}^n c_i\cdot d_i^p \leq \sum_{i = 1}^k d_i^p = \sum_{i = 1}^k \left[\sigma_i\left(D\right)\right]^p,$$ which completes the proof.