We present an Lq(Lp)-theory for the equation∂tαu=ϕ(Δ)u+f,t>0,x∈Rd;u(0,⋅)=u0. Here p,q>1, α∈(0,1), ∂tα is the Caputo fractional derivative of order α, and ϕ is a Bernstein function satisfying the following: ∃δ0∈(0,1] and c>0 such that(0.1)c(Rr)δ0≤ϕ(R)ϕ(r),0<r<R<∞. We prove uniqueness and existence results in Sobolev spaces, and obtain maximal regularity results of the solution. In particular, we prove‖|∂tαu|+|u|+|ϕ(Δ)u|‖Lq([0,T];Lp)≤N(‖f‖Lq([0,T];Lp)+‖u0‖Bp,qϕ,2−2/αq), where Bp,qϕ,2−2/αq is a modified Besov space on Rd related to ϕ.Our approach is based on BMO estimate for p=q and vector-valued Calderón-Zygmund theorem for p≠q. The Littlewood-Paley theory is also used to treat the non-zero initial data problem. Our proofs rely on the derivative estimates of the fundamental solution, which are obtained in this article based on the probability theory.