The recent gravitational wave transient GW190521 has been interpreted by the LIGO-Virgo collaboration (LVC) as sourced by a binary black hole (BH) merger. According to the LVC parameter estimation, at least one of these progenitors falls into the so-called pair-instability supernova mass gap. This raises the important question of how and when these progenitors formed. In this paper we use an accretion model with super-Eddington mass accretion rate obtained from General Relativity hydrodynamics simulations to analyse the scenario wherein the GW190521 original progenitors (OPs) formed at lower masses (and spins) and grew to their estimated LVC parameters by relativistic accretion. We consider that the environment wherein the binary is immersed has density gradients as well as a dependence on the Mach number of the gas. Taking the LVC parameter estimation at $z=0.82$ as the endpoint of the accretion evolution, we estimate the initial masses and spins of the OPs at three different red-shifts $z=100, 50$, and $20$. We found three distinct possible types of OPs: $(i)$ $10^{-4} M_{odot} - 3 M_{odot}$ almost non-rotating (with Kerr spin parameter $a_{star}< 10^{-2}$) primordial BHs; $(ii)$ $3 M_{odot} - 40M_{odot}$ slowly rotating ($ 10^{-2} < a_{star} < 0.5$) stellar mass BHs; $(iii)$ $40M_{odot} - 70M_{odot}$ BHs with a moderate spin parameter $a_{star}sim 0.5$, which could originate from the collapse of high mass Pop III stars. The mass spread is due to varying the density gradient and the relativistic Mach number of the cosmic plasma; the variation of the masses due to the origin at different red-shifts, on the other hand, is negligible, $sim 2%$ ...