I tried to solve the problem in 2 ways:

1) Using CombineGVCFs which produced a combined vcf file:

samples=$(find . | sed 's/.\///' | grep -E 'g.vcf$' | sed 's/^/--variant /') # place sample paths into variable
path/to/gatk --java-options "-Xmx36G" CombineGVCFs \
       $(echo $samples) \
       -O path/to/combined.vcf \
       -R path/to/ref.fa

2) Using GenomicsDBImport (after making bed file for my reference) which produced a my_database folder:

samtools faidx ref.fna # make .fai file 
awk '{print $1 "\t0\t" $2}' ref.fna.fai > ref.fna.bed # make .bed file

samples=$(find . | sed 's/.\///' | grep -E 'g.vcf$' | sed 's/^/--variant /') # place sample paths into variable
path/to/gatk --java-options "-Xmx36G" GenomicsDBImport \
          $(echo $samples)\
          --genomicsdb-workspace-path my_database \
          --intervals path/to/ref.fna.bed

After that, I run GenotypeGVCFs using "my_database" or "combined.vcf"

 path/to/gatk --java-options "-Xmx36G" GenotypeGVCFs \
           -R path/to/ref.fna \
           -V gendb://path/to/my_database OR -V path/to/combined.vcf \
           -O path/to/genotypes.vcf

In both cases, I got the same number of SNPs in genotypes.vcf.