Does the sample size influence the number of PCs needed to explain a fixed percentage of variance? In PCA, I have observed the following pattern: larger sample sizes from the same population(s) imply that a higher number of principal components (PCs) are required to explain a fixed percentage of variance. 
Using an anecdotal example, if I have a sample size of 100, then 10 PCs explain 90% of the variance. However, if I have a sample size of 200, then 20 PCs may be required to explain the same 90% of variance.
My questions are: 


*

*Does this pattern hold generically for all kinds of data?

*If the answer to the previous question is "yes", then is there any literature reference I can use to make such claim? (I've searched but didn't find any).

*If the answer to the previous question is "no", is there a simple rational I can use to explain why this pattern occurs?

 A: You are probably in the $n\ll p$ situation. As @Glen_b wrote in the comment above, such a behaviour is expected.
Illustration and intuition
Let $p=1000$. Let the population be a multivariate normal with mean zero and some arbitrary covariance matrix $\boldsymbol\Sigma$. Let us sample $n$ points from this population and compute the sample covariance matrix $\mathbf S_n$. Here I plotted the spectrum (sorted eigenvalues) of $\boldsymbol\Sigma$ (black line) and the spectra of $\mathbf S_n$ for various values of $n$:

The sum of all eigenvalues (trace of $\mathbf S_n$) remains approximately constant because it is equal to the sum of variances of each of the $p$ variables, and those can be reasonably well estimated already with small $n$: $$\operatorname{tr}(\mathbf S_n)\approx \operatorname{tr}(\boldsymbol\Sigma).$$ But if $n=100$, then certainly only $100$ eigenvalues can be non-zero, so the same trace has to be "spread" over only $100$ values, meaning that the leading eigenvalues will be much larger than the population ones. As $n$ grows, the leading eigenvalues will decrease and the tail will grow. Notice that even once $n>p$, the bias still remains and only with $n=10000$ (ten times the dimensionality) the spectrum starts to look like the population spectrum.
The dots mark the number of components that explain $90\%$ of the variance. The more horizontal the spectrum looks, the larger this number, so by now it should be clear that it will increase with increasing $n$.
For clarity, here is the same example with $\boldsymbol\Sigma = \mathbf I$. The same effect can still be clearly observed:

Theory
If $\mathbf x \sim \mathcal N (0, \boldsymbol \Sigma_{p\times p})$, then $$(n-1) \mathbf S_n \sim \mathcal W_p(\boldsymbol \Sigma, n-1),$$ where $\mathcal W$ is Wishart distribution. Wishart distribution is well studied, so I expect that there are results about the sampling properties of the eigenvalues of Wishart-distributed matrices. I am not familiar with this field so I cannot say much more, but this would be a starting point for further exploration.
References
I think you might want to cite Jollife, 2002, Principal Components Analysis, section 3.6 "Probability Distributions for Sample Principal Components". Here is what he writes there, page 48 ($l_i$ denote eigenvalues of $\mathbf S_n$  and $\lambda_i$ denote eigenvalues of $\boldsymbol \Sigma$):

One specific point [...] is that $E(l_1 ) > \lambda_ 1$ and
  $E(l_p ) < \lambda_p$. In general the larger eigenvalues tend to be overestimated and the smaller ones underestimated.

He also adds that

If a distribution other than the multivariate normal is assumed, distributional results for PCs will typically become less tractable.

The references given around are Jackson, 1991, A User’s Guide to Principal Components and Srivastava and Khatri, 1979, An Introduction to Multivariate Statistics. I am not familiar with these books. 
Note that Jolliffe does not explicitly comment on the fact the you need more PCs to explain a certain percentage of variance. But perhaps you can write something like that:

For smaller $n$ less PCs are needed to explain the same amount of variance because for smaller $n\ll p$ leading eigenvalues tend to be overestimated and trailing eigenvalues underestimated (Jollife 2002, Section 3.6).


Matlab code to produce these figures
clear all
p = 1000;
ns = [100 200 500 1000 5000 10000];

W = randn(p,p);
Sigma = transpose(W)*W;
%// alternatively: Sigma = eye(p);
spectrum_population = sort(eig(Sigma), 'descend');

figure('Position', [100 100 1000 400])
hold on
col = lines(length(ns));

for i = 1:length(ns)
    X = randn(ns(i),p);
    X = X * chol(Sigma);
    spectra(i,:) = sort(eig(cov(X)), 'descend');

    h(i) = plot(spectra(i,:), 'Color', col(i,:));
    ind = find(cumsum(spectra(i,:)) > 0.9*sum(spectra(i,:)), 1);
    plot(ind, spectra(i,ind), '.', 'MarkerSize', 20, 'Color', col(i,:))

    leg{i} = ['n = ' num2str(ns(i))];
end

h(length(ns)+1) = plot(spectrum_population, 'k', 'LineWidth', 2);
ind = find(cumsum(spectrum_population) > 0.9*sum(spectrum_population), 1);
plot(ind, spectrum_population(ind), 'k.', 'MarkerSize', 20)
leg{length(ns)+1} = 'Population';

legend(h, leg)
legend boxoff 
axis([0 p 0 max(spectra(:))])

