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};Omegaright}$, $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$.