I’m so glad I’ve finally solved this! Here are some tips for students like me (just a GISer without any knowledge about plate reconstruction).
We usually use ArcGIS or QGIS to create features in shapefile format, and we can easily add a new attribute column named “PlateID”. At first, I was confused about this “ID” - I mistakenly thought it was just like the FID in ArcGIS. So I didn’t properly connect my target points to plates in the GPlates GPML file, which is why my points didn’t rotate with the time slider. As a complete beginner, I struggled to understand the relationship between each point and its corresponding plate.
Then I read another post in the forum where John mentioned (pyGplates API to rotate shapefiles - #2 by john.cannon) that we can use pip to install Pygplates and import points from a text file while assigning PlateID. So I generated a txt file from my shapefile and used Python 3.9 (note the version requirement - I also got stuck on this minor issue). Then I used Zahirovic_etal_2022_OptimisedMantleRef_and_NNRMantleRef.rot and 250-0_plate_boundaries.gpml to construct points of interest (POI) for the 50 Ma to present time period. (Actually, I’m not entirely sure if my use of these two files is correct or not.) With this setup, I can now export the 50 Ma positions to use in PaleoDEM for extracting paleo-topographic DEMs.
Here I want to thank John and the entire GPlates team once again. Your responses are always timely and detailed, especially the introduction and sample code of pyGPlates!