To embrace big-data neuroimaging, harmonization of site effect in resting-state functional magnetic resonance imaging (R-fMRI) data fusion is a fundamental challenge. Comprehensive evaluation of potentially effective harmonization strategies, particularly with specifically collected data has been rare, especially for R-fMRI metrics. Here, we comprehensively assess harmonization strategies from multiple perspectives, including efficiency, individual identification, test-retest reliability and replicability of group-level statistical results, on widely used R-fMRI metrics across multiple datasets including data obtained from the same participants scanned at several sites. For individual identifiability (i.e., whether the same subject could be identified across R-fMRI data scanned across different sites), we found that, while most methods decreased site effects, the Subsampling Maximum-mean-distance based distribution shift correction Algorithm (SMA) outperformed linear regression models, linear mixed models, ComBat series and invariant conditional variational auto-encoder. Test-retest reliability was better for SMA and adjusted ComBat series than alternatives, while SMA was superior to the latter in replicability, both in terms of Dice coefficient and the scale of brain areas showing sex differences reproducibly observed across datasets. Moreover, we examined test-retest datasets to identify the best target site features to optimize SMA identifiability and test-retest reliability. We noted that both sample size and distribution of the target site matter and introduced a heuristic target site selection formula. In addition to providing practical guidelines, this work can inform continuing improvements and innovations in harmonizing methodologies for big R-fMRI data.