That is a lot of overlapping questions to ask, and that makes it difficult to answer coherently. Also, since you do not provide the commands you used for analysis, I have to guess at what you mean.

"Can the cluster size threshold output of clustsim then apply throughout the entire brain (e.g., even in those slices that don't have data for all subjects)?"

If you have such gaps in coverage (very common), then you can do either of 2 things:

a) use 3dttest++ with '-zskip' to tell it to skip value that are exactly 0 - for this to work you will have to mask each subject's data with the proper coverage mask -- afni_proc.py produces such a mask, named mask_epi_anat.SUBJECTID+tlrc -- but note that afni_proc.py does NOT apply this mask to the stats results file. If you want this approach, you will have to mask the stats dataset with a 3dcalc command.

b) use a mask in 3dttest++ that is the intersection of all the individual coverage masks, so that you actually do have complete data (from each subject) at each analyzed voxel -- you can create such a mask by using the 3dmask_tool program and applying it to all the mask_epi_anat.SUBJECTID+tlrc.HEAD datasets

From what you say, it seems like you are using choice "a" above. I believe that it is proper to use the cluster thresholds across the entire masked dataset -- when individual subjects are properly masked, as described above, and when a proper mask is supplied to 3dttest++. But then you have to decide what the proper mask to use in 3dttest++ (for ClustSim) will be. If you are using different masks for each subject, then logically the -mask option to 3dttest++ should be the union of the individual masks instead of the intersection. To be clear, when I saw "across the entire masked dataset", I mean that the cluster threshold applies when used within the mask supplied to 3dttest++. If you were to use the cluster thresholds across a larger volume, that would be incorrect.

My personal preference is to use choice "b" above and make the intersection mask of all subjects -- and thus analyze only voxels that are common to all subjects. If this isn't acceptable, then you can use "a" but need to be careful with your masking, as I described above.

"That is, if I want to visualize in a different program (I know, GASP!), would I need to threshold on z-stats (selecting as my threshold the z-stat that corresponds to the p-value that I want to use as my cluster defining threshold, e.g., p<.001) in order for me to use the cluster size threshold output from 3dclustsim to validly control for multiple comparisons? Or would it be sufficient to threshold on t-stats?"

I don't completely understand what your suggestion is in the above question. It sounds to me like you want to use "3dttest++ -clustsim" to produce the cluster threshold. 3dttest++ will produce a dataset with the threshold volume being made up of z-statistics. So to threshold on t-statistics, you must want to not only visualize in another GASP program, but also re-analyze the data in such a GASP program since you won't have t-statistics from 3dttest++. I cannot advise you on re-analyzing the data in another program.

"Second: Is this approach different from randomise in FSL? Or is the main difference that the underlying stat is a z-stat instead of a t?"

I do not pay much attention to how either SPM or FSL work, and it would not be proper for me to give advise about using these other software packages.