Answer. Yes, for appropriate boundary conditions (e.g., Dirichlet or Neumann) the Laplace operator on bounded domains with sufficiently smooth boundary has no positive eigenfuntions, except for those that belong to the leading eigenvalue.
Here are the details:
Part 0. Notation.
Consider an $L^p$-space over a $\sigma$-finite measure space. For an element of $f$ of an $L^p$-space over of $\sigma$-finite measure space, where $p \in [1,\infty)$, I'll use the notation $f \ge 0$ to mean $f(\omega) \ge 0$ for almost all $\omega$. I'll use the notation $f \gg 0$ to mean $f(\omega) > 0$ for almost all $\omega$ (some people use the notation $f > 0$ for that, but from an order theoretic point of view this notation is a bit unfortunate).
Part 1. A Perron-Frobenius (or Krein-Rutman) type result for operators which "improve" the positivity of functions:
Theorem. Let $(\Omega,\mu)$ be a $\sigma$-finite measure space and let $p \in [1,\infty)$ and assume that $L^p(\Omega,\mu)$ is non-zero. Let $T: L^p(\Omega,\mu) \to L^p(\Omega,\mu)$ be a compact linear operator with the following property: for non-zero $0 \le f \in L^p(\Omega,\mu)$ one has $Tf \gg 0$.
Then the spectral radius os $T$ satisfies $r(T) > 0$, it is an eigenvalue with an eigenvector $u \gg 0$, and it is the only eigenvalue of $T$ which has an eigenvector $v \ge 0$.
Proof. That fact that $r(T) > 0$ is a special case of a theorem of Ben de Pagter; see the main result of de Pagter's 1986 paper "Irreducible compact operators".
The fact that $r(T)$ is an eigenvalue with an eigenvector $u \ge 0$ is the classical Krein-Rutman theorem. The fact that $u$ actually satisfies $u \gg 0$ is an immediate consequence of the positivity improving assumption on $T$. The fact that $r(T)$ is the only eigenvalue of $T$ with an eigenvector $v \ge 0$ can, for instance, be found in Theorem V.5.2(iv) and its corollary (both on page 329) of Helmut H. Schaefer's 1974 book "Banach lattices and positive operators". $\square$
Part 2. Application to the Laplace operator with local boundary conditions.
Consider a bounded domain $\Omega$ in $\mathbb{R}^d$ and assume that the boundary of $\Omega$ is sufficiently smooth to ensure that the Laplace operator with the boundary conditions that we consider in the following has compact resolvent.
Let $\Delta: L^2(\Omega) \supseteq \operatorname{dom}(\Delta) \to L^2(\Omega)$ denote the Laplace operator, where the domain $\operatorname{dom}(\Delta)$ encodes the boundary conditions.
We assume the boundary conditions to be one of the classical choices, for instance Dirichlet or Neumann or mixed boundary conditions. Then we have the following result:
Corollary. If $\lambda$ is an eigenvalue of $\Delta$ with an eigenfunction $u \ge 0$, then $\lambda$ is the leading eigenvalue of $\Delta$.
Proof.
Let $\mu \in \mathbb{R}$ be a number that is strictly larger than the largest eigenvalue of $\Delta$.
Then the resolvent operator $(\mu I - \Delta)^{-1}: L^2(\Omega) \to L^2(\Omega)$ satisfies the assumptions of the operator $T$ in the theorem in part 1. Moreover, the spectral mapping theorem for the resolvent tells us that the mapping $z \mapsto \frac{1}{\mu - z}$ maps the eigenvalues of $\Delta$ $1$-to-$1$ to the eigenvalues of $(\mu I - \Delta)^{-1}$ and that the eigenvectors are thereby preserved. In particular, $\frac{1}{\mu -\lambda_0}$ is the spectral radius of $(\mu I - \Delta)^{-1}$.
So by the theorem in part 1, the claim follows. $\square$
Part 3. Further remarks.
The same arguments apply to elliptic operators with coefficients and with lower order terms.
For self-adjoint operators with compact resolvent (like the Laplace operator with Dirchlet or Neumann or mixed boundary conditions), there is a simpler argument that was pointed out in Christian Remling' answer.
An advantage of the argument sketched above is that it also works for non-self adjoint operators which one gets, for instance, by considering second-order operators with lower order terms.
The same result is no longer true for the Laplace operator with non-local boundary conditions. For instance, consider the one-dimenional domain $\Omega = (0,1)$ and the Laplace operator with the boundary conditions $u'(1) = \mathrm{e} u'(0)$. Then the constant function $1$ is an eigenvector for the eigenvalue $0$. But $1$ is also an eigenvalue (with $\exp$ as an eigenfunction).
The compactness assumption in the theorem in part 1 can, in the case $p=2$, be replaced with a self-adjointness assumption; see Theorem 2.2 in the 2015 paper "Note on basic features of large time behaviour of heat kernels" by Keller, Lenz, Vogt, and Wojciechowski.
(One can apply this reference to $-T$ to obtain the above theorem for self-adjoint but non-compact $T$. However, that's actually a bit of a detour, since the result in this reference is itself designed to deal with unbounded operators.)