Deriving evolution equations accounting for both anomalous diffusion and reactions is notoriously difficult, even in the simplest cases. In contrast to normal diffusion, reaction kinetics cannot be incorporated into evolution equations modeling subdiffusion by merely adding reaction terms to the equations describing spatial movement. A series of previous works derived fractional reaction-diffusion equations for the spatiotemporal evolution of particles undergoing subdiffusion in one space dimension with linear reactions between a finite number of discrete states. In this paper, we first give a short and elementary proof of these previous results. We then show how this argument gives the evolution equations for more general cases, including subdiffusion following any fractional Fokker-Planck equation in an arbitrary $d$-dimensional spatial domain with time-dependent reactions between infinitely many discrete states. In contrast to previous works which employed a variety of technical mathematical methods, our analysis reveals that the evolution equations follow from (i) the probabilistic independence of the stochastic spatial and discrete processes describing a single particle and (ii) the linearity of the integro-differential operators describing spatial movement. We also apply our results to systems combining reactions with superdiffusion.