We investigate the combination of Isogeometric Analysis (IGA) and proper orthogonal decomposition (POD) based on the Galerkin method for model order reduction of linear parabolic partial differential equations. For the proposed fully discrete scheme, the associated numerical error features three components due to spatial discretization by IGA, time discretization with the -scheme, and eigenvalue truncation by POD. First, we prove a priori error estimates of the spatial IGA semi-discrete scheme. Then, we show stability and prove a priori error estimates of the space-time discrete scheme and the fully discrete IGA--POD Galerkin scheme. Numerical tests are provided to show efficiency and accuracy of NURBS-based IGA for model order reduction in comparison with standard finite element-based POD techniques.