Short story: if $\log n \ll k < c n$ for some constant $c<1$ then $$ \|P_V A P_V\|_\mathrm{op} \asymp \operatorname{trace}[A]/n + \|A\|_\mathrm{op}k/n $$ up to multiplicative constants dependent on $c$.
By rotational invariance, assume without loss of generality that $A=\operatorname{diag}(a_1,\ldots,a_n)$ with $a_1\ge a_2 \ge \ldots \ge a_n$. The object of study is the operator norm $\|P_V A P_V\|_\mathrm{op}$ where $P_V$ is an orthogonal projector on the Grassmannian. Equivalently, the goal is to study $\|\sqrt{A} P_V\|_\mathrm{op}$.
It is easier to work with Gaussian matrices as more random matrix theory results are available. Realize $P_V$ as $X(X^TX)^{-1}X^T$ where $X\in R^{n\times k}$ has iid $N(0,1/n)$ entries.
I will focus on the situation $k\le c n$ for constant $c<1$. Then $(X^TX)^{-1}X^T$ has bounded singular values with overwhelming probability in the sense $(c_*\le \text{singular values} \le c^*$ for constants $c_*,c^*$ depending on $c$ only), see for instance Theorem II.13 in Davidson and Szarek (2000), although that reference has a missing $\sqrt n$ in the statement of the theorem. This means that up to multiplicative constants depending on $c$, we may as well characterize $\|\sqrt AX\|_\mathrm{op}$. This does not lose much for the approach below, which characterizes $\|\sqrt AX\|_\mathrm{op}$ up to multiplicative constants anyway.
The expected operator norm of inhomogeneous Gaussian matrices with independent entries has been fully characterized in [1]. In the symmetric case [1, Theorem 1.1], if $b_{ij}\ge 0$ and the matrix has entries $M_{ij} = b_{ij} g_{ij}$ for $g_{ij} \sim N(0,1)$ then $$E\|M\|_\mathrm{op}\asymp\max_i \sqrt{\sum_j b_{ij}^2} + \max_{ij} b_{ij}^*\sqrt{\log i}$$ up to multiplicative universal constants independent of the dimensions and the $\{b_{ij}\}$; above $b_{ij}^*$ are the entries of the matrix $(b_{ij})$ obtained by reordering the rows/columns such that $\max_i b_{1i}^* \ge \max_i b_{2i}^* \ge ...$.
In our case, $\sqrt A X$ has independent normal entries but is not symmetric, so we look instead at $$ \begin{pmatrix} 0 & \sqrt{A}X \\ X^T \sqrt{A} & 0 \end{pmatrix} $$ which is symmetric and has the same operator norm as $\sqrt{A} X$ up to multiplicative constants (this symmetrization device is explained in the later sections of [1]). Then $$\max_i \sqrt{\sum_j b_{ij}^2} \asymp \sqrt{\operatorname{trace}[A]/n} + \sqrt{\|A\|_\mathrm{op} k/n}, $$ and the term $\max_{ij} b_{ij}^* \sqrt{\log i}$ is dominated by $\sqrt{\|A\|_\mathrm{op} k/n}$ if $k\gg\log n$.
Conclusion: $\|P_V A P_V\|_\mathrm{op}$ is of order $\operatorname{trace}[A]/n + \|A\|_\mathrm{op} k/n$ up to multiplicative constants depending on the constant $c<1$ whenever $\log n \ll k \le c n$. I left under the rug the issue of concentration of $\|\sqrt{A} X\|_\mathrm{op}$ around its expectation which should be easy to obtain using the standard concentration of Lipschitz functions of independent normals (see for instance Theorem 5.6 in the book by Boucheron, Lugosi and Massart).
[1]: The dimension-free structure of nonhomogeneous random matrices (with Rafał Latała and Pierre Youssef) Invent. Math. 214, 1031-1080 (2018).