Let $H$ be a Hilbert space and $E$ a Banach space. We set up a theory of stochastic integration of ${\cal L}(H,E)$-valued functions with respect to $H$-cylindrical Liouville fractional Brownian motion with arbitrary Hurst parameter $0<\beta <1$. For $0<\beta <\frac 12$ we show that a function $\Phi \colon (0,T)\to {\cal L}(H,E)$ is stochastically integrable with respect to an $H$-cylindrical Liouville fractional Brownian motion if and only if it is stochastically integrable with respect to an $H$-cylindrical fractional Brownian motion. We apply our results to stochastic evolution equations $$ {\rm d}U(t) = AU(t) {\rm d}t + B {\rm d}W_H^\beta (t) $$ driven by an $H$-cylindrical Liouville fractional Brownian motion, and prove existence, uniqueness and space-time regularity of mild solutions under various assumptions on the Banach space $E$, the operators $A\colon \scr D(A)\to E$ and $B\colon H\to E$, and the Hurst parameter $\beta $. As an application it is shown that second-order parabolic SPDEs on bounded domains in $\mathbb R^d$, driven by space-time noise which is white in space and Liouville fractional in time, admit a mild solution if $\frac {1}{4}d<\beta <1$.