Firstly, the first derivative discontinuities you observe in your heightmap texture are just an artifact of flattening the spherical/angular domain into six gnomonic projection "charts". Because the projection changes abruptly at the edge of a cube, the gradients have an apparent change in direction there too. But once this is re-projected to a sphere, that kink vanishes, as the projection for reading the cubemap exactly cancels out the projection used for writing it. So as long as your ultimate intent is to render this texture on a sphere or otherwise in the angular domain (like a skybox), these projection artifacts are safe to ignore.

As for how to avoid problems when trying to convert this greyscale heightmap into normals, my advice is: don't.

Calculate your normals analytically as part of evaluating the noise function, rather than trying to infer them later by taking finite differences of the heightmap.

This is cheaper and gives higher-quality normals, and saves you from a world of edge- and corner-case grief in trying to get sample-based values to agree across all the tangent space seams.

You can find examples of how to compute derivatives of a gradient noise field simultaneously as you sample it in this previous answer (for Simplex noise) or in this article from Inigo Quilez.

This will give you the direction in 3D space in which the noise field is getting brighter fastest. Zero-out the radial component to get the rate of change along the sphere's surface at that point. The magnitude tells you the slope of the height function "If you slide one unit (run) in this direction, height value increases by rise = length(gradient)". Multiply that by your height scale to control the intensity of the normals.

Here's some example code showing how we can use these analytical gradients to compute a normal vector in world space:

vec3 sampleDir = normalize(worldPos - center);
vec4 noise = noised(sampleDir); // Using IQ's "noise-with-derivative"

// IQ's version returns noise value in x, gradient in yzw.
float height = noise.x;
vec3 gradient = noise.yzw;

// Zero out component perpendicular to the sphere.
float3 onSphere = gradient - dot(sampleDir, gradient) * sampleDir;

// This is a scalar measure of incline in the familiar "rise over run" sense.
float slope = length(onSphere) * heightScale;

float3 normal = sampleDir;

// If the slope is zero, "straight outward" is our normal. Otherwise...
if (slope > 0.0) {
    // Unit vector in sphere's tangent plane, pointing to 'uphill' direction.
    vec3 run = onSphere/slope;
    // Perpendicular vector as large as the change in height over that run.
    vec3 rise = slope * sampleDir;

    // Vector pointing uphill along the hill.
    vec3 alongSlope = rise + run;

    // Bitangent vector - perpendicular to sphere normal and hill gradient.
    float3 perpendicular = cross(onSphere, sampleDir);

    // We use this bitangent to rotate the vector up the slope 90 degrees
    // to give us our world space normal.
    normal = normalize(cross(perpendicular, alongSlope));
}

The advantage of forming this normal in world space is that all adjacent cube faces will agree on the value of the normal where they meet. The downside is that you lose some precision when storing object/world-space normals, you can't compress the texture as easily, and it's harder to blend the result with detail normal maps or scale the intensity of the normals after the fact.

If you prefer to work in tangent space, you can use the TBN matrix at that point of the cube-sphere to transform the world space result into the tangent space for the plane. As long as you use the same tangent space for writing and for reading, you'll get back to the same "net" normal after the fact, and you should still have seamless agreement between adjacent faces. But if you store the normal map in tangent space, don't store/read it as a cubemap (adjacent faces don't always agree about the orientation of the tangent space, causing ugly interpolation where they meet).