Strictly speaking, MCMC means constructing a Markov chain which has some desired property (usually either a prespecified equilibrium distribution or a prespecified expectation) and then numerically drawing samples from it. If you are given a Markov chain it is really just called Monte Carlo simulation.
Anyway, at each fixed time $t$, the value of the chain will have some given distribution, given by $p_0 P^t$, where $p_0$ is the initial distribution (written as a row vector) and $P$ is the transition probability matrix (in the row-stochastic convention i.e. $\sum_{j=1}^n P_{ij}=1$).
So the long term behavior of the distributions at each fixed $t$ is determined by the behavior of $P^t$ as $t \to \infty$. Any stochastic matrix has all eigenvalues of modulus at most $1$. It also has $1$ as an eigenvalue. The ergodicity assumption additionally tells us that the eigenvalue $1$ has multiplicity $1$ and that all of its other eigenvalues will have modulus strictly less than $1$. (Without irreducibility it is possible that the eigenvalue $1$ has a higher multiplicity; without aperiodicity it is possible to have an eigenvalue like $-1$ or $e^{2\pi i/3}$.) The reversibility assumption tells you (among other things) that $P$ is diagonalizable with real eigenvalues (because reversibility furnishes a similarity transformation which turns $P$ into a symmetric matrix).
As a result of these three facts you can write
$$p_0 P^t = \pi + \sum_{i=2}^n c_i \lambda_i^t q_i$$
where $\lambda_i$ are eigenvalues all less than $1$ in absolute value, $q_i$ are the left eigenvectors, and $c_i$ are coefficients depending on $p_0$. This is just what you get from diagonalizing $P^t$ (without going through and computing the $c_i$ since they do not really matter for this presentation).
You can see that as $t \to \infty$, the remaining eigenvalues will decay because of the power of $t$. Assuming the $c_i$ are all significant (the typical situation for a "random" initial distribution), the slowest decay is by $\lambda_{\max}$, since it is closer to $1$ in absolute value. So the larger the spectral gap $\xi$, the faster the convergence. Note that the negative eigenvalues do matter, which is why you wrote it as $\lambda_{\max}=\max \{ |\lambda_2|,|\lambda_n| \}$.
The precise relationship between the spectral gap and the dynamical properties of the chain is more subtle than the rest of what I've said here. Loosely speaking when the spectral gap is small, the eigenvector corresponding to $\lambda_{\max}$, which I might call a "metastable eigenvector", will be significant and positive for one set of states and significant and negative for another set of states. It might also be smaller on a third set of states.
The "metastable eigenvector" then corresponds to the slow exchange of probability between the first two sets of states (because the decay of its contribution to the distribution results in a flow from one set to the other, with which set is which depending on the sign of the corresponding $c_i$). You can look up metastability in Markov chains for more about this. To get a feel for it, consider a matrix like
$$P=\begin{bmatrix} 1/3 & 2/3 & 0 & 0 \\
1/4-\epsilon/2 & 3/4-\epsilon/2 & \epsilon & 0 \\
0 & \epsilon & 1/5-\epsilon/2 & 4/5-\epsilon/2 \\
0 & 0 & 1/6 & 5/6 \end{bmatrix}$$
where $0<\epsilon \ll 1$ is a small parameter. For $\epsilon=10^{-6}$ for example, this has relatively rapid equilibration within the sets $\{ 1,2 \}$ and $\{ 3,4 \}$ (driven by eigenvalues of about $1/12$ and $1/30$ respectively) and slow equilibration between the two of them (driven by an eigenvalue of about $1-10^{-6}$).
Best Answer
But suppose we are still guaranteed that $\big \vert \lambda_1\big \vert \gt \big \vert \lambda_2\big \vert \geq \ldots \geq \big \vert \lambda_n\big \vert$ where the $\lambda_i$'s are all the complex eigenvalues of $P$, then can we say that the given irreducible chain is aperiodic
Yes. This is a basic spectral result, or if you prefer a basic Perron Frobenius theory result. Since we know $\lambda_1 = 1$, i.e. the Perron root, suppose you look at your transition matrix $P$ and 'remove' the perron vectors, i.e. consider the matrix $B := P- \mathbf {1\pi}^{\top}$. The matrix $B$ is what determines the difference between the distribution in your chain and the stationary distribution.
$B$ has maximal modulus eigenvalue $\big \vert \lambda_2 \big \vert \lt 1$. So $\lim_{k \to \infty} B^k \to \mathbf 0$ and $\lim_{k \to \infty} P^k =\mathbf {1\pi}^{\top}$. This limit cannot exist if $P$ is periodic. So $P$ is aperiodic. Finite state chains always have at least one positive recurrent class, with an associated stationary distribution; so the non-existence of a limit is equivalent to a given chain having periodicity.
Equivalently a finite state (time homogenous) markov chain with a single communicating class is periodic iff it only has one eigenvalue on the unit circle.
Addendum
here is a sketch of a different route to the above result, one which relies on some probabilistic and analytic results and very little linear algebra. Everything below assumes there is a single communicating class.
1.) Prove that for any Markov chain with at least one non-zero entry on the diagonal of the n x n transition matrix $A$, you have $A^k \gt 0$ for all $k \geq K$ -- where this inequality is interpreted component-wise -- i.e. all components of $A^k$ are strictly positive. You can tighten the bound but an easy proof would be to select say $K = 5 n^2$ -- all that matters is K is some constant that depends only on $n$.
2.) Since looking at the diagonal gives useful information, consider that
$\text{trace}\big(P^r\big)$
$ = \lambda_1^r + \lambda_2^r + ... +\lambda_n^r $
$= 1 + \sum_{j=2}^n \lambda_j^r $
$= \big \vert 1 + \sum_{j=2}^n \lambda_j^r \big \vert$
$\geq \Big \vert \big \vert 1\big \vert -\big \vert \sum_{j=2}^n \lambda_j^r \big \vert \Big \vert $
$\geq 1 -\big \vert \sum_{j=2}^n \lambda_j^r \big \vert $
$\geq 1- \sum_{j=2}^n \big \vert\lambda_j \big \vert^r $
$\geq 1- (n-1)\cdot \big \vert\lambda_2\big \vert^r$
by application of triangle inequality
and for large enough $r$ we have
$\text{trace}\big(P^r\big) \geq 1- (n-1)\cdot \big \vert\lambda_2\big \vert^r \gt 0$, i.e. $P^r$ has at least one positive component on its diagonal for all $r\geq R$.
In combination with (1) this implies that $P^m \gt 0$ for all $m \geq M$.
3.) Since for large enough powers of $P$, the transition matrix is strictly positive, we can show that $P^m \to \mathbf {1\pi}^{\top}$. For several different elementary takes consider section 11.4 Fundamental Limit Theorem for Regular Chains in the free book by Grinstead and Snell
-- in particular (i) the main proof that starts the section, (ii) the Doeblin Proof, and (iii) the exercise 9 suggested by Doyle. All three are quite nice.
https://math.dartmouth.edu/~prob/prob/prob.pdf