This paper introduces a fractional heat equation, where the diffusion operator is the composition of the Bessel and Riesz potentials. Sharp bounds are obtained for the variance of the spatial and temporal increments of the solution. These bounds establish the degree of singularity of the sample paths of the solution. In the case of unbounded spatial domain, a solution is formulated in terms of the Fourier transform of its spatially and temporally homogeneous Green function. The spectral density of the resulting solution is then obtained explicitly. The result implies that the solution of the fractional heat equation may possess spatial long-range dependence asymptotically.