Exporting layer tree in QGIS

Question

I need to create a layer tree for other user researchers to use. In the final dataset there will be ~20 sub-groups, each with dozens of layers. (It's possible that this isn't the best way to do it -- the reason for so many layers is to have individually-selectable lines, where the end user can hide all layers but one.)

I want to export the top level group in a way that:

• preserves the organization of layers into sub-groups and the top-level group
• includes the symbology
• is compatible with being loaded into existing QGIS projects.

What I’ve tried so far:

• Programmatically saving a layer: I can use QgsVectorFileWriter.writeAsVectorFormat to save a single layer, but I haven’t figured out how to get multiple layers into the same file or to save the layer groups.
• Using the Processing Toolbox’s Database -> Package layers function. This makes it straightforward to export as many layers as desired, but I can’t figure out how to preserve the hierarchy.
• [23 Nov update] Right-click on top level group -> Export -> Save as Layer Definition File. This saves the hierarchy, but when imported, none of the layers have any data. If instead I use layers with a CSV as the provider, this kind of works, but the data isn't packaged with the exported file -- instead it is a relative path on my computer, which isn't good for distribution.

Is there a way to do this?

I'm using QGIS 3.22 on a Mac.

---

23 Nov Update:

Here's a screenshot of my demo QGIS project; I'm trying to export the tree rooted at "Project Index":

And the code used to generate that layer tree:

import random

root = QgsProject.instance().layerTreeRoot()
top_level_group = root.insertGroup(0, "Project Index")  # Inserts at specified index

institutions = {1: ['A', 'B', 'C'], 2: ['D', 'E'], 3: ['F']}

for institution, projects in institutions.items():
    institution_group = top_level_group.addGroup("Institution {}".format(institution))
    for project in projects:
        project_group = institution_group.addGroup("Project {}".format(project))
        # Add a point layer containing a single point
        pt_x = random.randint(-10000, 10000)  # 20km wide box about origin
        pt_y = random.randint(-10000, 10000)
        pt_uri = "point?crs=epsg:3031"
        pt_layer = QgsVectorLayer(pt_uri, "Point", "memory")
        pt_feature = QgsFeature()
        pt_feature.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(pt_x, pt_y)))
        pt_prov = pt_layer.dataProvider()
        pt_prov.addFeature(pt_feature)
        pt_layer.updateExtents()
        QgsProject.instance().addMapLayer(pt_layer, False) 
        project_group.addLayer(pt_layer)
        
        # Add a polygon layer containing a 5km square
        poly_x0 = random.randint(-10000, 10000)
        poly_y0 = random.randint(-10000, 10000)
        vertices = [(poly_x0, poly_y0), (poly_x0, poly_y0 + 5000),
                    (poly_x0+5000, poly_y0+5000), (poly_x0+5000, poly_y0)]
        poly_points = [QgsPointXY(xx, yy) for xx,yy in vertices]
        poly_uri = "polygon?crs=epsg:3031"
        poly_layer = QgsVectorLayer(poly_uri, "Polygon", "memory")
        poly_feature = QgsFeature()
        poly_feature.setGeometry(QgsGeometry.fromPolygonXY([poly_points]))
        poly_prov = poly_layer.dataProvider()
        poly_prov.addFeature(poly_feature)
        poly_layer.updateExtents()
        QgsProject.instance().addMapLayer(poly_layer, False) 
        project_group.addLayer(poly_layer)
        
        # Add a polyline layer
        npts = random.randint(5,10)
        line_points = []
        for _ in range(npts):
            xx = random.randint(-10000, 10000)
            yy = random.randint(-10000, 10000)
            line_points.append(QgsPoint(xx, yy))
        print(line_points)
        line_uri = "linestring?crs=epsg:3031"
        line_layer = QgsVectorLayer(line_uri, "LineString", "memory")
        line_feature = QgsFeature()
        line_feature.setGeometry(QgsGeometry.fromPolyline(line_points))
        line_prov = line_layer.dataProvider()
        line_prov.addFeature(line_feature)
        line_layer.updateExtents()
        QgsProject.instance().addMapLayer(line_layer, False) 
        project_group.addLayer(line_layer)