This paper establishes a Feynman-Kac formula to represent the solution to general time inhomogeneous stochastic parabolic partial differential equations driven by multiplicative fractional Gaussian noises in bounded domain where L_t is a second order uniformly elliptic operator whose coefficients can depend on time and generates a time inhomoegenous Markov process. The idea is to use the Aronson bounds of fundamental solution of the associated heat kernel and the techniques from Malliavin calculus. The newly obtained Feynman-Kac formula is then applied to establish the Holder regularity in the space and time variables. The dependence on time of the coefficients poses serious challenges and new results about the stochastic differential equations are discovered to face the challenge. An amazing application of the Feynman-Kac formula is about the matching upper and lower bounds for all moments of the solution, an critical tool for the intermittency. For the latter result we need first to establish new small ball like bounds for the diffusion associated with the parabolic differential operator in the bounded domain which is of interest on its own.