It depends on the true covariance structure of population. If multiple features have the same population variance, sampling variance can actually mix up the observed features arbitrarily, whereas if all features have different population variances this cannot happen.
Let me show this through a derivation. Assume your true population covariance matrix is $A$ and your observed is $A + \epsilon B$, where $\epsilon$ is a small positive number. Basically, your can think of $B$ as the direction that the covariance matrix is perturbed in (perturbation meaning 'error' due to sampling), and $\epsilon$ is the magnitude of that perturbation.
Remember that PCA is basically looking at the eigenvectors of the covariance matrix.
If the true covariance is $A = aI$, the eigenvectors of $A + \epsilon B$ are just eigenvectors $v$ of $B$. In other words, running PCA on your observed data will return features determined entirely by your sampling variation!
The reason this happens is that PCA does not have a unique solution for your true covariance matrix in the first place. That is, whenever you try to diagonalize a matrix with duplicate eigenvalues, you run into the issue that there are multiple valid ways to pick unit eigenvectors (because you have eigenspaces that have dimension greater than one). Thus even if $A$ is not of the form $aI$, but has duplicate eigenvalues, there will be some set of eigenvectors of $A + \epsilon B$ which will deviate from the original features in a way dictated entirely by the eigenvectors of $B$.
Fortunately, in reality it is rare to find variables with exactly the same population variance. In this case, we can show that the features are robust against sampling variation. This is actually true in general (i.e. regardless of the assumption that your population variables are uncorrelated).
To see this, we can just approximate the eigenvectors of $(A + \epsilon B)$ under some mild assumptions. Basically we want to show that they are just a perturbation of the eigenvectors of $A$ by a term that scales roughly linearly with $\epsilon$.
The first assumption is of course that all eigenvalues of $A$ are unique. However, second, we need to assume that the eigenspaces of $A$ and $B$ are disjoint, meaning that they have different eigenvectors. Again we can just appeal to reality -- a random matrix will almost never have the same eigenvectors as a fixed matrix (unless that fixed matrix is a multiple of $I$).
A warning that this "proof" is both very long and also not exactly valid. The key weakness is where I assume that the solution is analytic, which is sort of what we actually want to prove. Regardless, hopefully it will provide some insight.
[start of proof]
Let $u$ denote the $i^{th}$ eigenvector of $A$ for some $i$, and have eigenvalue $\alpha$. We can arbitrarily write the eigenvector of $(A + \epsilon B)$ as $(u+v)$ for some $v$ by taking an eigenvector of the matrix and then letting $v$ be that eigenvector minus $u$.
Now, what is $(A + \epsilon B)(u+v)$? It should be $\lambda (u + v)$ for some $\lambda$.
For simplicity, we can rewrite this as $(\alpha + \epsilon \gamma)(u+v)$ for $\gamma = \epsilon^{-1}(\lambda - \alpha)$. This makes sense because the eigenvalues are a continuous function of the matrix entries, which can be proven by using that the determinant is a continuous function (since it is a polynomial) and then using that the eigenvalues are defined using the determinant.
On the other hand, we also have
$$ (A + \epsilon B)(u+v) = Au + Av + \epsilon B (u+v) = \alpha u + Av + \epsilon B (u+v). $$
Thus
\begin{align*}
(A + \epsilon B)(u+v)
&= (\alpha + \epsilon \gamma)(u+v) = \alpha u + \epsilon \gamma u + \alpha v + \epsilon \gamma v \\
Av + eBu + eBv
&= \epsilon \gamma u + (a + \epsilon \gamma)v \\
\epsilon B u + (A + \epsilon B)v
&= \epsilon \gamma u + (\alpha + \epsilon \gamma)v \\
\epsilon(B - \gamma I)u + ((A - \alpha I) + \epsilon (B - \gamma I))v
&= 0
\end{align*}
Now, there are two possible ways this equation could hold. One involves having $(A - \alpha I)v = 0$. However, this is would mean that $v$ is an eigenvector of $A$ with eigenvalue $\alpha$, and thus so is $u+v$. The problem with this is that it would imply tht $(u+v)$ is an eigenvector of $B$ as well, and by assumption $B$ cannot share eigenvectors with $A$. Thus certainly $(A - \alpha I)v$ is nonzero.
However, looking at the equation, the $(A - \alpha I)v$ is the only term that does not appear proportional to $\epsilon$. Why is this interesting? Well, because if it were not dependent on $\epsilon$ at all, the equation could not have a solution -- since $\epsilon$ is essentially arbitrary, taking the limit as $\epsilon \to 0$ would create an inconsistent equation. The fix is to recognize that $v$ depends on $\epsilon$, and really it must be at least proportional. More precisely, it should be
$$ v = w_{0} + \epsilon w_{1} + \epsilon^{2} w_{2} + (\mathrm{higher \ order \ terms}). $$
I am making a bold assumption that $v$ is essentially analytic in $\epsilon$, but since we are basically solving a polynomial equation ($\gamma$ is also analytic in $\epsilon$), it seems reasonable. This is where the approximation comes in anyway.
Now, note that the $w_{0}$ must be zero, since as $\epsilon \to 0$, we must have $u+v \to u$.
Basically, if $\epsilon$ is quite small, then $v \approx \epsilon w$. Furthermore,
\begin{align*}
\epsilon(B - \gamma I)u + ((A - \alpha I) + \epsilon (B - \gamma I))v
&= \epsilon(B - \gamma I)u + ((A - \alpha I) + \epsilon (B - \gamma I))\epsilon w \\
&= \epsilon(B - \gamma I)u + \epsilon (A - \alpha I) w
+ \epsilon^{2} (B - \gamma I) w \\
&\approx \epsilon (B - \gamma I)u + \epsilon (A - \alpha I)
\end{align*}
because $\epsilon^{2} \approx 0$.
Therefore, dividing out by $\epsilon$, in the end we are just solving
$$ (B - \gamma I)u + (A - \alpha I) w = 0 $$
or really
$$ (B - \gamma I)u = -(A - \alpha I) w. $$
Now, $-(A - \alpha I) w$ is some unknown vector in the columnspace of $(A - \alpha I)$.
Knowing this allows to solve for $\gamma$. Once we have done this, the term $(B - \gamma I)u$ becomes a known vector and we are just solving a linear equation.
The value of $\gamma$ is unique because the columnspace of $(A - \alpha I)$ has smaller dimension than the whole space. In other words, it is 'difficult' to get the image of $u$ exactly into the subspace, so there is only one way to do it. However, the solution for $w$ is not unique, in the sense that $(\alpha I - A)w = y$ does not have a unique solution ($y = (B - \gamma I)u$), with the reason being that $(\alpha I - A)$ has a nontrivial nullspace, in particular containing the eigenvectors of $A$ with eigenvector $\alpha$. However, we can pick the vector $w$ solving our linear equation which minimally violates the equation we actually wanted to solve. In other words, we can pick it so the norm of $(B - \gamma I)w$ is minimal. This is unique.
[end of proof]
Ok, so what have we actually shown with all of this? Basically, for small $\epsilon$ it is possible to get a unique approximation to the eigenvectors of $(A + \epsilon B)$ that has minimal 'error'. These eigenvectors are perturbed from the eigenvectors of $A$ by a small amount which is proportional to $\epsilon$. Therefore, so long as the features of the true population would be uniquely chosen by PCA (i.e. they have distinct population variances), the features of the observed data also can be uniquely chosen by PCA and they are perturbed from the true features by an amount roughly proportional to the size of the sampling error (assuming the sampling error is small!).
