Land surface temperature (LST) is a crucial parameter driving the dynamics of the thermal state on land surface. In high-latitude cold region, a long-term, stable LST product is of great importance in examining the distribution and degradation of permafrost under pressure of global warming. In this study, a generalized additive model (GAM) approach was developed to fill the missing pixels of the MODIS/Terra 8-day Land Surface Temperature (MODIS LST) daytime products with the ERA5 Land Skin Temperature (ERA5ST) dataset in a high-latitude watershed in Eurasia. Comparison at valid pixels revealed that the MODIS LST was 4.8–13.0 °C higher than ERA5ST, which varies with land covers and seasons. The GAM models fairly explained the LST differences between the two products from multiple covariates including satellite-extracted environmental variables (i.e., normalized difference water index (NDWI), normalized difference vegetation index (NDVI), and normalized difference snow index (NDSI) as well as locational information. Considering the dramatic seasonal variation of vegetation and frequent snow in the cold region, the gap-filling was conducted in two seasons. The results revealed the root mean square errors (RMSE) of 2.7 °C and 3.4 °C between the valid MODIS LST and GAM-simulated LST data in the growing season and snowing season, respectively. By including the satellite-extracted land surface information in the GAM model, localized variations of land surface temperature that are often lost in the reanalysis data were effectively compensated. Specifically, land surface wetness (NDWI) was found to be the greatest contributor to explaining the differences between the two products. Vegetation (NDVI) was useful in the growing season and snow cover (NDSI) cannot be ignored in the snow season of the study region. The km-scale gap-filled MODIS LST products provide spatially and temporally continuous details that are useful for monitoring permafrost degradation in cold regions in scenarios of global climate change.