Soybean protein content is a valuable quantitative trait controlled by multiple genes. The epistatic interaction of these genes can increase protein content observably. In this study, we used the multifactor dimensionality reduction method and a soybean high-density genetic map including 5308 markers to identify stable loci controlling protein content in soybean across 23 environments. In total, 31 897 046 single nucleotide polymorphism (SNP) – protein interaction pairs were detected. Among these, 46 stable SNP interaction pair associations with soybean protein content were identified under multiple environments, with 2 and 44 SNP pairs stably detected across four and three environments, respectively. Hot spot regions for interaction pairs were detected on linkage groups Gm17, Gm06, and Gm03, consistent with previous quantitative trait locus mapping. The epistatic effects and contributions of the stable interaction pairs ranged from 0.0008 to 0.5483 and 0.0003 to 0.5126, respectively. Eight SNP epistatic interaction subnets were constructed. Ten candidate genes from these interaction subnets showed a relationship with seed protein storage or amino acid biosynthesis and metabolism. The results of this study provide insights into the genetic architecture of soybean protein content and can serve as a basis for marker-assisted selection breeding.