ABSTRACT We study the gravitational fragmentation of circumstellar discs accreting extremely metal-poor ($Z \le 10^{-3}\, \mathrm{Z}_{\odot }$) gas, performing a suite of 3D hydrodynamic simulations using the adaptive mesh refinement code enzo. We systematically follow the long-term evolution for 2 × 103 yr after the first protostar’s birth, for the cases of Z = 0, 10−5, 10−4, and $10^{-3}\, \mathrm{Z}_{\odot }$. We show that evolution of number of self-gravitating clumps qualitatively changes with Z. Vigorous fragmentation induced by dust cooling occurs in the metal-poor cases, temporarily providing ∼10 self-gravitating clumps at Z = 10−5 and $10^{-4}\, \mathrm{Z}_{\odot }$. However, we also show that the fragmentation is a very sporadic process; after an early episode of the fragmentation, the number of clumps continuously decreases as they merge away in these cases. The vigorous fragmentation tends to occur later with the higher Z, reflecting that the dust-induced fragmentation is most efficient at the lower density. At $Z = 10^{-3}\, \mathrm{Z}_{\odot }$, as a result, the clump number stays smallest until the disc fragmentation starts in a late stage. We also show that the clump mass distribution depends on the metallicity. A single or binary clump substantially more massive than the others appear only at $Z = 10^{-3}\, \mathrm{Z}_{\odot }$, whereas they are more evenly distributed in mass at the lower metallicities. We suggest that the disc fragmentation should provide the stellar multiple systems, but their properties drastically change with a tiny amount of metals.