We consider harmonic Toeplitz operators T-V = PV :H(Omega) -> H(Omega) where P : L-2(Omega) -> H(Omega) is the orthogonal projection onto H(Omega) = {u is an element of L-2 (Omega)) vertical bar Delta u = 0 in Omega}, Omega subset of R-d, d >= 2, is a bounded domain with boundary partial derivative Omega is an element of C-infinity and V : Omega -> C is an appropriate multiplier. First, we complement the known criteria which guarantee that T-V is in the pth Schatten-von Neumann class S-p, by simple sufficient conditions which imply T-V is an element of S-p(,w), the weak counterpart of S-p. Next, we consider symbols V >= 0 which have a regular power-like decay of rate & nbsp;gamma > 0 at partial derivative Omega, and we show that T-V is unitarily equivalent to a classical pseudo-differential operator of order-gamma, self-adjoint in L-2 (partial derivative Omega). Utilizing this unitary equivalence, we obtain the main asymptotic term of the eigenvalue counting function for T-V, and establish a sharp remainder estimate. Further, we assume that Omega is the unit ball in R-d, and V = (V) over bar is compactly supported in Omega, and investigate the eigenvalue asymptotics of the Toeplitz operator T-V. Finally, we introduce the Krein Laplacian K, self-adjoint in L-2 (Omega), perturb it by a multiplier V is an element of C((Omega) over bar; R), and show that sigma(ess)(K + V) = V (partial derivative Omega). Assuming that V >= 0 and V-vertical bar partial derivative Omega = 0, we study the asymptotic distribution of the discrete spectrum of K +/- V near the origin, and find that the effective Hamiltonian which governs this distribution is the Toeplitz operator T-V.