We consider a system of $d$ non-linear stochastic fractional heat equations in spatial dimension $1$ driven by multiplicative $d$-dimensional space-time white noise. We establish a sharp Gaussian-type upper bound on the two-point probability density function of $(u(s, y), u (t, x))$. From this result, we deduce optimal lower bounds on hitting probabilities of the process $\{u(t, x): (t, x) \in [0, \infty[ \times \mathbb{R}\}$ in the non-Gaussian case, in terms of Newtonian capacity, which is as sharp as that in the Gaussian case. This also improves the result in Dalang, Khoshnevisan and Nualart [\textit{Probab. Theory Related Fields} \textbf{144} (2009) 371--424] for systems of classical stochastic heat equations. We also establish upper bounds on hitting probabilities of the solution in terms of Hausdorff measure.