In order to apply Nyquist-Shannon sampling theory, we need to know the maximum frequency that will be present in the signal we intend to measure. We will do this by rewriting a time-dependent observable in terms of the frequencies present in the Hamiltonian $H$.
Consider an arbitrary observable $O$ which will time evolve under application of $H$ in the Heisenberg picture (with $\hbar = 1$) according to
\begin{align}
O(t) &= e^{i H t} O e^{-i H t} \tag{1}
\\&= \sum_{j, k} e^{i(\omega_j - \omega_k) t} \langle \omega_j|O |\omega_k\rangle |\omega_j \rangle \langle \omega_k | \tag{2}
\end{align}
where for simplicity I consider the finite-dimensional case such that $H$ admits a spectral decomposition
\begin{equation}
H = \sum_{k} \omega_k |\omega_k\rangle\langle\omega_k |\tag{3}
\end{equation}
Take any initial state $|\psi_0\rangle$ and define $c_{jk} = \langle \omega_j|O |\omega_k\rangle \langle \psi_0|\omega_j \rangle \langle \omega_k | \psi_0\rangle$ for housekeeping, and we find that the expected value of $O(t)$ is
\begin{align}
\langle O(t) \rangle &= \sum_{j,k} c_{jk} e^{i(\omega_j - \omega_k) t}\tag{4}
\end{align}
Then defining the discrete set of frequency differences
\begin{align}
\Delta = \{\omega_j - \omega_k: \omega_k,\omega_j \in \text{spec}(H)\}\tag{5}
\end{align}
we can gather terms of the same frequency difference to rewrite (4) as
\begin{align}
\langle O(t) \rangle &= \sum_{\delta \in \Delta} \left(\sum_{\substack{j,k: \\\omega_j - \omega_k = \delta}}c_{jk}\right) e^{i\delta t} \tag{6}
\\&:= \sum_{\delta \in \Delta} a(\delta)e^{i\delta t}\tag{7}
\end{align}
This provides the expected value as a Fourier decomposition, and we can identify the highest frequency present as $\max (\Delta) = \omega_{max}(H) - \omega_{min}(H)$ (no need to worry about the sign since by construction $\delta \in \Delta$ implies $-\delta \in \Delta$). Then applying Nyquist-Shannon sampling theorem we find that in order to reconstruct $\langle O(t) \rangle$ it is sufficient to measure $O$ at a fequency $\omega_s$ satisfying
\begin{equation}
\omega_s \geq 2\max(\Delta) \tag{8}
\end{equation}
which of course will need to be performed over many different trials of preparing $|\psi_0\rangle$ and measuring $O$.
Applying this to your question, if you are trying to reconstruct the complete wavefunction $|\psi (t)\rangle$ then yes, you should be able to do so by constructing each element in the set
\begin{equation}
\{|\psi (0)\rangle, |\psi(\Delta t)\rangle, |\psi(2 \Delta t)\rangle, \dots |\psi(N \Delta t)\rangle\}\tag{9}
\end{equation}
at a time interval $\Delta t \leq 1 / f_s = 2 \pi / \omega_s$. This follows because determining the elements of $|\psi (t)\rangle$ (using, e.g. quantum state tomography) involves a series of experiments measuring different observables. None of these observables can have time dependence with a frequency greater than $\max (\Delta)$ and therefore no composition of these observables needed for the state tomography can have a frequency greater than $\max (\Delta)$ either. Of course, tomography experiments are generally expensive as they require an experimental overhead scaling like $O(\text{dim}(H))$. Though maybe there might be savings for that scaling in the context of this experiment where each state in the sequence is somehow related to the previous.
On the other hand I'm not sure if a sampling rate $1 / f_s$ is necessary for every set of observables; it may be that some sets of observables exhibit different maximum frequencies in their time dependence (all upperbounded by $\max(\Delta)$ of course) and so this might be interesting to look into further.
Finally, yes you will need to be concerned with shot noise for the estimator $\tilde{O}(t)$ of $\langle O(t)\rangle$, as the deviations of $\tilde{O}(t)$ around the true mean can possibly introduce much higher frequencies to the spectrum of $\tilde{O}(t)$ compared to $\langle O(t)\rangle$.
