The subject of this work is the application of fully discrete Galerkin finite element methods to initial-boundary value problems for linear partial integro-differential equations of parabolic type. We investigate numerical schemes based on the Padé discretization with respect to time and associated with certain quadrature formulas to approximate the integral term. A preliminary error estimate is established, which contains a term related to the quadrature rule to be specified. In particular, we consider quadrature rules with sparse quadrature points so as to limit the storage requirements, without sacrificing the order of overall convergence. For the backward Euler scheme, the Crank-Nicolson scheme, and a third-order (1,2) Padé-type scheme, the specific quadrature rules analyzed are based on the rectangular, the trapezoidal, and Simpson’s rule. For all the schemes studied, optimal-order error estimates are obtained in the case that the solution of the problem is smooth enough. Since this is important for our error analysis, we also discuss the regularity of the exact solutions of our equations. High-order regularity results with respect to both space and time are given for the solution of problems with smooth enough data.