Apply -1 x scale (mirror) to object mesh and keep custom normals with python
Alright, I found the solution with some help of the answer of batFINGER here.
There's surprisingly little information on the internet about manipulating face normals at the vertex positions, confusingly called 'split normals' or 'loop normals'.
To work with split normals with python you first need to use the method calc_normals_split() on the mesh or else you won't get accurate vector values for the split normals or loop normals.
It took me hours finding that out.
In stead of setting x scale to -1 and applying transforms at the object level, you should edit the vertices, face normals and split normals at the mesh level.
code:

import bpy

## get active object and mesh
ob = bpy.context.object
me = ob.data
## Very important to include line below for working with custom normals!
me.calc_normals_split()
me.use_auto_smooth = True
## go through each polygon, reverse second to last loop,
## multiply x vector of loops with -1, store mirrored vector in
## mirror_split_normals for later use.
mirror_split_normals = []
for iter_poly in me.polygons:
    reverse_normals_indices = [iter_poly.loop_indices[0]] + [i for i in reversed(iter_poly.loop_indices[1:])]
    for i in reverse_normals_indices:
        mirror_split_normals.append((me.loops[i].normal[0] * -1, me.loops[i].normal[1], me.loops[i].normal[2]))

## mirror vertices in x by multiplying x by -1
for v in me.vertices:
    v.co = v.co[0] * -1, v.co[1], v.co[2]

## flip face normals and apply mirrored split normals
me.flip_normals()
#me.update()
me.normals_split_custom_set(mirror_split_normals)

Original mesh:

Mirrored mesh: