The origins of the Galilean satellites – namely Io, Europa, Ganymede, and Callisto – is not fully understood yet. Here we use N-body numerical simulations to study the formation of Galilean satellites in a gaseous circumplanetary disc around Jupiter. Our model includes the effects of pebble accretion, gas-driven migration, and gas tidal damping and drag. Satellitesimals in our simulations first grow via pebble accretion and start to migrate inwards. When they reach the trap at the disc inner edge, scattering events and collisions take place promoting additional growth. Growing satellites eventually reach a multi-resonant configuration anchored at the disc inner edge. Our results show that an integrated pebble flux of ≥2 × 10−3 MJ results in the formation of satellites with masses typically larger than those of the Galilean satellites. Our best match to the masses of the Galilean satellites is produced in simulations where the integrated pebble flux is ∼10−3 MJ. These simulations typically produce between three and five satellites. In our best analogues, adjacent satellite pairs are all locked in 2:1 mean motion resonances. However, they have also moderately eccentric orbits (∼0.1), unlike the current real satellites. We propose that the Galilean satellites system is a primordial resonant chain, similar to exoplanet systems as TRAPPIST-1, Kepler-223, and TOI-178. Callisto was probably in resonance with Ganymede in the past but left this configuration – without breaking the Laplacian resonance – via divergent migration due to tidal planet–satellite interactions. These same effects further damped the orbital eccentricities of these satellites down to their current values (∼0.001). Our results support the hypothesis that Io and Europa were born with water-ice rich compositions and lost all/most of their water afterwards. Firmer constraints on the primordial compositions of the Galilean satellites are crucial to distinguish formation models.