Divide and conquer!

I think the main point for performance is to have smallest objects as soon as possible.

For that, the idea is to use a dichotomy: split the object in half, then each half in half, etc. So cutting an object is faster at each iteration.

The following code is doing it using the bounding box of the object. Boundings are in local coordinates in Blender.

So we get the boundings. Find the spatially largest axis then cut in half on it.

Once vertices are selected, we use 'select linked' to avoid the faces or edges to be cut.

import bpy
import bmesh
import time
from mathutils import Vector


def cut_object( obj ):
    #Gets the bounds
    bounds = [b[:] for b in obj.bound_box]
    #0 is the min
    min_bounds = Vector( bounds[0] )
    #6 is opposite corner to 0
    max_bounds = Vector( bounds[6] )
    delta_bounds = max_bounds - min_bounds
    
    #Get the axis which corresponds to the largest area
    axis = 0
    if ( delta_bounds.x < delta_bounds.y ): axis = 1
    #Commented below as using Z axis has unwanted effects in your situation
    #if ( delta_bounds.y < delta_bounds.z ): axis = 2

    #The center is the limit    
    limit = min_bounds[axis] + 0.5 * delta_bounds[axis]

    #Set the object active
    bpy.ops.object.mode_set(mode = 'OBJECT')
    bpy.context.scene.objects.active = obj
    #Set selection mode
    bpy.ops.object.mode_set(mode="EDIT")
    bpy.ops.mesh.select_mode(type="VERT")
    
    #Select the wanted vertices
    bpy.ops.object.mode_set(mode = 'OBJECT')
    
    for v in obj.data.vertices:
        v.select = v.co[axis] < limit

    #Extend to linked parts (so wont cut the faces or edges)    
    bpy.ops.object.mode_set(mode="EDIT")
    bpy.ops.mesh.select_linked()
    
    #Get the selected vertex count
    selected_count = sum( 1 for v in obj.data.vertices if v.select )
    
    print( "Axis: " + str( axis ) )
    print( "Selected: " + str( selected_count ) + "/" + str( len( obj.data.vertices ) ) )
    
    result = True
    #Check that will separate something
    if selected_count == 0 or selected_count == len( obj.data.vertices ):
        result = False
    else:
        bpy.ops.mesh.separate(type='SELECTED')
    
    bpy.ops.object.mode_set(mode="OBJECT")
    
    return result
    
print( '-----------------' )

max_vertices = 10000
max_iterations = 20
iteration = 0

found = True

start_time = time.time()

discarded_objects = {}

while found and iteration < max_iterations:
    
    found = False

    #Get all objects that have more than the wanted vertex amount
    for obj in [x for x in bpy.context.scene.objects if len(x.data.vertices) > max_vertices and discarded_objects.get( x ) != True]:
        if cut_object( obj ):
            found = True
        else:
            #We can not cut it, so discard it
            discarded_objects[obj] = True

    iteration += 1            
    print( "Iteration: " + str( iteration ) + "   Objects: " + str( len( bpy.context.scene.objects ) ) )
    
end_time = time.time()

print( "Elapse: " + str( end_time - start_time ) )

The code is commented. So I think all is said about that.

Results

For a limit of 25k per chunk that gives 33s here.

For 10k, this is 50s.

Comments

You may prefer cutting differently. For instance defining a virtual cube surrounding the object, then cut it in half and use its part to cut the mesh.

The given code version does not included cutting along Z. It works, but is far longer when the wanted amount of vertices is low (10k for instance). The reason is we can only separate few faces at each iteration.

I'm also wondering why you need to start with such large mesh? Why the process does not allow to start from small parts instead this big one?

Here is the blend file with the script included: