We use online data assimilation to combine information from a linear inverse model of coupled atmosphere-ocean dynamics with proxy records to create a new annual-resolution reconstruction of atmosphere and ocean fields over the last millennium. Instrumental validation of reconstructed sea-surface temperature and 0-700 m ocean heat content shows broad regions of positive spatial correlations, and high correlations (~0.6-0.9) for global averages and indices of large-scale modes of atmospheric variability. Compared to previous reconstructions, the online reconstructions show global and hemispheric averages with little-to-no millennial-scale trend and global-mean temperatures ~0.25-0.5 K cooler during early periods (1000-1400 C.E.). The spatial anomaly differences of average temperature between an early (1000-1250 C.E.) and later (1400-1700 C.E.) period show warm anomalies over high-latitude Europe and cool tropical conditions in partial agreement with previous assessments. The addition of online data assimilation, which provides dynamical memory to climate proxy information, is shown to be crucial for adequately characterizing decadal-to-centennial-scale variability of 0–700 m ocean heat content. Furthermore, the climate forecasts provide model-based physical constraints for atmosphere-ocean interaction, which become increasingly important during early periods when less proxy information is available for assimilation.