Many of your questions are answered elsewhere, so I added additional links below but to summarize in one place relative to your questions I offer the explanations below that together with the links should answer your questions:
Required Sampling Frequency:
If the signal is real, then yes the sampling frequency must be greater than twice the highest frequency for a signal bandwidth that extends down to DC. However if the signal is complex then the sampling frequency should be greater than the highest frequency for same case. This is because a complex signal has a unique spectrum over the Nyquist frequency range of $\pm f_s/2$ where $f_s$ is the sampling rate, while for a real signal the negative frequencies are the complex conjugate of the positive frequencies (identical magnitude, opposite phase).
Note that prior to sampling an analog anti-alias filter is required on the signal, otherwise energy if present in higher Nyquist zones (frequencies above the Nyquist frequency range defined above) will alias into the Nyquist frequency range being observed. (see this post which explains how aliasing occurs: Aliasing after downsampling). It is not feasible to implement a brick-wall filter, so additional margin to the sampling rate is usually added to allow for the characteristics of the anti-alias filter. For this you need to understand the properties of the signal signal of interest prior to sampling (as once sampled there is no way to distinguish an alias from a signal in the first Nyquist zone!). For example in your case, where you are interested in observing the signal out to 20 Harmonics of 60 Hz: The anti-alias filter should pass the 20th harmonics without distortion and then reject the higher harmonics so that you can minimize the sampling rate needed. You would then choose the sampling rate depending on your filter implementation such that no energy folds in to distort the 20th harmonic.
How Many Point FFT is Needed:
This has been answered several times, so see the links given below for further details and caveats, but the main point is that the frequency resolution of each "point" is 1/T Hz where T is the length of your signal in seconds. Further the bins in the DFT (the FFT is just an algorithm that computes the DFT) as typically returned extend from DC to one bin less that $f_s$. Due to aliasing the DC is equivalent to $n f_s$ where n is any integer, and likewise for all other frequencies in the DFT, so the DFT that extends frmo DC to $f_s$ equivalently extends from DC to one bin less than $f_s/2$ and then from $-f_s/2$ to one bin less than $f_s$ (see fftshift in Matlab). As detailed in the links given, due to the relationship of the sampling rate and the total length in time T of your signal, each bin will have a frequency resolution of $1/f_s$ when no windowing is used.
Note importantly that each bin of the DFT is the result of a Sinc filter, so if your sampling rate is not an integer multiple of 60 Hz, you will also have the effect of "scalloping loss" that will distort the measurement of each harmonic. (see the "Bank of Filters" explanation at this link: Intuition for sidelobes in FFT). By choosing a sampling rate that is an integer multiple of 60 Hz, each 60 Hz harmonic will be properly centered in each bin.
Suggestion for your specific application:
With relationship to choosing the sampling rate and number of bins for your specific application, if the 60 Hz and its harmonics are the dominant signal and you don't have other signal energy of significance then there is a simple solution that will work quite elegantly in accurately estimating each harmonic with minimum measurement distortion. As described above, choose a sampling rate that is an integer multiple of 60 Hz, high enough to allow for the anti-alias filter used, and then choose the number of points such that each bin passes a harmonic of 60 Hz and rejects all the other harmonics. This means the number of samples in the FFT is $f_s/60$ where $f_s$ is an integer multiple of 60 Hz > 2400 Hz (assuming your signal is real). For example 4800 Hz may be a good choice as you suggest in your options (allowing for a fairly relaxed anti-alias filter!), and in this case $N =4800/60= 80$ points. So in this case $f_s=4800 Hz$ and with 80 points each point is 60 Hz, so 0, 60, 120, 180 etc centered on each harmonic and perfectly nulling the other Harmonics. The actual magnitude measured for each point is the total area under the Sinc function frequency response for each bin, so if there is other energy present at non-harmonic frequencies, this would distort the magnitude measured in each bin. But if the harmonics are dominant, this would be an excellent and simple approach. This is somewhat similar to the objective in this link in removing harmonic signals (Different way to separate a particular frequency from a signal) as each bin in the DFT is a rotated moving average filter just as done in this referenced link when only one bin is of interest. Note from that link how every other harmonic is nulled!
Please see these other responses which will further help give insight:
Specific frequency resolution
FFT frequency resolution
Number of FFT points required for a specific frequency resolution for an oversampled signal
FFT and number of samples relations
What happens when N increases in N-point DFT