The adaptive evolution of the SARS-CoV-2 spike glycoprotein represents a dynamic trajectory through high-dimensional sequence space, where viral fitness is constrained by epistatic couplings across structural domains. While single-site substitutions are well-characterized, the collective dynamics and long-range interactions between non-adjacent amino acid residues remain difficult to map. Analyzing time-resolved evolutionary trajectories from over 15 million viral sequences, we introduce a punctuated tensor decomposed pair clustering method (PunTeD PCM) to extract nonlocal residue couplings. We identify connected co-mutation networks spanning disparate functional domains and show that residue coevolution exhibits punctuated, non-equilibrium dynamics synchronized with major variant sweeps. Mapping these epistatic networks onto adaptive sequence landscapes quantifies how long-range epistatic interactions shape viral fitness, establishing a statistical-physics framework for modeling evolutionary dynamics in complex biological systems.