We give the upper and the lower estimates of heat kernels for Schrödinger operators $H=-Δ+V$, with nonnegative and locally bounded potentials $V$ in $\mathbb{R}^d$, $d \geq 1$. We observe a factorization: the contribution of the potential is described separately for each spatial variable, and the interplay between the spatial variables is seen only through the Gaussian kernel - optimal in the lower bound and nearly optimal in the upper bound. In some regimes we observe the exponential decay in time with the rate corresponding to the bottom of the spectrum of $H$. Our estimates identify in a fairly informative and uniform way the dependence of the potential $V$ and the dimension $d$. The upper estimate is more local; it applies to general potentials, including confining and decaying ones, even if they are non-radial, and mixtures. The lower bound is specialized to confining case, and the contribution of the potential is described in terms of its radial upper profile. Our results take the sharpest form for confining potentials that are comparable to radial monotone profiles with sufficiently regular growth - in this case they lead to two-sided qualitatively sharp estimates. In particular, we describe the large-time behaviour of nonintrinsically ultracontractive Schrödinger semigroups - this problem was open for a long time. The methods we use combine probabilistic techniques with analytic ideas. We propose a new effective approach which leads us to short and direct proofs.