We consider harmonic Toeplitz operators $T_V = PV:{\mathcal H}(\Omega) \to {\mathcal H}(\Omega)$ where $P: L^2(\Omega) \to {\mathcal H}(\Omega)$ is the orthogonal projection onto ${\mathcal H}(\Omega) = \left\{u \in L^2(\Omega)\,|\,\Delta u = 0 \; \mbox{in}\;\Omega\right\}$, $\Omega \subset {\mathbb R}^d$, $d \geq 2$, is a bounded domain with $\partial \Omega \in C^\infty$, and $V: \Omega \to {\mathbb C}$ is a suitable multiplier. First, we complement the known criteria which guarantee that $T_V$ is in the $p$th Schatten-von Neumann class $S_p$, by sufficient conditions which imply $T_V \in S_{p, {\rm w}}$, the weak counterpart of $S_p$. Next, we assume that $\Omega$ is the unit ball in ${\mathbb R}^d$, and $V = \overline{V}$ is radially symmetric, and investigate the eigenvalue asymptotics of $T_V$ if $V$ has a power-like decay at $\partial \Omega$ or $V$ is compactly supported in $\Omega$. Further, we consider general $\Omega$ and $V \geq 0$ which is regular in $\Omega$, and admits a power-like decay of rate $\gamma > 0$ at $\partial \Omega$, and we show that in this case $T_V$ is unitarily equivalent to a pseudo-differential operator of order $-\gamma$, self-adjoint in $L^2(\partial \Omega)$. Using this unitary equivalence, we obtain the main asymptotic term of the eigenvalue counting function for the operator $T_V$. Finally, we introduce the Krein Laplacian $K \geq 0$, self-adjoint in $L^2(\Omega)$; it is known that ${\rm Ker}\,K = {\mathcal H}(\Omega)$, and the zero eigenvalue of $K$ is isolated. We perturb $K$ by $V \in C(\overline{\Omega};{\mathbb R})$, and show that $\sigma_{\rm ess}(K+V) = V(\partial \Omega)$. Assuming that $V \geq 0$ and $V{|\partial \Omega} = 0$, we study the asymptotic distribution of the eigenvalues of $K \pm V$ near the origin, and find that the effective Hamiltonian which governs this distribution is the Toeplitz operator $T_V$.